Atmospheric Neutrino Oscillations for Earth Tomography

Modern proposed atmospheric neutrino oscillation experiments, such as PINGU in the Antarctic ice or ORCA in Mediterranean sea water, aim for precision measurements of the oscillation parameters including the ordering of the neutrino masses. They can, however, go far beyond that: Since neutrino oscillations are affected by the coherent forward scattering with matter, neutrinos can provide a new view on the interior of the earth. We show that the proposed atmospheric oscillation experiments can measure the lower mantle density of the earth with a precision at the level of a few percent, including the uncertainties of the oscillation parameters and correlations among different density layers. While the earth's core is, in principle, accessible by the angular resolution, new technology would be required to extract degeneracy-free information.


Introduction
Using neutrinos for Earth tomography is a dream much older than modern oscillation physics, see Ref. [1] for a review: Early proposals exploit the increase of the neutrino cross sections with energy, leading to significant neutrino absorption over the earth's diameter for energies larger parts atmospheric neutrino oscillations are most sensitive within this scenario. We will use proposed experiments such as PINGU ("Precision IceCube Next Generation Upgrade") [28] in the Antarctic ice or or ORCA ("Oscillation Research with Cosmics in the Abyss") [29] in Mediterranean sea water, which are modern megaton-sized neutrino oscillation experiments designed for neutrino oscillation precision measurements with leading sensitivity to the neutrino mass ordering -and thus the Earth matter effect; see the Appendix for the simulation techniques.
Earlier discussions in that direction include the matter effect sensitivity [30] and the sensitivity to the core composition [31].

