Cosmic-ray interactions with the atmosphere produce a flux of neutrinos in all directions with energies extending above the TeV scale1. The Earth is not a fully transparent medium for neutrinos with energies above a few TeV, as the neutrino–nucleon cross-section is large enough to make the absorption probability non-negligible2. Since absorption depends on energy and distance travelled, studying the distribution of the TeV atmospheric neutrinos passing through the Earth offers an opportunity to infer its density profile3,4,5,6,7. This has never been done, however, due to the lack of relevant data. Here we perform a neutrino-based tomography of the Earth using actual data—one-year of through-going muon atmospheric neutrino data collected by the IceCube telescope8. Using only weak interactions, in a way that is completely independent of gravitational measurements, we are able to determine the mass of the Earth and its core, its moment of inertia, and to establish that the core is denser than the mantle. Our results demonstrate the feasibility of this approach to study the Earth’s internal structure, which is complementary to traditional geophysics methods. Neutrino tomography could become more competitive as soon as more statistics is available, provided that the sources of systematic uncertainties are fully under control.

A reliable estimate of the density profile of the Earth is essential to solve a number of important problems in geophysics, such as the dynamics of the core and mantle, the mechanism of the geomagnetic dynamo or the bulk composition of the Earth9. Most of our knowledge about the internal structure of the Earth and the physical properties of its different layers comes from geophysics and, in particular, from seismological data. Moreover, information from geodesy, geomagnetic and geodynamical data, solid state theory and high-temperature/pressure experimental results is also used.

The determination of the density distribution of the Earth from the bulk sound velocity of seismic waves in combination with normal modes is a well-established method with statistical uncertainties in the lower mantle and the outer core at the percent level and below for 250–300 km resolving intervals, with larger errors as radial resolution increases10,11. The reconstruction of a three-dimensional profile is, however, a very demanding nonlinear inversion problem of different seismic data10,11,12. Moreover, as wave velocities also depend on composition, temperature, pressure and elastic properties, this necessarily introduces uncertainties in the density estimate. Most studies of Earth’s radial structure are based on empirical relations between seismic wave velocities and density, such as Birch’s law, which may fail at the higher densities of Earth’s core, and the Adams–Williamson equation13. A good understanding of the Earth’s interior, aiming at simultaneously determining the density variations and the origin of such waves in terms of temperature and composition variations, cannot be done from seismic velocity variations alone, and another, independent piece of information is needed. Therefore, a precise modelling of the different layer compositions which are crossed by seismic waves is required. Even though several million earthquakes occur in the Earth every year, only of the order of hundred of them have magnitudes larger than 6 (ref. 14). Most of them do not occur on the surface, and the origin of the wave must be inferred by comparing time delays from different seismographs. Finally, only a small fraction of the registered seismic waves cross the Earth’s core. For all these reasons, using other complementary and independent methods to infer the density profile of the Earth is important.

Neutrinos can be used to study the Earth’s interior in several ways. First, experiments such as KamLAND and Borexino are currently measuring the so-called geo-neutrino flux (that is, neutrinos produced by the decay of radioactive elements in the Earth’s interior15,16), which provides information that can be used to understand its composition. On the other hand, a good knowledge of neutrino propagation through the Earth may give relevant information about the Earth’s density profile. Neutrino propagation does, indeed, depend on the details of the matter structure between the source and the detector. For neutrinos with energies below 1 TeV, the matter profile affects the neutrino oscillation pattern, whereas for neutrinos with energies in the multi-TeV range, the neutrino flux observed at the detector depends on the number of nucleons along its path, as neutrinos can undergo inelastic scattering and become absorbed17. Indeed, the idea of performing absorption radiographies of the Earth with neutrinos dates back to more than four decades ago. To our knowledge, the first mention of this possibility was advanced in an unpublished CERN preprint in October 1973 by Placci and Zavattini18 and by Volkova and Zatsepin in a talk in 197419, considering man-made neutrinos. The idea of combining Earth’s neutrino radiographies (that is, performing a neutrino tomography) is based on studying the attenuation of neutrinos crossing the Earth from different angles with respect to the position of the detector. The column depth traversed by a neutrino that has passed through the entire Earth’s diameter is 11 kton cm2 (1.1 × 1010 cm water equivalent). For neutrinos with an energy of ~40 TeV, the absorption length in the Earth becomes comparable to its diameter, \(\left( {n{\kern 1pt} \sigma } \right)^{ - 1}\sim 2R_ \oplus\), where n is the average nucleon number density, σ the neutrino–nucleon total cross-section and R = 6,371 km the mean radius of the Earth. Therefore, for few-TeV neutrinos there is a non-negligible probability for the incoming neutrino flux to be suppressed, enσL < 1, where \(L = 2R_ \oplus {\kern 1pt} {\mathrm{cos}}\,\theta _{\rm z}\) is the path length in the Earth as a function of the zenith angle θz (Fig. 1a).

Atmospheric neutrinos offer a large range of baselines (from a few to thousands of kilometres) and energies (from MeV to tens of TeV), with an energy spectrum that falls as ~E−3.7. Therefore, they represent a suitable source for neutrino tomography. Although neutrino interactions are rare, with the operation of kilometre-cube detectors such as IceCube, a large event sample can be harvested. In this work we use the publicly available IceCube one-year up-going muon sample, collected during 2011–2012 and referred to as IC86 (IceCube 86-string configuration), which contains 20,145 muons detected over a live time of 343.7 days8 (a preliminary attempt using IceCube data with very limited event statistics was presented in 201220). These muons are produced by up-going neutrinos and antineutrinos which, after crossing the Earth, interact via charged-current processes in the bedrock or ice surrounding the detector. While propagating inside the detector at a speed higher than the speed of light in ice, these muons emit Cherenkov light, which is detected by the digital optical modules of the IceCube array. Details about the data sample and the modelling of the predicted event rate are provided in Methods.

The energy and zenith distributions of the IC86 sample are shown in Fig. 1b. Since the atmospheric neutrino spectrum is a steeply falling function of the energy and, for the lowest energies, the neutrino absorption length is much larger than the Earth’s diameter, most of the neutrinos in the sample do not undergo significant absorption. Therefore, the distribution of the full sample is very similar to the atmospheric neutrino distribution at the Earth’s surface, which is more peaked towards the horizon1. For higher energies, however, the observed event spectrum corresponding to up-going neutrinos with the longest trajectories through the Earth (\({\mathrm{cos}}\,\theta _{\rm z}^{{\mathrm{rec}}}\sim - 1\)) is suppressed with respect to the zenith-symmetric flux corresponding to down-going neutrinos that propagate only a few tens of kilometres without crossing the Earth \(\left( {{\mathrm{cos}}\,\theta _{\rm z}^{{\mathrm{rec}}}\sim 1} \right)\). The effect is more pronounced for neutrinos with higher energies and for those with longer propagation paths in the Earth, as they have a larger probability of interaction. Hence, by studying the zenith and energy distributions of the atmospheric neutrino flux and comparing them with the flux without attenuation, information on the Earth’s density profile can be extracted. However, all events are useful: the events with the lowest energies or more horizontal trajectories serve to fix the overall normalization and zenith distribution of the unattenuated atmospheric neutrino flux.

To illustrate how to remove the intrinsic zenith dependence on the atmospheric neutrino flux when comparing with the observed data, we depict the ratio of the observed number of events to the expected one in the case of no absorption, Ndata/Nno_att, as a function of the zenith angle. If all energies in the IC86 sample are considered (Fig. 2a), statistics is dominated by the low-energy events and the maximum observed suppression is at the 10% level or below. For events with energies above 5 TeV (Fig. 2b), however, the suppression in some of the most vertical angular bins \(\left( {{\mathrm{cos}}\,\theta _{\rm z}^{{\mathrm{rec}}} < - 0.6} \right)\) is up to 50%. For all energies, the suppression is larger for more vertical trajectories, which imply a longer propagation path. As an indication, we also show the expectations for the central value and the 1σ statistical error of this ratio using the one-dimensional Preliminary Reference Earth Model (PREM)21.