Model and Methods
We propose a whole-Earth model with seven different density layers adopted from the Preliminary Reference Earth Model (PREM) profile [26], which is shown in Fig. 1, to identify the regions with highest sensitivity. We split the PREM profile into seven layers at depths d, where the characteristic density jumps occur (cf., solid curves in Outer core (6), 2860 km d 5151 km, Inner core (7), 5151 km d R E = 6371 km (R E : Earth radius). Note that compared to seismic waves, which tend to be reflected or refracted at density jumps, neutrino oscillations are not very sensitive to structures or even strong gradients shorter than the oscillation length [23], and therefore cannot resolve the density jumps precisely.
Therefore it is reasonable to adopt this knowledge from geophysics.
Each baseline (see rays in Fig. 1) is separated into sections going through the density layers. Within each density layer, we follow the PREM profile [26], where the matter profile is discretized into a sufficient number of steps with constant density. The oscillation probabilities are then evaluated with the evolution operator method (see e.g. Ref. [32]): the initial state |ν α is propagated through all matter density slices with thicknesses x j and constant densities ρ j through all crossed layers by as the Hamiltonian within each layer H is not explicitly time-dependent. The transition probability then reads P αβ = | ν β |V(x n , ρ n ) . . . V(x 1 , ρ 1 )|ν α | 2 .
Note that in general which means that the different operators do not commute and the probability will depend on the order the layers are traversed. This is an important difference to X-ray or absorption tomography, which is only sensitive to the path-integrated attenuation.
Our measured quantity in each density layer is actually a factor linearly re-scaling the density profile in this layer, as the actual density profiles for different baselines are slightly different even if they cross the same layer. Although this model is an approximation, it yields similar results for the relative matter precision compared to alternatives (such as choosing the density within each layer to be constant), but maintains accuracy of the oscillation probabilities and the oscillation measurements for the more realistic PREM profile. For convenience, we call the measured scaling factor for layer i "ρ i /ρ i " , and depict it as error on the average matter density.
Note that since neutrino oscillations are not sensitive to structures or changes shorter than the oscillation length [23], additional parameters, such as multiple layers or gradients in the layers, cannot be resolved anymore beyond that level. It is clear that similar arguments apply to individual geophysical techniques, such as using the earth's free oscillation modes, see Refs. [33,34]. As a consequence, "structural" information from neutrino oscillation tomography has to rely on strong density jumps (leading to interference in the probabilities) or different baselines, and "average" information has to rely on some knowledge from geophysics on scales shorter than the oscillation length (smoothing the density profile). Since new ways to combine neutrino oscillations with -or compare them to -geophysical results require further research, and atmospheric oscillation tomography is limited by the complexity from the number of parameters (oscillation parameters, systematics, and geophysics parameters), we choose the approach introduced above.
Furthermore, note that we do not include constraints on the total mass and rotational inertia of the earth, which means that (technically speaking) some of our variations would violate these important constraints. However, in order to include these, one needs to define a correction scheme, i.e., which layers are corrected for density variations to maintain these constraints. One possibility has been discussed in Ref. [20]: Since changes of the innermost densities of the earth (e.g., inner core) influence mass and rotational inertia less that the outermost parts (where the volume is much larger), one can use small adjustments of the outer densities to compensate for large density changes in the innermost earth in spite of the higher densities there. Since it is clear that the final result would depend on that correction scheme, and additional constraints would rather improve our result than deteriorate it (in a similar way as the free oscillation result [33]), we do not consider the total mass and rotational inertia constraints in this work. An alternative (but computationally more expensive) approach would be to generate very different fit density profiles from the very beginning, and define a measure how well they fit neutrino oscillations and other potential constraints [33,23].
The precision on ρ i /ρ i is obtained by minimizing the ∆χ 2 over all oscillation parameters, auxiliary systematics parameters, and the other ρ j /ρ j (j = i) simultaneously. We also impose a 30% external constraint on ρ j /ρ j for j = i, i.e., we assume that there is some crude knowledge on the other layer densities from geophysics and whole-Earth constraints. From the geophysics perspective, this is a very coarse constraint. From the particle physics perspective, it has the advantage that it prevents the n-dimensional line minimization techniques used for the analysis from falling into unphysical solutions, such as negative densities (here the penalty χ 2 would exceed nine). It does not have any significant consequences for the result, except from the outer core density measurement which suffers from correlations with the inner core density.
In order to illustrate the underlying physics, consider a simple example using neutrino oscillations in constant matter density. The oscillation probability P µe = P (ν µ → ν e ) can (neglecting contributions from solar terms) be approximated as P µe sin 2 θ 23 sin 2 (2θ 13 ) sin 2 ∆m 2 31 L 4E .
This probability is (apart from the factor sin 2 θ 23 ) just a two-flavor oscillation probability, where the fundamental parameters ∆m 2 31 and θ 13 are replaced by effective parameters in matter ∆m 2 31 = ξ · ∆m 2 31 and sin(2θ 13 ) = sin(2θ 13 )/ξ with the mapping parameter ξ ≡ sin 2 (2θ 13 ) + (cos(2θ 13 ) −Â) 2 and the matter potentialÂ ≡ ±2 √ 2G F n e E/∆m 2 31 ; the different signs refer to neutrinos (plus) and antineutrinos (minus). Here the quantity of interest is the electron density in Earth matter n e , which can be converted into the matter density by n e = Y e ρ/m N using the electron fraction Y e (number of electrons per nucleon) and the nucleon mass m N . While one has for hydrogen Y e = 1, heavier materials prefer Y e 0.5 because of approximately equal numbers of protons and neutrons. We fix Y e = 0.5 in this study, but one should keep in mind that one actually measures the product of Y e × ρ. 1 It is easy to see that the conditionÂ → cos(2θ 13 ) minimizes ξ, leading to effective maximal mixing. This case is often referred to as "matter resonance", and can be re- From this figure, one can immediately see that excellent sensitivity is expected to the Lower Mesosphere (layer 5). Although inner core and outer core can be, in principle, resolved, the corresponding data will be smeared over direction, the covered solid angle (the event rate is proportional to) is smaller, and the relevant energies are close to the experiment threshold. While 1 The allowed range for Y e is actually small for typically used geophysical composition models -which implies that the composition is much harder to measure than the matter density. The reason is that heavier stable nuclei typically contain similar numbers of protons and neutrons -as long as there is no significant hydrogen content. A well studied example in that context is the outer core, see Ref. [31], Table 1: The values of Y e vary at the level of one percent -which is beyond the relative precision we find in this study.  these points can be illustrated with the simple constant matter approach, the realistic matter profile of the earth leads to interesting interference effects and a parametric enhancement coming from the oscillation length matching the mantle-core-mantle structure of the earth [35,36], see also Ref. [37], which are treated numerically. Additional complications are the composition of the atmospheric neutrino flux, containing both electron and muon flavors, and the inability of the detectors to discriminate neutrinos from antineutrinos; see e.g. Refs. [38,39] for details.
We use two event samples (muon track-and cascade-like) for the analysis, including all these effects; for analysis details, see the Appendix.
We point out that a "proof of principle" for the independent extraction of the layer densities requires an experiment simulation including systematics, correlations with the oscillation parameters, and correlations among the layer densities in a self-consistent framework, see the Appendix, which is novel in this work. The simulation techniques are based on Ref. [40] using an extended version of the GLoBES ("General Long Baseline Experiment Simulator") software [41,42], which can handle the required level of complexity. Figure 2: Experiment sensitivity to matter density. Projected experiment precision (1σ error bars) for PINGU (left) and ORCA (right) after ten years of data taking for the matter density layers corresponding to Fig. 1. Here the normal mass ordering best-fit values are assumed, and correlations (with systematics, oscillation parameters, and other layer densities) are taken into account. The solid curves correspond to the PREM matter density profile [26].

Results
For the matter density measurement, one can adopt two viewpoints: a) Tomography approach: what parts of the earth are atmospheric neutrino oscillations most sensitive to? b) Precision approach: suppose that better geophysical information exists on some layers, with what precision can a specific density be extracted? To address a), we show in Fig. 2 (see also Tab. 1) the expected precision including systematics and correlations with oscillation parameters and other matter densities. Since the zenith angle resolution (Fig. 1) prohibits a resolution of layers 1, 2, no sensitivity can be obtained, and the sensitivity to layers 3 and 4 is weak.
The best precision is found in the lower mantle with 5% and 4% for PINGU and ORCA, respectively. The corresponding ∆χ 2 is shown in the left panel of Fig. 3 for ORCA: it is wellbehaved Gaussian and correlations with other density layers are not important. This result may, at a first glance, not be too exciting compared to the collective constraints from geophysics in- Figure 3: Impact of parameter degeneracies. Here the log-likelihood ∆χ 2 is shown as a function of the relative error on the average density for three different layers (in different panels) for ORCA. Solid curves include correlations among different matter density layers, with systematics, and oscillation parameters, whereas dashed curves do not include the matter density layer correlations. The horizontal lines correspond to 1σ, 2σ, and 3σ for a Gaussian ∆χ 2 .
cluding free oscillations, total mass, and moment of inertia of the earth, which are believed to constrain the mantle density at the per cent level [33,34] -although the statistical interpretation of these precisions (confidence level of the error) seems less straightforward than in the present case. We have nevertheless demonstrated that neutrino oscillations can contribute at a similar level with an independent technique and different systematics. They may even be competitive to Ref. [33] if the whole-Earth constraints (mass, rotational inertia) are included, and our method does not rely on the assumption of linearized perturbation theory as Ref. [34]. Future tests of neutrino tomography may use similar techniques for better comparisons, which are, however, currently subject to computational constraints. Further applications may include the test of ambiguities and structures, such as the seismic wave-inferred low shear velocity provinces (LLSVPs) or ultra-low velocity zones (ULVZs) in the lower mantle. lation with the inner core density (the event rates mix within the zenith angle resolution). While the 1σ precision (solid curve) roughly corresponds to the one obtained for the core composition estimate in [28], it is clear that the shown degeneracies prohibit a self-consistent extraction of the outer core density up to higher confidence levels. This result applies to the chemical composition measurement as well as, in comparison to Ref. [31], detector setups closer to the experimental proposals are used, and the densities of the other layers are left free.
If, however, viewpoint b) is adopted, the dashed curve will represent the core density measurement, and the impact of correlations is reduced. The intrinsic oscillatory structure in Fig. 3 remains, as illustrated in Fig. 4 for three different values of ρ/ρ for P µe , corresponding to the dot marks in Fig. 3 (middle panel). The oscillation peak at E 6 GeV is almost perfectly rematched by ρ/ρ = 1.55, while it is very different for ρ/ρ = 1.25. The low energy differences are more difficult to resolve due to the smaller effective mass and poorer directional and energy resolutions for lower energies. On the other hand, the lower bound on the outer core density is robust, as lower densities correspond to higher resonance energies where the effective masses of the detectors increase.
We chose the NO best-fit earlier in this study; however, the actual oscillation parameters chosen by Nature may be different. We therefore show the result for ORCA and the IO in Note, however, as both the mass ordering sensitivity and the matter effect sensitivity scale with terms ∝ sin 2 θ 23 in the appearance oscillation channels, the performance for the earth density measurements will scale in a similar way to Fig. 10 with this parameter. This means that the actual result could be much better depending on the oscillation parameter values chosen by Nature. An example is shown in the right panel of Fig. 5 for ORCA for a parameter set within the 3σ currently allowed range. Here a precision of better than 3% is obtained for the lower mantle density. In this case, the outer core density can be actually measured with a precision better than 5%, and the degeneracies can be resolved at almost 2σ.
While the inner core density may be the prime target from the geophysics perspective, as it is the most difficult to access, currently planned instruments do not allow for a high confidence level extraction even if all the other densities were known (see dashed curve in right panel of Fig. 3). This measurement operates close to the detection threshold, where also energy and zenith angle resolutions are weaker, and it suffers from a very small solid angle covered by the inner core. A more densely instrumented detector, such as proposed in [43,44], would have a lower threshold and potentially better low energy directional and angular resolutions helping both the inner and outer core density measurements. Especially in combination with geophysical data on the outer core, an extraction of the inner core density may then become possible.