We have parametrized the Earth’s density with a one-dimensional five-layer profile with constant density in each of the layers (Fig. 1a). One of the edges is chosen at the core/mantle boundary and another one at the inner core/outer core boundary, so that we select three layers in the core (one for the inner core and two for the outer core) and two layers in the mantle. We have checked that, with this number of layers, current data are not yet sensitive to the particular profile within a given layer (see Supplementary Figs. 3,4 and Supplementary Table 1) and, therefore, there is no expected gain when using more layers or a more realistic density profile. We fit the average density of each of the layers, which is allowed to vary freely, and obtain our main result, a one-dimensional density profile of Earth measured by means of weak interactions (Fig. 3). With one-year statistics the uncertainties are large, but still compatible with results from geophysical methods within a 68% credible interval. Notice that these results are obtained from one-dimensional marginalized posterior probability distributions and correlations among all the parameters in the fit (five densities and four nuisance parameters) are not shown here. Therefore, they give a conservative representation of the allowed ranges for the density of individual layers. For the interested reader, technical details regarding the fit procedure are given in Methods.

From the results of the fit, we compute the mass of the Earth as weighted by neutrinos and obtain \(M_ \oplus ^\nu = \left( {6.0_{ - 1.3}^{ + 1.6}} \right) \times 10^{24}\) kg (Fig. 4a), to be compared to the most precise gravitational measurement to date22,23 of \(M_ \oplus ^{{\mathrm{grav}}} = \left( {5.9722 \pm 0.0006} \right) \times 10^{24}\,{\mathrm{kg}}\). Clearly, albeit within large uncertainties, both results are in very good agreement.

We can also estimate the mass of the Earth’s core, a parameter that may be useful (as soon as statistical errors decrease) as an input for geophysical measurements of the Earth’s density profile. The result for this quantity is \(M_{{\mathrm{core}}}^\nu = \left( {2.72_{ - 0.89}^{ + 0.97}} \right) \times 10^{24}\,{\mathrm{kg}}\), which is slightly larger than the result from geophysical density models that estimate the mass of the core to be ~33% of the total mass of the Earth (see Fig. 4b).

From our measurement of the one-dimensional density profile we can determine the Earth’s moment of inertia, for which we get \(I_ \oplus ^\nu = \left( {6.9 \pm 2.4} \right) \times 10^{37}\,{\mathrm{kg}}\,{\mathrm{m}}^{\mathrm{2}}\) (Fig. 4c), in agreement with the current (gravitationally inferred) measurement24 of the mean moment of inertia, \(I_ \oplus ^{{\mathrm{grav}}} = \left( {8.01736 \pm 0.00097} \right) \times 10^{37}\,{\mathrm{kg}}\,{\mathrm{m}}^2\). The smaller moment of inertia from neutrino data, as compared to gravitational measurements, implies a central value with a larger departure from homogeneity, as shown in Fig. 4c (even though they are fully compatible between each other due to the still large uncertainties).

Another piece of information regarding the Earth’s interior that we can extract from the currently available data is to detect that the core is denser than the mantle. Notice that, implicitly, this is a strong assessment in favour of a non-homogeneous Earth (something that was expected to be possible to prove at 3σ after ten years of IceCube data3 and seems to be already established at more than 2σ using IC86 alone). We determine the difference between the average density within the two layers we divide the mantle into, \(\bar \rho _{{\mathrm{mantle}}}\), and the average density within the three layers corresponding to the core, \(\bar \rho _{{\mathrm{core}}}\). The result for this difference, measured by weak interactions, is \(\left( {\bar \rho _{{\mathrm{core}}}^\nu - \bar \rho _{{\mathrm{mantle}}}^\nu } \right) = 13.1_{ - 6.3}^{ + 5.8}\,{\mathrm{g}}\,{\mathrm{cm}}^{-3}\) (Fig. 4d). From this result, a denser Earth’s mantle has a p-value of 0.011 for our default model of the atmospheric neutrino flux.