Summary and Conclusions
We have demonstrated that atmospheric neutrino oscillations measured by planned detectors can provide excellent sensitivities to the lower mantle density and give a robust lower bound on the outer core density. The obtained information is complementary to that of seismic waves, as different quantities (electron density versus seismic wave velocity) and different propagation paths (straight lines versus refracted curves) are tested, and the underlying systematics are very different.
Finally, neutrino oscillation tomography is yet a very young discipline which only has become feasible after the discovery of a non-zero value of θ 13

Simulation Methods and Mass Ordering Sensitivity
The primary physics target for the PINGU [28] and ORCA [29] experiments is the mass ordering determination [38], i.e., the question if the mass eigenstate m 3 is lighter or heavier than m 1 and m 2 . The mass ordering can be measured by matter effects, as the resonance condition A → cos(2θ 13 ) can be only implemented for neutrinos and sgn(∆m 2 31 ) = +1 (normal ordering, NO) or antineutrinos and sgn(∆m 2 31 ) = −1 (inverted ordering, IO); see definition ofÂ in main text. Therefore, the normal mass ordering will lead to an enhancement of the oscillation effect for neutrinos and suppression for antineutrinos, and the inverted mass ordering to an enhancement for antineutrinos and suppression for neutrinos. Note that PINGU and ORCA cannot distinguish between neutrinos and antineutrinos directly, but have to rely on flux and cross section differences; the India-based Neutrino Observatory (INO) [46,47] is a different proposal which can discriminate between neutrinos and antineutrinos by magnetization of an iron detector, but the detector mass is much smaller. Although we propose a spin-off of the main physics target in the main text, we need to establish the mass ordering determination in a self-consistent framework in order to demonstrate that the matter profile sensitivity is guaranteed even without extra equipment.
We therefore show that we can reproduce the mass ordering sensitivities of the experimental collaborations, including the uncertainties of all oscillation parameters, such as δ CP . We treat all 6 oscillation parameters, 7 matter densities, and, using the pull method, 12 auxiliary systematics parameters equally. The precision for one parameter, such as a matter density or an oscillation parameter, can be obtained by projecting the resulting 25-dimensional fit manifold onto a one-dimensional sub-space by minimization of the ∆χ 2 over all parameters not shown.
Most importantly, this framework is fully self-consistent in the sense that any measurement of the matter density is consistent with the measurement of the oscillation parameters, which can