As a test of consistency and as a matter of accounting for further systematic uncertainties, all observables have also been computed for other atmospheric neutrino fluxes, as well as using different modelling of the inner structure of the Earth. In all cases, the results are compatible with the ones presented here (Supplementary Figs. 14 and Supplementary Table 1).

At high enough energies (few TeV), the passing of neutrinos through the Earth is sensitive to the number density of nucleons and, therefore, this test represents an effective counting of nucleons in the Earth. Unlike gravitational methods, the estimation of the Earth’s mass with neutrinos relies purely on weak interactions and on the nucleon masses. Conceptually, this is a completely different method from gravitational ones. We have shown that, using the publicly available data from the IceCube neutrino telescope, this method starts being feasible. Future data will greatly improve the measurements presented here (we remind the reader that more data already collected by IceCube in the same energy range are not yet publicly available in the format required to perform this analysis, but hopefully will be released soon), including data from the future KM3NeT detector in the Mediterranean Sea25. For this reason, we have also estimated the projected sensitivity with future IceCube data (Supplementary Figs. 5,6).

As a final comment, it is important to stress that a non-gravitational measurement of Earth’s mass, as is the one presented here, could also probe that all the matter that contributes to Earth’s gravitational field is baryonic matter (protons, neutrons and electrons). With current neutrino data, however, a small fraction in the form of (non weakly-interacting) dark matter, which would not attenuate the passage of neutrinos, cannot yet be fully excluded.

The 2011–2012 up-going muon sample of the IceCube detector (IC86)8 was selected because of its very good angular information (better than a degree for TeV energies) and its good statistics around TeV energies. These 20,145 muons, detected over 343.7 days, were produced by atmospheric muon neutrino and antineutrino interactions in the medium surrounding the detector. In turn, these neutrinos originated from decays of atmospheric pions and kaons (and with a contamination from other sources below 0.1%) in the Northern Hemisphere, and have all traversed the Earth. The sample covers a solid angle of 2π, making it particularly suitable for the kind of study performed here.

The energy of the muons in the IC86 sample lies between ~400 GeV and ~20 TeV and is reconstructed, based on energy losses along the track, with a resolution of \(\sigma _{{\mathrm{log}}\left( {E_\mu /{\mathrm{GeV}}} \right)}\sim 0.5\). Since the median opening angle with respect to the parent neutrino direction is \(0.7{\kern 1pt} (E_\nu /{\mathrm{TeV}})^{ - 0.7}\) degrees8, the muon direction is a very good proxy for the original neutrino direction. The muon zenith angle can be reconstructed with a resolution in cos θz between 0.005 and 0.015.

Models of atmospheric neutrino fluxes

The atmospheric neutrino flux is characterized in terms of the cosmic-ray primary spectrum entering the atmosphere and the hadronic-interaction model that controls the development of the shower that finally produces the flux of neutrinos. Several choices for the model of the atmospheric neutrino flux are currently compatible with all available data. In this Letter, our default choice for the primary cosmic-ray spectrum is the combined Honda–Gaisser primary cosmic-ray spectrum with the Gaisser–Hillas H3a correction (HG-GH-H3a)26 and our default hadronic-interaction model is the QGSJET-II-04 hadronic-interaction model27. Nevertheless, we also considered the Zatsepin–Sokolskaya (ZS) cosmic-ray spectrum28 and the SIBYLL2.3 hadronic-interaction model29, and combined them to obtain a set of four different models for the atmospheric neutrino fluxes. The results for all these cases are shown in Supplementary Figs. 1,2 and Supplementary Table 1. This overall systematic uncertainty results in shifts of the allowed range for the fitted and derived quantities of about 20–30%, well within the statistical uncertainties. In all cases we always account for the variation of four continuous nuisance parameters, which we describe below.

Nuisance parameters

To account for systematic uncertainties, we consider four of the continuous nuisance parameters described in the IC86 paper8, to which we refer the reader for further details. First, the overall flux normalization (N) is allowed to vary within a factor of two of the central value of each model, which is larger than current uncertainties30,31. However, the low-energy component of the observed neutrino spectrum is extremely effective in substantially reducing the normalization uncertainty. Second, the pion-to-kaon ratio (π/K) determines the relative contribution to the neutrino flux from pion or kaon decays. A smaller value of this parameter implies a harder atmospheric neutrino spectrum. We normalize it to one and use a 10% Gaussian prior. Third, the uncertainty on the spectral shape of the atmospheric neutrino spectrum, Δγ, is accounted for by a tilt in the energy spectrum, with a pivot energy close to the median of the neutrino energy distribution. We add it as a Gaussian prior with a 5% error. Fourth, the uncertainty in the efficiency of the digital optical modules (DOMeff) affects the determination of the reconstructed energy: a smaller efficiency implies a shift to lower energies. We allow this parameter to vary freely between 0.9 and 1.19, given that its central value is 0.99.

We have estimated the contribution of these nuisance parameters to the error budget of the four derived quantities presented in the main text and in Supplementary Table 1 (Earth’s mass, core mass, moment of inertia and the difference in average density of the core and mantle) by comparing our results with the outcome of a fit where the four nuisance parameters have been fixed to their corresponding best-fit values. We have found that these systematic errors contribute to approximately 30% of the error on the derived quantities, \(M_ \oplus ^\nu\), \(M_{{\mathrm{core}}}^\nu\), \(I_ \oplus ^\nu\) and \(\bar \rho _{{\mathrm{core}}}^\nu - \bar \rho _{{\mathrm{mantle}}}^\nu\).

Finally, note that the optical properties of the ice surrounding the detector also represent an important source of systematic uncertainties, which are comparable to current uncertainties32. This may limit the potential of the method once statistical uncertainties decrease. Nevertheless, with the available information we cannot further investigate this issue further. This should be done with results from the full IceCube Monte Carlo.

Neutrino–nucleon cross-sections

Overall, in the energy range relevant for this analysis (that is, neutrino energies between a few hundred GeV and a few tens of TeV), the uncertainty on the neutrino–nucleon and antineutrino–nucleon cross-sections are about 2–3% and 4–10%, respectively33. This can be parametrized by using different parton distribution functions (PDFs). In this work, our default PDFs are the ones from the HERAPDF set34. Given the degeneracy between the interaction cross-section and the Earth’s density, we expect the propagated uncertainties in the derived quantities to be of the same order. We have checked that this is so and given the much larger statistical uncertainties in the data, the errors in the cross-sections have no significant effect on our results.

Note that taking a complementary approach to the one presented in this Letter, one could try to confirm the value of the neutrino–nucleon cross-section at these energies, assuming the Earth’s density profile to follow the PREM32,35.

Propagation of neutrinos through the Earth

The transport equations for neutrinos traversing the Earth, which we solve using the ν-SQuIDs package36, consist of four main ingredients (see, for example37): the standard evolution Hamiltonian in matter, which includes the vacuum mass-mixing terms and the effect of coherent forward scattering off electrons of the medium, given by the matter interaction potential; the attenuation effect caused by neutrino inelastic interactions with matter, either via charged-current or neutral-current processes; the redistribution of neutrinos from higher to lower energies after neutral-current interactions; and, finally, the neutrino regeneration term from tau lepton decays. Neutrino flavour oscillations in matter, given by the first term, represent the dominant effect for neutrino energies below a few hundred GeV. On the other hand, the other terms become dominant for neutrinos with higher energies. Since the neutrino–nucleon cross-section increases with energy, at these energies the neutrino flux gets attenuated2. In the case of neutral-current interactions neutrinos are degraded in energy38. In the case of charged-current interactions, neutrinos are absorbed and a lepton of the same flavour is produced. Whereas in the case of electron and muon neutrinos (and antineutrinos), the associated lepton (electron or muon) is rapidly absorbed in the Earth and does not contribute to the high-energy neutrino flux, the tau leptons produced after tau neutrino charged-current interactions decay before losing too much energy. In these decays, a new tau neutrino (or antineutrino) with lower energy is produced and thus, they get regenerated39. Moreover, secondary electron and muon neutrinos (and antineutrinos) are also produced after tau lepton decays40. For the energies we consider here, neutrino oscillations are suppressed and the effects of tau neutrino regeneration and secondary production of electron and muon neutrinos are negligible. On one hand, tau neutrinos are rarely produced in the atmosphere, and on the other hand, this effect is important only for spectra much harder than the atmospheric neutrino one. Therefore, for the sake of saving computational time, we have not included the regeneration or secondary production terms in this work. We stress the corrections are much smaller than the precision on the determination of the Earth’s profile achieved with current data.