Common Simulation Framework
We simulate PINGU and ORCA for the first time within an identical oscillation framework, the same binnings, the same systematics implementation and parameters, the same definition and computation of the performance indicators as for long-baseline experiments such as LBNF-DUNE [52], and the same Earth density profile, using an extended version of the GLoBES software [41,42]. That GLoBES version allows for user-defined, channel-based systematics treatment across experiment boundaries, which was first applied in Ref. [48]. The zenith angle smearing (redistribution of events) between incident θ z and reconstructed θ z is computed by an oscillation-channel dependent migration matrix R ijk ≡ R(E i , (θ z ) j , θ k z ) (i, j, k: bin indices) as a function of neutrino energy E i . This implies that upgoing events exceeding θ z = 180 • are reconstructed in the corresponding zenith angle bin under a different azimuth. The advantage of this method is that the zenith angle resolution does not become asymmetric at the boundaries, and Gaussian behavior is better reproduced. We interpret the zenith angle resolutions ∆θ z in terms of a normalized Gaussian in order to compute the migration matrix integrated over the θ z range covered by zenith angle bin j. The energy smearing between incident E and reconstructed E energy is described by an energy smearing matrix S ij ≡ S(E i , E j ).
We parameterize this matrix with a Gaussian similar to Eq. (4), such that integrated over the E range of the ith energy bin, unless noted otherwise. Note that we interpret the median zenith angle and energy resolutions given by the experimental collaborations as ∆θ z (E) and ∆ E (E), which reproduces their sensitivities sufficiently well. This interpretation has limitations depending on the structure of the actual event migration matrices, such as non-Gaussianities, which can be only addressed by the experimental collaborations. However, it has the advantages that it can be identically applied to both experiments and allows for comparability and transparency of the assumptions.
For the simulation, we use two event samples (muon tracks and cascades), which are a combination of a number of channels. For instance, one "channel" corresponds to ν e → ν e (CC), including a specific source flux, cross sections, efficiencies, and fiducial volume; altogether there are 32 such oscillation channels (from the two initial flavors ν e , ν µ into the three final flavors ν e , ν µ , ν τ makes six, times two for neutrinos and antineutrinos plus four non-oscillating neutral current channels, makes 16, times two because separate channels with separate efficiencies for the two event samples are needed, makes 32). Each neutrino event is counted either in the "right" channel as signal, or in the "wrong" channel as background, depending on the (mis-)identification probabilities.
Systematics are treated exactly in the same way as in long-baseline simulations. It is important to note that not only the values for the systematical errors are important, but also how systematics are correlated among different channels and bins. For example, one may not know a certain cross section, but one does know that the same value has to be used everywhere the same interaction is measured. These kind of correlations are self-consistently implemented, see Ref. [48] for details. We list the considered systematics in Tab. 2, together with the assumed values; each of these corresponds to one auxiliary parameter. For example, the cross section × fiducial mass errors for different event types are assumed to be known (externally measured) to about 5% in the considered energy range, see e.g. discussion in Ref. [48]. Note that at this point we adopt the identical systematics implementation and errors for both experiments, following physical arguments (e.g., cross section uncertainties); however, one may think about alternative, more inclusive concepts, see e.g. Ref. [53].
For the atmospheric neutrino flux, we use updated versions [54], where we use the azimuth- δ CP = 254 • . Note that these solutions lie in different octants, and that δ CP is slightly different.
We compute the ∆χ 2 for a certain true mass ordering as the minimal ∆χ 2 over all oscillation parameters with the other mass ordering. This definition is exactly the same as the one used for long-baseline experiments, such as LBNF-DUNE. It has two implications: First of all, the symmetrical (∆m 2 31 ) eff for the ν µ disappearance channels [56,57,58,59,30,60] (∆m 2 31 ) eff = ∆m 2 31 − ∆m 2 21 (cos 2 θ 12 − cos δ CP sin θ 13 sin 2θ 12 tan θ 23 ) , clearly indicates that the mass ordering sensitivity must, in general, depend on the fit value of δ CP . This has been demonstrated for PINGU, see Fig. 3 in Ref. [40]. Second, it is well known that for non-maximal atmospheric mixing, the minimal χ 2 may be either found in the wrong ordering-right octant (fit θ 23 similar to true θ 23 ) or in the wrong ordering-wrong octant (fit θ 23 π/2− true θ 23 ) region. Note that some of the published documents of the experimental collaborations do not yet include all of these effects, which are however essential for a fair assessment of the matter profile sensitivity.