Earth modelling

We have considered two models for the Earth, fixing the position of the core/mantle boundary and the transition from the inner to the outer core. On one hand, we have parametrized the Earth’s density with a one-dimensional five-layer model, R1,…,R5, with constant density in each of the layers (Fig. 1a). The layers are defined as follows: R1[0,0.195]R, R2[0.195,0.3725]R, R3[0.3725,055]R, R4[0.55,0775]R, R5[0.775,1]R, where R is Earth’s mean radius, R = 6,371 km. The first layer corresponds to the inner core, the second and third layers (of equal thickness) to the outer core, and the final two layers (of equal thickness) to the mantle. The density of each of these layers is allowed to float freely and independently. On the other hand, we have also considered a model with five layers, again, but with a density profile within each layer that follows that of the PREM. The density in each of these layers is multiplied by a factor which is also allowed to vary freely and independently of the others. For the current data set, we do not expect that a larger number of layers would change the results presented here. Indeed, this is partly explained by the similarity of the results of the flat-layer model and the PREM-based model with five layers. Since our aim in this work is to evaluate the sensitivity of the neutrino attenuation effect on the Earth’s density profile, throughout this work, we have not imposed any external constraint on the Earth’s mass or moment of inertia, which are (gravitationally) known much more precisely than what can currently be achieved with neutrinos22,23,24. The results for the different cases are shown in Supplementary Figs. 3,4 and Supplementary Table 1. In order to understand the importance of these external constraints, we have also performed an analysis of the present statistical sample. including the total mass of the Earth and its moment of inertia as external priors. Since approximately 70% of the Earth’s mass is located in the mantle, this procedure, in practice, corresponds to fixing the mantle density to the PREM value within small errors. The 68% credible interval for the upper mantle density changes from ρ5[1.22,4.78] g cm3 for the unconstrained fit to ρ5[4.43,4.79] g cm3 for the constrained one. On the other hand, these priors have a rather small impact on the inner and outer core densities. The 68% credible interval for the average core density from an unconstrained fit is \(\bar \rho _{{\mathrm{core}}}^\nu \in [10.2,20.8]\) g cm3, whereas for the constrained fit we get \(\bar \rho _{{\mathrm{core}}}^\nu \in [9.65,18.6]\) g cm3.

Parameter estimation

To quantitatively assess the power of the one-year up-going muon IC86 sample to determine the Earth’s density profile, we performed a likelihood analysis using all the events in the data sample and characterizing each event by its reconstructed muon energy and zenith angle. The full likelihood is defined as the bin product of the Poisson probability of measuring \(N_i^{{\mathrm{data}}}\) for the expected value \(N_i^{{\mathrm{th}}}\) times the product of Gaussian probabilities for the pulls of the nuisance parameters. The log-likelihood (up to a constant) is given by

$${{\mathrm{ln}}{\cal L}({\mathbf{\rho }};{\mathbf{\eta }})} = {\mathop {\sum}\limits_{i\, \in \,{\mathrm{bins}}} \left( {N_i^{{\mathrm{data}}}{\mathrm{ln}}\,N_i^{{\mathrm{th}}}({\mathbf{\rho }};{\mathbf{\eta }}) - N_i^{{\mathrm{th}}}({\mathbf{\rho }};{\mathbf{\eta }})} \right)} - \mathop {\sum}\limits_j \frac{{\left( {\eta _j - \eta _j^0} \right)^2}}{{2\sigma _j^2}}$$