Experiment-Dependent Specifications
Experiment-specific assumptions for this study include the effective mass (for different event types), the mis-identification probabilities between muon tracks and cascades, and the zenith angle and energy resolutions.  in consistency with Ref. [27]. The corresponding quantities are shown in Fig. 7. While the effective mass (upper left panel) for all cascade channels is directly obtained from Ref. [63], the effective mass for muon tracks for this configuration has been adopted from Refs. [29,61]. The flavor identification information (upper right panel), comes from Ref. [29]. The zenith angle resolution in the lower left panel can for muon tracks be directly obtained from Refs. [29,61], whereas the resolutions for the electromagnetic cascade channels have been adopted from Ref. [63]; these resolutions are also shown in Ref. [27]. The assumed resolutions for neutral currents and hadronic cascades are assumed to be similar to muon tracks for low energies, and about a factor of two weaker for high energies. Note that the zenith angle resolutions are not parameterized, but instead interpolating functions are used to pre-compute the migration matrices.
The energy resolutions have been obtained in a similar way from Refs. [29,61] for muon tracks, and have been adopted from Ref. [63] for cascades, see lower right panel. We assume that they Comparing the lower rows between Fig. 6 and Fig. 7, our assumptions imply slightly better zenith angle and energy resolutions for ORCA than PINGU at least for ν (ν e ) cascades and lower energies, which means that one may expect slightly better sensitivities for the matter density measurement -which we find in the main text.