where the subindex i refers to a bin in reconstructed muon energy \(\left( {E_\mu ^{{\mathrm{rec}}}} \right)\) and cosine of the reconstructed zenith angle \(\left( {{\mathrm{cos}}\,\theta _{\rm z}^{{\mathrm{rec}}}} \right)\); \(N_i^{{\mathrm{th}}}({\mathbf{\rho }};{\mathbf{\eta }})\) is the expected number of events for a given value of the densities in each layer (parameterized by ρ) and the nuisance \(\left( {{\mathbf{\eta }} \equiv \{ N,\pi /K,{\mathrm{\Delta }}\gamma ,{\mathrm{DOM}}_{{\mathrm{eff}}}\} } \right)\) parameters in the ith bin; and \(N_i^{{\mathrm{data}}}\) is the number of data events in the same ith bin. The index j corresponds to the nuisance parameters with Gaussian prior (π/K and Δγ) and σj is the Gaussian error. To compute the likelihood for a given value of the parameters, we first propagate the neutrino fluxes from the atmosphere to the detector for both neutrinos and antineutrinos, then we weigh the events from the IceCube Monte Carlo with the propagated flux, which is a function of the true neutrino energy Eν and the zenith angle θz, and we construct two-dimensional histograms as a function of the reconstructed variables: \(E_\mu ^{{\mathrm{rec}}}\) and \(\theta _{\rm z}^{{\mathrm{rec}}}\) (using 10 bins in muon energy and 60 angular bins).

All the credible intervals we indicate in the main text and in the Supplementary Information correspond to the highest posterior density interval (that is, the shortest interval on a posterior probability density for a given confidence level) for one-dimensional marginalized distributions. As a consequence, these intervals always include the point with the highest posterior density, which we also indicate, as a reference, for each quantity. Unless otherwise stated, the credible intervals are all provided for an integrated 68% posterior probability.

Forecast for ten years of data

It is interesting to get an idea of how the measurements shown in this paper may improve as soon as more data become available. For this reason, we also compare our results with the outcome of an analysis performed assuming ten years of data. For this forecast analysis, we consider the combination of the Honda–Gaisser primary cosmic-ray spectrum with the Gaisser–Hillas H3a correction (HG-GH-H3a) and the QGSJET-II-04 hadronic-interaction model for the atmospheric neutrino flux, and we simulate the future data assuming the PREM density profile. Simulated data are subsequently fitted (using the same atmospheric neutrino flux) with a five-layer model, as in the current analysis using the IC86 sample, albeit with a density following the PREM profile within each layer. This approach is used to avoid coarse binning with higher statistics, in the case of a piecewise flat profile. Although with current data, considering five layers with constant density is equally as good as assuming a more realistic profile such as the PREM model (see Supplementary Figs. 3,4 and Supplementary Table 1), with more data, a finer modelling of the Earth with more than five layers or a more accurate profile within layers would be certainly needed. The results of the comparison are shown in Supplementary Figs. 5,6.

In our default forecast analysis, as described above, we assume that future data will come along with a better determination of the atmospheric neutrino flux model31,41,42,43,44,45,46,47 and that, therefore, a fit of the forecast data can be performed using only the flux model used to generate the data themselves. We recall that uncertainties in the flux models will also be reduced by other complementary future measurements (such as, for example, the measurement of the atmospheric muon flux, improved cosmic-ray measurements, better understanding of hadronic interactions, measurements of atmospheric neutrino fluxes at lower energies, and even at similar energies for down-going neutrinos of different flavours). For the forecast to properly take into account potential improvements on the ingredients of the analysis, other different type of data would certainly have to be included, going beyond the scope of this Letter. Nevertheless, we have studied the impact of the discrete choices of primary cosmic-ray spectrum and hadronic-interaction model on the forecast. We have used different combinations to generate and to fit mock data. We have found that, for some fluxes, the results of the fit show some disagreement with the gravitational measurement of the Earth’s mass and of the Earth’s moment of inertia (whereas the Earth’s core mass and the core/mantle density jump are little affected by the choice of flux model), which translates into systematic uncertainties. In particular, independently from the combination of primary cosmic-ray spectrum and hadronic-interaction model used to generate the data, we have found that knowledge of the latter would reduce the error on the Earth’s mass and moment of inertia by 30–40%. On the other hand, for a given hadronic-interaction model, additional knowledge about the primary cosmic-ray spectrum may reduce uncertainties by only, at most, a few percent. Therefore, we can conclude that a good understanding of the hadronic-interaction model is important to reduce errors significantly, whereas uncertainties due to our current knowledge of the cosmic-ray primary spectrum have a rather small impact.

We have also performed more detailed ten-year forecast analyses, considering different density profiles within each layer (either flat or following the PREM) and several configurations of layers. From these analyses we have verified that: the statistical error in the outer mantle layer could go down to a few percent; the statistical error in the inner mantle layer will get reduced down to around 10%; a finer description (more layers) of the one-dimensional Earth’s density profile than the one used in the present work will be needed. It is not yet clear if with ten years of data it will be possible to determine the location of the core/mantle boundary just by looking at high-energy neutrino data, but what is clear from the forecast analysis is that a simple five-layer Earth’s model would not be the optimal one to analyse the data, and more layers would represent a better description of the density profile. For example, we have checked that the results of a five-layer fit would be affected by the choice of the profile within layers, such as, for instance, flat layers (constant density within each layer) versus layers with a density profile that follows the PREM.

We stress that the study in this Letter is based on an event sample with maximum muon energies of ~20 TeV. For future updates, adding events with slightly higher energies would certainly be useful: since for energies of ~40 TeV the neutrino mean-free path corresponds to the Earth’s diameter, extending the sample up to ~100 TeV could help improving the results. Actually, the IC86 sample already includes some events for which substantial Earth absorption takes place, as the reconstructed muon energy is smaller than the true neutrino energy and, therefore, most high-energy neutrino events end up into lower muon energy bins. However, due to the steeply falling atmospheric neutrino spectrum, going to much higher energies would significantly reduce the already low neutrino flux at the detector. Moreover, the higher the energies, the more likely the observed events are dominated by the astrophysical neutrino flux, which has, in principle, a (yet unknown) zenith distribution different from the atmospheric flux, and this introduces an additional uncertainty. Incidentally, events with higher energies than the ones used in this Letter have already been considered to determine the neutrino–nucleon cross-section32,36.

Code availability

The codes used in this work are publicly available. For the propagation of neutrino fluxes through the Earth we use the ν-SQuIDs package36, which can be downloaded from Once the neutrino fluxes are converted into reconstructed energy and zenith angle distributions using the Monte Carlo detector simulation provided along with the data, the statistical analyses are performed by means of the MultiNest nested sampling algorithm48,49,50, which can be freely downloaded from The Multinest library implements a multi-modal nested sampling algorithm which is intended to compute the Bayesian evidence for a given likelihood function and priors. The algorithm produces a weighted Markov chain that can be used to reconstruct the marginalized posterior probability of the parameters. For that reconstruction we use the GetDist Python library (, which is one of the most efficient and robust tools for Bayesian analyses for problems with a large number of parameters, as in our case.

The IceCube data we consider in this paper are the same sample used by the collaboration to search for resonant matter effects induced by light sterile neutrinos8. The Monte Carlo results used to simulate the detector characteristics and all data are publicly available and can be downloaded from

A.D. thanks G. Cultrera, D. Errandonea, A. Kavner, C. Piromallo and G. Soldati for useful discussions. A.D. and J.S. were supported by the Generalitat Valenciana under grant PROMETEO II/2014/050 and by the Spanish MINECO grants FPA2014-57816-P and FPA2017-85985-P. S.P.-R. is supported by the Generalitat Valenciana under grant PROMETEOII/2014/049, by the Spanish MINECO grants FPA2014-54459-P and FPA2017-84543-P, by a Ramón y Cajal contract, and also partially by the Portuguese FCT through the CFTP-FCT Unit 777 (PEst-OE/FIS/UI0777/2013). The authors also acknowledge support by the Spanish MINECO under grant SEV-2014-0398. J.S. is also supported by the Spanish MINECO grant FPA2016-76005-C2-1-P, María de Maetzu program grant MDM-2014-0367 of ICCUB and research grant 2017-SGR-929. All authors are supported by the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreements No. 690575 and 674896.