Sensitivity to Mass Ordering
We show the sensitivity for the NO and IO best-fits in Fig. 8. In order to compare to the existing literature, we have fixed δ CP to the respective best-fit value for the dashed curves; these reproduce the official versions from the experimental collaborations [29,64] very well. For the solid curves, we fully take into account the minimization over δ CP , which can somewhat affect Figure 9: Mass ordering sensitivity as a function of θ 23 . Number of σ (assumed to be χ 2 ) for the PINGU (left) and ORCA (right) as a function of the true sin 2 θ 23 . Solid curves correspond to the normal ordering, dashed curves to the inverted ordering. The best-fit (BF) values are marked as well. Three years exposure assumed, matter densities are assumed to be known.
the sensitivities depending on the parameter point [40]. The final sensitivities for the NO are in fact surprisingly similar for PINGU and ORCA, leading to a 3σ discovery after about three to four years of operation within the identical systematics and oscillation framework; for the IO the required times are somewhat longer.
While the solid and dashed curves in Fig. 8 assume the matter profile of the earth to be precisely known, we also perform a self-consistent simulation with the neutrino Earth model described above (dotted curves). Interestingly, we find almost identical performances for PINGU and ORCA for the normal mass ordering. The PINGU sensitivity for the inverted ordering is most affected by the unknown matter density. An inspection of the systematics pulls reveals an 18% increase of the lower mantle density and a 6% decrease of the outer core density (compared to the input densities) for the three year sensitivity. Here the external knowledge from geophysics may in fact be essential to improve the mass ordering sensitivity.
An important cross-check is the dependence of the mass ordering sensitivities on the true θ 23 , see Fig. 9. These figures reproduce the qualitative behavior of the experimental collabora- Figure 10: Required running time to establish the mass ordering at 3σ. The contours show the approximate required running time [yr] to reach 3σ as a function of true δ CP and true θ 23 within the currently allowed parameter space at 3σ [55]. Matter densities are assumed to be fixed.
tions [27] but are somewhat more conservative for ORCA because the correlation with δ CP is fully included. In some cases jumps between the octants where the minimal χ 2 found, see e.g. right panel around sin 2 θ 23 = 0.43, which are somewhat sensitive to the experiment implementation for small θ 23 .
Finally we show in Fig. 10 the required running time (contours, in years) to establish the mass ordering at 3σ as a function of true δ CP and true θ 23 . For this figure, the performance has been computed as a function of δ CP and θ 23 simultaneously for three years of operation. Then we have linearly scaled the χ 2 to obtain an estimate for how long it takes to reach 3σ. As a result, a 3σ discovery is guaranteed for PINGU and ORCA within the anticipated timescale of the matter density measurement even in remote (allowed) regions of the parameter space, whereas in the most optimistic case, ORCA can find the NO already after one year. As a consequence, the mass ordering will be resolved in either case at the timescale of the matter density measurement discussed in the main text, and the corresponding degeneracy cannot affect these results.