# Unveiling the outer core composition with neutrino oscillation tomography

Lukas Maderer  Université Paris Cité, CNRS, Astroparticule et Cosmologie,
F-75013 Paris, France Edouard Kaminski Contact<EMAIL_ADDRESS>Université
Paris Cité, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris,
France  João A. B. Coelho  Université Paris Cité, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France Université Paris-Saclay, CNRS/IN2P3, IJCLab,
F-91405 Orsay, France  Simon Bourret  Université Paris Cité, CNRS,
Astroparticule et Cosmologie, F-75013 Paris, France Véronique Van Elewyck
Contact<EMAIL_ADDRESS>Université Paris Cité, CNRS, Astroparticule et
Cosmologie, F-75013 Paris, France Institut Universitaire de France, F-75005
Paris, France

###### Abstract

In the last 70 years, geophysics has established that the Earth’s outer core
is an FeNi alloy containing a few percent of light elements, whose nature and
amount remain controversial today. Besides the classical combinations of
silicon and oxygen, hydrogen has been advocated as the only light element that
could account alone for both the density and velocity profiles of the outer
core. Here we show how this question can be addressed from an independent
viewpoint, by exploiting the tomographic information provided by atmospheric
neutrinos, weakly-interacting particles produced in the atmosphere and
constantly traversing the Earth. We evaluate the potential of the upcoming
generation of atmospheric neutrino detectors for such a measurement, showing
that they could efficiently detect the presence of 1 wt% hydrogen in an FeNi
core in 50 years of concomitant data taking. We then identify the main
requirements for a next-generation detector to perform this measurement in a
few years timescale, with the further capability to efficiently discriminate
between FeNiH and FeNiSixOy models in less than 15 years.

##  1 Introduction

The determination of the composition of the Earth core is a long-standing
debate in Earth sciences (e.g. [1, 2, 3, 4]). Seismology combined with
experimental petrology shows that the presence of light elements is required
to account for the density and seismic velocity profiles in the core, e.g.,
PREM [5]. But the precise nature of the light elements remains elusive, and
various combinations of, e.g., Si and O [6, 7], can equally fit PREM.
Hydrogen, the most abundant element in the proto-solar nebula, is the only
candidate that could account both for the density and velocity profiles in the
core without the need for any additional light element [8, 9, 10]. An FeNiH
core, rather than a more classical FeNiSixOy composition, would bear profound
consequences for Earth formation scenarios and modeling of convection and
generation of the magnetic field in the outer core. The aim of the present
study is to show how the development of high-performance detectors of
atmospheric neutrinos opens a new path to independently constrain the
composition of the core and test the FeNiH hypothesis in particular.

Neutrinos are neutral elementary particles that exist in three flavors:
electron ($e$), muon ($\mu$) and tau ($\tau$). Because they feel only the weak
interaction, they can traverse and emerge from very dense media, including the
Earth itself. The interactions of cosmic rays with the atmosphere generate an
abundant and ubiquitous flux of energetic neutrinos, which undergo flavor
oscillations while crossing the Earth [11]. The probability of their flavor
transition depends on the neutrino energy and path length, but also on the
density of electrons $n_{e}$ in the traversed medium [12, 13]. Because of its
dependence in $n_{e}$, this effect is key in accessing chemical properties of
the traversed medium through neutrino oscillation physics. Assuming an Earth’s
radial structure and composition, $n_{e}$ as a function of $r$, the radial
distance from the center of the Earth, is given by

$n_{e}(r)=\frac{Z}{A}(r)\times\rho(r)\times\mathcal{N}_{A}$ (1)

where $\mathcal{N}_{A}$ is the Avogadro number, $\rho$ is the mass density,
and $Z/A$ is the proton-to-nucleon ratio defined as:

$\frac{Z}{A}=\sum_{i}w_{i}\frac{Z_{i}}{A_{i}}$ (2)

where $Z_{i}$, $A_{i}$ and $w_{i}$ are respectively the local atomic number,
standard atomic weight and weight fraction of element $i$ of the material.
Chemical elements with different Z/A will therefore generate distinct signals
in the neutrino oscillation pattern for a given mass density. The signature of
these oscillations in the atmospheric neutrino signal observed at a detector
can thus reveal the nature of the matter they interacted with.

While this method of Earth tomography has been conceptually explored in
different contexts since the 1980s (see e.g. [14, 15, 16, 17, 18, 19, 20] and
references therein), only now does the upcoming generation of atmospheric
neutrino detectors provide a concrete opportunity to evaluate its capability
of probing the density and/or composition of the deep Earth. Oscillation
tomography with atmospheric neutrinos has been recently discussed by several
authors [21, 22, 23, 24, 25, 26, 27], but its relevance and potential for the
characterization of the Earth’s core composition has not been systematically
explored and discussed yet from the Earth sciences point of view. We propose
to do so in this study, with a specific focus on the capacity of this method
to tackle the question of Hydrogen in the core, thanks to the very large Z/A
of H compared to Fe and other light elements.

Section 2 presents our main results on neutrino tomography of the Earth’s
outer core. In section 2.1 we start by quantifying the expected theoretical
reach of the method under the hypothesis of perfect detector capabilities,
showing how different assumptions for the outer core composition will indeed
modify the number of neutrinos of different flavours interacting at a detector
site. In section 2.2 we move to investigating the realistic performances of
the two main families of atmospheric neutrino detectors presently under
construction, and discuss their ability to detect the presence of 1wt%
Hydrogen in the outer core. Finally, we illustrate in section 2.3 the
evolution required for the next generation of detectors to be able to
discriminate FeNiH vs. FeNiSixOy models of the Earth’s core. The methods and
inputs used for the computations that support the present study are described
in detail in section 3. The context and implications of our findings are
further discussed in section 4.

##  2 Neutrino tomography of the Earth’s outer core

###  2.1 Theoretical approach

Our first step is to evaluate the intrinsic sensitivity of the method by
conducting detailed computations of the propagation of neutrinos through the
Earth and comparing the expected rate of events in a perfect detector under
different assumptions for the outer core composition. In this section we limit
ourselves to the key theoretical concepts of the method, required to
understand its application to the tomography of the core; a more detailed
presentation of neutrinos oscillations and propagation in matter is provided
in section 3.

When they cross the Earth, atmospheric neutrinos experience resonant
oscillations at a depth that depends on their energy and on the electron
density $n_{e}(r)$ as defined in equation 1. For core-crossing neutrinos, this
resonance happens at energies around 3 GeV, as further discussed in section
3.1. For a given flavor $\alpha\,(=e,\mu,\tau)$ observed at the detector, the
information on $n_{e}(r)$ is encoded in the rate of interacting events
$R^{int}_{\alpha}(E,\theta_{z})$:

$R^{\text{int}}_{\alpha}(E,\theta_{z})\equiv\frac{{\mathrm{d}}^{4}N^{\text{int}}_{\alpha}(E,\theta_{z})}{{\mathrm{d}}E\,{\mathrm{d}}\theta_{z}\,{\mathrm{d}}t\,{\mathrm{d}}M}=\left(\sum_{\beta=e,\mu}\frac{{\mathrm{d}}^{2}\Phi_{\nu_{\beta}}}{{\mathrm{d}}E\,{\mathrm{d}}\theta_{z}}(E,\theta_{z})\cdot
P_{\nu_{\beta}\to\nu_{\alpha}}(E,\theta_{z})\right)\cdot\frac{\sigma^{\text{int}}_{\nu_{\alpha}}(E)}{m_{N}}$
(3)

where we have made explicit the dependence on the neutrino energy $E$ and on
incoming zenith angle $\theta_{z}$, as defined in figure 1 (left), under the
assumption of a spherically symmetric mass density profile of the Earth. The
angle $\theta_{z}$ is directly related to the path length $L$ through the
Earth: $L\approx-2R_{\oplus}\cos{\theta_{z}}$ with $R_{\oplus}$ the radius of
the Earth. Extracting tomographic information on $n_{e}(r)$ in a given Earth
layer therefore requires a detailed knowledge of the flavor, path length and
energy distribution of neutrinos having crossed that layer.

In this equation, $R^{\text{int}}_{\alpha}(E,\theta_{z})$ represents the
differential rate of $\nu_{\alpha}$ interactions (with target nucleons of mass
$m_{N}$) at the detector location, as a function the energy and zenith angle
and per unit exposure (defined as the product of running time and target mass
of the detector). It is obtained as a product of the incident (differential)
fluxes of atmospheric neutrinos $\Phi_{\beta}$ (with $\beta=e,\mu$ as the
atmospheric $\nu_{\tau}$ component is negligible), the flavor oscillation
probability along each neutrino path $P_{\nu_{\beta}\to\nu_{\alpha}}$, and the
neutrino-nucleon cross-section $\sigma^{int}_{\nu_{\alpha}}$, that quantifies
the probability for a neutrino to interact (hence to generate a potentially
detectable signal). A quantitative example is shown in the right panel of
figure 1, illustrating both the oscillatory effects and the attenuation of the
flux at high energies, which is due to the energy power spectrum $\propto
E^{-3}$ of the cosmic ray flux that produce the neutrinos in the atmosphere
[28].

From equation (3), one can generate the two-dimensional (2D) distributions of
expected number of neutrino interactions of a given flavor as a function of
$(E,\theta_{z})$, or oscillograms, for any specific exposure of a perfect
detector that would observe atmospheric neutrinos from all directions. Because
of the dependence of neutrino oscillation on $n_{e}$ along the neutrino path,
changes in the outer core composition (i.e., different Z/A) induce variations
in the expected signal, depending on the neutrino energy and exact trajectory
(or equivalently, $\theta_{z}$), that get imprinted in the oscillograms.

We have quantified this effect by computing the relative difference
($\Delta\\!N/N$) in expected number of neutrino events between oscillograms
generated with different models of outer core composition. We considered five
different models of outer core composition (see table 1), ranging from a pure
FeNi alloy to a Hydrogen-rich (FeNiH ) one. figure 2 illustrates the expected
impact of the composition in terms of $\Delta\\!N/N$ for muon- and electron-
neutrino interactions, as a function of the neutrino energy, for two different
incoming directions. The predicted numbers of interactions for those
particular configurations are shown to differ by up to 4% for electron
neutrinos, and up to 20% for muon neutrinos, depending on the core composition
considered. Furthermore, because the shape of the signal changes a lot
depending on the zenith angle, it appears important to consider the signal in
two dimensions (i.e. as a function of $E$ and $\theta_{z}$). This is further
illustrated in figure 3 showing the full 2D distributions of $\Delta\\!N/N$
for muon- and electron-neutrino interactions, for a specific comparison
between FeNiH and FeNi core compositions. Such results confirm that a detailed
study of the neutrino event rate in the energy range 1 - 10 GeV has the
potential to constrain the Earth’s outer core composition.

###  2.2 Sensitivity of upcoming neutrino detectors to the outer core
composition

The distributions shown in the upper panels of figure 3 implicitly relate to a
detector with 100% detection efficiency, perfect energy and $\theta_{z}$
resolution, and infinite statistics. In reality, the measurement accuracy will
be limited by experimental effects, resulting in some blurring and attenuation
of the signal, as illustrated in figure 3 (lower panels). The detector
technology and specifications (size, detection efficiency, resolution and
particle identification) affect its ability to reconstruct the neutrino
properties (energy, direction and flavor). Additionally, the intrinsically
probabilistic nature of neutrino interactions as a quantum process induces a
statistical uncertainty on the final measurement (in terms of observed number
of events). The detector thus needs to be scalable to a sufficiently large
volume for this uncertainty to be smaller than the intrinsic signal.

Neutrinos are observed only indirectly, through the signal deposited in the
detector by the byproducts of their interaction. The ability to detect,
reconstruct an identify the different types of neutrino events in the energy
range of interest for Earth tomography is essentially driven by the detector
size, technology and density of sensors. We consider here two main
experimental approaches that are currently pursued for the upcoming generation
of neutrino detectors at the GeV scale, and that rely on different target
materials and observation strategies. Liquid Argon (LAr) detectors, such as
DUNE [29], are a type of Time Projection Chamber (TPC) [30] that detects the
electrons released upon ionization of Argon atoms along the paths of the
secondary charged particles emerging from the neutrino interaction. Water-
Cherenkov (wC) detectors instrument large volumes of water (or ice) with 2D –
such as Hyperkamiokande [31, 32] – or 3D – such as ORCA [33, 34]– arrays of
photosensors that detect the Cherenkov light induced by the secondary charged
particles traveling faster than light in water [35, 36]. As further discussed
in section 3.3, LAr-detectors reach the highest reconstruction precision,
while wC-detectors are scalable to larger volumes, especially when deployed in
a natural environment as is the case of ORCA.

To assess the impact of these experimental limitations on the sensitivity of
neutrino oscillation tomography, we have modeled the detector response with
four key features: the effective mass – defined as the product of the
instrumented target mass and the detection efficiency –, the direction and
energy resolutions – which characterize the mapping of the neutrino true
energy and zenith angle $(E,\theta_{z})$ into the respective reconstructed
quantities $E_{\text{reco}}$ and $\theta_{\text{reco}}$ –, and the efficiency
in identifying the neutrino flavor based on the secondary particles produced
in the interaction. In this study we assume that all detectors under
consideration have the capability to classify events into two main
observational classes: track-like events, which are mostly associated to
$\nu_{\mu}$ charged-current interactions producing a long muon track, and
cascade-like events, when there is no such single muon or the muon energy is
too small for it to be distinguished among the many other secondary particles
created in the event. The cascade channel therefore comprises mainly $\nu_{e}$
and (subdominantly) $\nu_{\tau}$ interactions, plus a small contribution from
$\nu_{\mu}$ interactions without an observable muon. More detail about the
mapping of neutrino interaction types into the observational track and cascade
classes is provided in section 3.3.

To investigate the influence of the detector characteristics on their
potential to constrain the core composition, we consider four benchmark
parametrizations: three with the expected capabilities of the detectors under
construction (DUNE, ORCA, and HyperKamiokande), and one additional,
hypothetical “Next-Generation" (NextGen) detector with improved performances.
The corresponding modelling, and the physical parameters used for each
specific detector, are described in section 3.3 and table 3. By folding the
detector response with the expected rate of neutrino interactions given by
equation (3), and summing over all relevant interaction types, we obtain the
predicted rate of observable events of a given observational class,
$R^{\text{obs}}_{\text{track}}(E_{\text{reco}},\theta_{\text{reco}})$ or
$R^{\text{obs}}_{\text{cascade}}(E_{\text{reco}},\theta_{\text{reco}})$ in a
given bin of neutrino reconstructed energy and zenith angle. By integrating
these quantities over time and effective mass of the detector, we generate 2D
histograms of the expected number of events as a function of $E_{\text{reco}}$
and $\theta_{\text{reco}}$, for a given detector exposure. Such oscillograms
are the cornerstone of our discussion of the detectors capability to constrain
the Earth’s core composition, as illustrated in the lower panels of figure 3.

For a given detector, our projected experimental sensitivities are obtained
from the detailed comparison of the full 2D oscillograms of expected neutrino
events, generated with different assumptions for the core composition. A
statistically meaningful way of quantifying the detector performance,
described in section 3.4, is to apply a $\chi^{2}$ hypothesis test to the 2D
histograms of expected events in bins of $(E_{reco},\theta_{reco})$. The total
$\Delta\chi^{2}$ associated to a given pair of composition models, say A and
B, which under the assumption of the validity of Wilk’s theorem is directly
related to the significance $\sigma=\sqrt{\Delta\chi^{2}}$, giving the
confidence level (C.L.) by which model A can be discriminated from model B.
Examples of signed $\Delta\chi^{2}$ maps for different upcoming detectors are
presented in figure 4 for the discrimination between FeNi and FeNiH models.
The three detectors achieve a comparable statistical significance in the
measurement of the outer core Z/A after 20 years of data taking, although
their sensitivity does not necessarily come from the same region of the
$(E,\theta)$ plane, nor from the same observational channel.

Figure 5 compares more directly the performances of the benchmark detectors in
discriminating outer core composition models. Despite the differences in
technology, size and reconstruction performances, we find that ORCA and DUNE
reach a similar precision of 0.016 (at 1$\sigma$ C.L.) on the absolute Z/A
measurement after 20 years of data taking. These results are comparable with
the ones previously published [24] based on a full simulation of the ORCA
detector. They confirm the capability of these detectors to measure the Z/A in
the Earth’s outer core with a 1$\sigma$ relative precision of 3 to 4 percent.
HyperKamiokande performs slightly better, with a relative precision of 2.5%,
yielding a discrimination power of approximately 0.5$\sigma$ between FeNi and
FeNiH models for 25 years of measurements. A higher sensitivity can be
achieved by combining the data from the three detectors into one single
measurement; in that case, a 1$\sigma$ C.L. discrimination power for FeNiH vs.
FeNi could be reached in 50 years of concomitant data taking, as can be seen
in figure. 6.

###  2.3 A next-generation detector to identify the light elements in the
outer core

While the previous result can be seen as a promising proof of concept of the
method, achieving the next level of sensitivity requires both beating the
statistical limitation by scaling up the detectors in size, and achieving
excellent reconstruction performances (flavor, $E_{\text{reco}}$,
$\theta_{\text{reco}}$), in order to better resolve the patterns in the 2D
$(E_{reco},\theta_{reco})$ event distributions.

Based on the results obtained for the upcoming detectors, we have investigated
possible realizations of such a next-generation, or NextGen, detector, that
would combine the best performances of current-generation instruments, taking
as benchmarks a 1 GeV detection threshold and almost perfect (99%)
discrimination between track-like and cascade-like events.

We performed a systematic scan of the discrimination potential of such
detector over a large range of sizes, energy and angular resolutions, as
presented in figure 7. A separation power of at least $1\sigma$ between
FeNiSi2O4 and FeNiH models can be obtained with a rather wide combination of
zenith/energy resolution parameters. A good compromise is found in the region
around $\sigma(\theta)=7^{\circ}$ angular resolution and $\sigma(E)/E=10\%$
energy resolution, not far from the performances associated with
HyperKamiokande, but for a significantly larger, 10-Mton, detector. We provide
a specific parametrization of such a NextGen detector in table 3.

This scan also illustrates how the small size and limited scalability of
Liquid Argon detectors is a serious obstacle to developing this technique to
achieve the required exposure levels in a reasonable running time, despite
their superior energy and zenith resolution. Even an asymptotic evolution of
the DUNE experiment, with perfect identification of all neutrino flavors and
interaction channels, would only reach a $0.5\text{\,}\sigma$ discrimination
power for FeNiSi2O4 vs. FeNiH after 50 years running. We conclude that better
perspectives for a core composition measurement arise from the water Cherenkov
approach, provided that event reconstruction and identification performances
similar or better to those of HyperKamiokande can be achieved in a much larger
(hence likely more sparsely instrumented) detector.

The performance of the NextGen detector is illustrated in figure. 8 in terms
of the signed $\Delta\chi^{2}$ maps for the discrimination between a FeNiH and
a FeNiSi2O4 composition of the outer core, for the track and cascade channels,
for 20 years of data taking. In this example, the total $\Delta\chi^{2}$ value
is 2.71, respectively 1.68 for track and 1.03 for cascade channel.

The qualitative jump in performance that could be achieved with the NextGen
detector is also evident in figure 6 where its sensitivity is compared to the
one obtained from combining ORCA, DUNE and Hyperkamiokande data. The
characteristics of the NextGen detector yield a sub-percent precision on the
Z/A measurement, sufficient to distinguish between FeNiSi2O4 and FeNiH
compositions at $>1\sigma$ in a running time of about 10 years, well within
the typical lifetime of a neutrino experiment (a discrimination at 90% C.L.
would be achieved in about 20 years).

We further illustrate such capabilities in figure. 9, where we address the
separability (at the $1\sigma$ level) between pairs of realistic core
composition models, as a function of the detector running time. Although a
discrimination between composition models whose Z/A are closer than 0.001
(i.e. FeNiSi7O2 vs FeNiSi2O4 , or FeNiSi2O4 vs. FeNiSiH ) appears out of
reach, our results suggest that the NextGen detector would be able to exclude
or confirm the FeNiH hypothesis against all the other realistic models
considered in this study. This result can be obtained in less than 20 years if
the actual outer core composition is in the family of models with an admixture
of Silicium and Oxygen as only light elements.

##  3 Materials and methods

###  3.1 Neutrino oscillations in matter

Neutrino flavor oscillations are a quantum-mechanical phenomenon that arises
because the neutrino flavor eigenstates $\nu_{e},\nu_{\mu},\nu_{\tau}$ – which
take part in weak interactions and are therefore the observable states – are
not identical to the neutrino mass eigenstates $\nu_{1},\nu_{2},\nu_{3}$ –
which describe their propagation in vacuum. The flavor eigenstates are a
quantum superposition of the mass eigenstates, whose relative phases change
along the neutrino propagation path. This evolving mix of states leads to an
oscillatory pattern of the detection probability of a given neutrino flavor,
which depends on the neutrino energy and traveled distance.

Neutrino oscillation probabilities are calculated by solving the Schrödinger
equation $i\partial_{t}|\nu(t)\rangle=H|\nu(t)\rangle$. The neutrino
propagation Hamiltonian $H$ is approximated as the sum of two terms:

$H=U\left(\begin{array}[]{ccc}0&0&0\\\ 0&\frac{\Delta{m^{2}_{21}}}{2E}&0\\\
0&0&\frac{\Delta{m^{2}_{31}}}{2E}\end{array}\right)U^{\dagger}+\left(\begin{array}[]{ccc}V_{e}&0&0\\\
0&0&0\\\ 0&0&0\end{array}\right).$ (4)

The first term represents intrinsic energy levels of the system in vacuum.
Because neutrinos are generally observed in an ultra-relativistic context, the
energy difference between mass states is approximately proportional to their
difference in mass-squared $\Delta m_{ij}^{2}=(m_{i}^{2}-m_{j}^{2})$. The
Hamiltonian is shown in the flavor basis characterized by eigenstates of the
weak interaction. The mixing matrix $U$ represents the unitary transformation
between the mass and flavor bases and is usually parametrized in terms of
three mixing angles ($\theta_{12}$, $\theta_{13}$, $\theta_{23}$) and a
complex phase $\delta$, which are fundamental constants of physics (see e.g.
[37] for a general discussion of oscillation formulae).

The second term in equation (4) arises from coherent interactions of neutrinos
with electrons in the medium in which they propagate [12, 13]. The effective
potential $V_{e}=\pm\sqrt{2}G_{F}n_{e}$ induces a change in the electron
flavor energy level that is directly proportional to the electron number
density $n_{e}$ and the Fermi constant $G_{F}$, which characterizes the
strength of the weak interaction. The sign of the matter potential is positive
for neutrinos and negative for antineutrinos.

The solution of the Schrödinger equation for such a system in matter of
constant density amounts to computing the eigensystem of the full Hamiltonian.
While in vacuum the eigenvalues and eigenvectors are given by
$\Delta{m^{2}_{ij}}/2E$ and the matrix $U$ directly, the matter potential
leads to effective eigenvalues and eigenvectors that depend on the electron
density. This process can be interpreted as a modification of the fundamental
parameters $\Delta{m^{2}_{ij}}$ and $\theta_{ij}$ induced by the neutrino
environment.

For example, the effective mixing angle $\tilde{\theta}_{13}$ and the
effective mass-squared splitting $\Delta\tilde{m}_{31}^{2}$ in matter are
related to the corresponding parameters in vacuum via:

$\displaystyle\Delta\tilde{m}^{2}_{31}\approx\xi\cdot\Delta
m^{2}_{31},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\sin^{2}2\tilde{\theta}_{13}\approx\frac{\sin^{2}2\theta_{13}}{\xi^{2}},$ (5)

with a mapping parameter

$\displaystyle\xi$ $\displaystyle=\sqrt{\sin^{2}2\theta_{13}+\left(\cos
2\theta_{13}-\frac{2\,E\,V_{e}}{\Delta m_{31}^{2}}\right)^{2}}.$ (6)

At the energy and distance scales relevant to atmospheric neutrino
oscillations, the probability of observing a transition between $\nu_{e}$ and
$\nu_{\mu}$ flavors is proportional to $\sin^{2}2\theta_{13}$, which is small
in vacuum, as observed in oscillation measurements of reactor antineutrinos
[38, 39, 40] that are insensitive to the matter potential. However, the
$\nu_{e}\leftrightarrow\nu_{\mu}$ transition probability may become large or
even maximal when $\xi^{2}\rightarrow\sin^{2}{2\theta_{13}}$, which translates
into a resonance condition for the neutrino energy:

$E\rightarrow\pm\frac{\Delta
m_{31}^{2}\cos{2\theta_{13}}}{2\sqrt{2}\,G_{F}\,n_{e}}.$ (7)

Neutrinos traversing the Earth will experience resonant oscillations in the
outer core for energies $\sim$$3\text{\,}\mathrm{G}$eV, and in the mantle for
$\sim$$7\text{\,}\mathrm{G}$eV, with the exact resonant energy depending on
the electron density of the material.

###  3.2 Computation of atmospheric neutrino propagation through the Earth and
determination of the interacting event rate

An almost isotropic neutrino flux is constantly produced by the interaction of
cosmic rays with the atmosphere, mainly coming from the decays of charged
mesons ($\pi^{\pm}$’s and $K^{\pm}$’s) and muons. For the present study, we
use as input the flux calculation produced by [28] for the Gran Sasso site
(without mountain over the detector), averaged over all azimuth angles and
assuming minimum solar activity. The flux is dominated by muon- and
electron-(anti)neutrino components, with an approximate ratio of 2:1 between
$\nu_{\mu}$ and $\nu_{e}$ flavors, and we have neglected the small
contribution of $\nu_{\tau}/\bar{\nu}_{\tau}$.

For a given zenith angle of incidence $\theta_{z}$, as defined in figure 1a,
the neutrino trajectory across the Earth is modelled along the corresponding
baseline through a sequence of steps of constant electron density which are
inferred from a radial model of the Earth with 42 concentric shells of
constant $n_{e}$, where mass density values are fixed and follow PREM. These
shells are grouped into three petrological layers (inner core, outer core, and
mantle+crust), whose composition, hence Z/A factor, is assumed to be uniform
and provided in table 2). In each shell, the electron density is determined
from the mass density and Z/A according to equation (2). The probabilities of
neutrino flavor transitions along their path through the Earth are computed
using the OscProb 111J. Coelho et al., https://github.com/joaoabcoelho/OscProb
package. The values of the parameters that enter the oscillation probability
computation are taken from the global fit of neutrino data performed with the
NuFIT analysis Version 5.0 [41]. We have assumed normal ordering of the
neutrino mass states, i.e. $m_{3}>m_{1}$, which is currently favoured by
global fits of neutrino data.

Neutrinos with energies in the range $\sim 1$-$100\text{\,}\mathrm{G}$eV
interact with matter mostly via scattering off nuclei, by exchanging either a
neutral ($Z^{0}$) or a charged ($W^{\pm}$) weak-force boson that triggers a
hadronic cascade. In neutral-current (NC) interactions, the neutrino survives
in the scattering products and escapes. In charged-current (CC) interactions,
the neutrino gives rise to a charged lepton counterpart of the same flavor
($e$, $\mu$, or $\tau$). While electrons immediately induce a secondary short
electromagnetic cascade, muons travel relatively long distances, depending on
their energy, before they are stopped. Because of the high mass of the $\tau$
lepton ($\approx$ $1.7\text{\,}\mathrm{GeV}$), the contribution of
(anti-)$\nu_{\tau}$-induced CC events is small in the range of energies
relevant for tomography studies. For completeness, this flavor is nevertheless
included in our computations.

The rate of neutrinos interacting in a given volume of target matter is then
computed, for each interaction channel (NC/CC; $e,\mu,\tau$;$\nu/\bar{\nu}$)
according to equation 3, using neutrino-nucleon cross-sections weighted for
water molecules, obtained with the GENIE Monte Carlo neutrino generator [42].
We have neglected here the small difference in cross-section for interactions
on Argon nuclei in the case of the DUNE detector.

###  3.3 Detector modelling

A neutrino interaction occurring within the detection volume will generate
secondary signals that constitute the event recorded by the detector. The
ability to detect, reconstruct an identify the different types of neutrino
events in the energy range of interest for Earth tomography is essentially
driven by the detector size, technology and density of sensors.

The performance of water-Cherenkov detectors is typically a trade-off between
the target volume for neutrino interactions and the density of photosensors
that sample the Cherenkov signal emitted by charged byproducts of the neutrino
interaction. The HyperKamiokande experiment [31, 32] will consist in two $\sim
200$-kton water tanks overlooked by photosensors covering its internal walls,
providing a sub-GeV threshold for neutrino detection and excellent
discrimination between electron-like and muon-like signatures. To access even
larger target volumes, the ORCA [33, 34] and PINGU [43] experiments propose to
instrument several Mtons of seawater or polar ice with a much sparser 3D array
of photosensors. The gain in statistics comes at the expense of a higher
detection threshold (few GeV) and worse energy and direction reconstruction
capabilities, which also limit their performances for neutrino flavor
identification. Water-Cherenkov detectors also cannot distinguish in first
approximation between neutrino and anti-neutrino events.

Liquid Argon Time Projection Chambers, on the other hand, are able to
reconstruct highly detailed 3D images of the neutrino event, providing
excellent flavor identification capabilities, even at low (sub-GeV) energies,
and energy and angular resolutions typically superior to water-Cherenkov
detectors. The size of these detectors is however limited by their much higher
cost. The first large-scale LArTPC detector was the ICARUS T600 detector [44],
with an active mass of 476 tons. A new generation of LArTPCs is being
developed for the DUNE experiment [29], that will include 4 detectors of 10
kton active mass each.

To simulate the detector response, we use a set of parametrized analytical
functions modeling the main performance features relevant for the detection,
reconstruction and classification of neutrino events:

  * •

The effective mass $M_{\mathrm{eff}}(E)$ is the product of the instrumented
target mass $M$ of the detector and its detection efficiency, i.e. the
probability for a neutrino interaction to be successfully detected as an
event. $M_{\mathrm{eff}}$ typically increases with the neutrino energy, until
reaching a plateau that saturates the fiducial mass of the detector. We have
conservatively neglected here a potential increase in $M_{\mathrm{eff}}$ for
$\nu_{\mu}$ CC events at high energies, corresponding to through-going muons
created in neutrino interactions outside of the detector target volume. The
threshold for detection is mainly driven by the density (and intrinsic
efficiency) of sensors. We approximate $M_{eff}(E)$ by a sigmoid function of
$\log(E)$ with two adjustable parameters $E_{th}$ and $E_{pl}$, which
correspond to energies where the detection efficiency reaches 5% and 95%
respectively.

  * •

The energy resolution parameterizes the relative error on the reconstructed
neutrino energy in the form of a Gaussian probability distribution function
(p.d.f.) with energy-dependent width: $\sigma_{E}(E)/E=A_{E}+B_{E}/\sqrt{E}$.

  * •

The angular resolution parameterizes the error on the measured zenith angle in
the form of a von Mises-Fisher p.d.f. on a sphere marginalized with respect to
azimuth. For wC detectors we take into account the dependence of
$\sigma_{\theta}$ on the energy as
$\sigma_{\theta}(E)=A_{\theta}+B_{\theta}\sqrt{E}$.

  * •

The classification efficiency $\varepsilon_{class}(E)$ describes the
probability for a neutrino event of energy $E$ to be correctly classified into
one of the topological channels observable by the detector. We model it with a
sigmoid as a function of $\log(E)$ with adjustable threshold
($E^{class}_{th}$) and plateau ($E^{class}_{pl}$) energies, maximal
identification probability $P^{class}_{max}$ and minimum probability of 50%
(corresponding to the case of no separation power).

Events are classified according to their topological features, which depend on
the signature of each interaction channel (NC/CC; $e/\mu/\tau;\nu/\bar{\nu}$).
In practice, the relative sparseness of instrumentation limits the signal
sampling in the detector, and thus the reconstruction and identification
performances. All experiments discussed here have a basic capability to
classify events into two main observational classes: tracks, when a long muon
track is observed, most of the time originating from a $\nu_{\mu}$ CC event;
and cascades, the latter comprising $\nu_{e}$ CC events, but also most
$\nu_{\tau}$ CC interactions and all NC interactions. For NC interactions, a
fraction of the energy of the event is carried away by an invisible outgoing
neutrino, resulting in a lower reconstructed energy and degraded angular
resolution. Because NC-induced events are blind to neutrino flavor, they only
decrease the experiment sensitivity. Fine-grained detectors like DUNE and
HyperKamiokande have some capability to separate and reject the NC-induced
events from the $\nu_{e}$ CC events, thereby increasing the flavor purity of
the cascade channel. To account for a possible related improvement in the
detector sensitivity, we have considered here both the baseline case that
includes the NC-induced events in the cascade sample, and the optimistic case
where the NC-induced contributions are completely removed from the sample.

Table 3 provides the specific inputs chosen for different benchmark detectors
representative of the upcoming generation of GeV neutrino detectors:
HyperKamiokande, ORCA, and DUNE. The DUNE-like model accounts for the superior
resolution and reconstruction capabilities of the Liquid Argon detection
technique, at the cost of a much smaller instrumented volume. At the other
extreme, the ORCA-like model reflects the higher detection threshold and
degraded reconstruction performances that result from a sparsely instrumented
- although larger - volume. The HyperKamiokande-like model is a middle-way
option combining a medium-size detector with a low detection threshold and
good reconstruction and identification capabilities. Some examples of the
corresponding response functions are presented in figure 10.

###  3.4 Computation of the expected signal

The rate of interacting neutrinos is calculated as in equation (3); in the
following we use the subscript “true" for the quantities (energy and zenith
angle) that refer to the incoming neutrino properties. This rate is then
weighted with the appropriate detection efficiency at each incoming neutrino
energy $E_{\text{true}}$. The expected signal in a realistic experiment is
computed by a discretized convolution of the interacting rate at energy
$E_{\text{true}}$ and incident at zenith angle $\theta_{\text{true}}$ over the
E-$\cos{\theta}$ plane according to the energy and zenith resolution p.d.f.s.
The reconstructed events for each interaction channel are distributed into the
two observational channels (tracks and cascades) according to the
classification efficiency function $\varepsilon_{\text{class}}({E})$. Every
bin in the final, $(E_{\text{reco}},\theta_{\text{reco}})$ event oscillogram
therefore (i) contains a certain fraction of misreconstructed events coming
from other $(E_{\text{true}},\theta_{\text{true}})$ bins (ii) misses some
events that end up misreconstructed into different
$(E_{\text{reco}},\theta_{\text{reco}})$ bins.

The final event rate expected in a given channel, in a given bin of
reconstructed energy and zenith angle, is thus obtained as follows for the
track channel:

$\displaystyle
R^{\text{obs}}_{\text{tracks}}(E_{\text{reco}},\theta_{\text{reco}})\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\sum_{E_{\text{true}},\theta_{\text{true}}}\left[R^{\text{int}}_{\text{tracks}}(E_{\text{true}},\theta_{\text{true}})\,\varepsilon_{\text{class}}(E_{\text{true}})\
+R^{\text{int}}_{\text{casc}}(E_{\text{true}},\theta_{\text{true}})\,\left(1-\varepsilon_{\text{class}}(E_{\text{true}})\right)\right]$
(8) $\displaystyle\hskip 45.52458pt\times
PDF_{\text{angle}}(\theta_{\text{reco}};E_{\text{true}},\theta_{\text{true}})\times
PDF_{\text{energy}}(E_{\text{reco}};E_{\text{true}})$ $\displaystyle\hskip
45.52458pt\times\Delta E_{\text{true}}\times\Delta\theta_{\text{true}}\times
M_{\mathrm{eff}}(E_{\text{true}})$

where $PDF_{\text{angle}}$ and $PDF_{\text{energy}}$ represent the probability
distribution functions of the angular and energy resolutions as defined above.
A similar expression can be straightforwardly derived for the event rate in
the cascade channel.

Independently of any classical measurement error, the intrinsically
probabilistic nature of neutrino interactions as a quantum process induces a
statistical uncertainty on the number of events observed in each bin of
$(E_{\text{reco}},\theta_{\text{reco}})$. The probability to detect $N$
neutrino events of a given topology, e.g. tracks, in a predetermined interval
of reconstructed energy and zenith angle is distributed according to a Poisson
law, parametrized by the expected number of events
$N_{\text{exp}}=R^{\text{obs}}_{\text{tracks}}(E_{\text{reco}},\theta_{\text{reco}})\times
T$, where $T$ corresponds to the duration of data taking. The observed number
$N$ will randomly deviate from this expectation with a relative standard
deviation $\sigma(N)/N_{\text{exp}}=1/\sqrt{N_{\text{exp}}}$. This uncertainty
on the final measurement can only be mitigated by increasing the size of the
event sample, which is proportional to the detector exposure.

Once these fluctuations are taken into account, we compute the log-likelihood
ratio $\lambda_{LLR}$ of the resulting event histograms, which yields a
measure of the significance to reject model B (the "test hypothesis") if model
A (the "null hypothesis") is true. For large enough statistics, the Poisson-
distributed events converge towards a normal distribution, so that
$-2\lambda_{LLR}$ can be approximated by:

$-2\lambda_{LLR}\approx\Delta\chi^{2}_{i}=\frac{(n_{A}^{i}-n_{B}^{i})^{2}}{n_{A}^{i}}$
(9)

The total $\Delta\chi^{2}$ (integrated over all bins of energy and zenith
angle) is related to the significance $\sigma=\sqrt{\Delta\chi^{2}}$, that
gives the confidence level (C.L.) by which the test hypothesis can be excluded
under the assumption of the null hypothesis. For illustration purposes, we
present the maps in terms of the signed $\Delta\chi^{2}$, ie. by mutiplying
the value of $\Delta\chi^{2}_{i}$ in each bin by $sign(n_{A}^{i}-n_{B}^{i})$.

##  4 Discussion

This paper presents a detailed study and comparison of the potential of
upcoming atmospheric neutrino detectors to constrain the outer core
composition to precision levels of relevance for geophysics. By computing the
propagation and oscillation of atmospheric neutrinos through the Earth with
well-established tools from the neutrino physics community, we first confirm
the theoretical capability of neutrino oscillation tomography to discriminate
between different realistic core composition models. We note that the largest
contribution to the theoretical signal is provided by atmospheric neutrinos
reaching the detector in the $\nu_{\mu}$ flavor, which may explain why many of
the neutrino oscillation tomography studies to date have focused on this
detection channel only.

The ultimate reach of this method, however, significantly depends on how the
expected signal is effectively extracted in a realistic experiment. In order
to address this question in a flexible way and with a view to the future, we
have developed our analysis based on analytical parametrisations of the
detector response functions, similar in concept to the studies by Winter [20]
and Rott et al. [21]. This approach has allowed us to investigate for the
first time different classes of experiments - namely Liquid Argon TPCs and
water-Cherenkov detectors - within a unified framework. Its implementation in
terms of generic detection parameters would easily allow for further
comparisons with other neutrino detection techniques that could be proposed in
the future.

One notable conclusion to be drawn from our study, is that the respective
contribution of the different neutrino flavours reaching the detector to the
total signal (in terms of discrimination power for different core composition
models) can vary significantly for different choices of the parametrized
response functions, hence for different detectors, even within the same
experimental approach. This observation emphasizes the importance of
considering the neutrino signal in all its components - here, namely the two
observational channels tracks and cascades, in order to accurately assess the
expected reach of a given experiment.

Taking advantage of the versatility of our framework, we have further
investigated the desirable detection performances for a NextGen detector that
- contrarily to the upcoming generation of experiments - would be optimized
for neutrino tomography. We find that such performances stem from two
complementary aspects: (i) the large statistics of events provided by an
experiment with 10 Mton target mass and a relatively low ($\sim$ 1 GeV) energy
threshold that fully covers the energy range where resonant oscillation
effects are expected for core-crossing neutrinos; and (ii) the state-of-the-
art event reconstruction and classification capabilities. Limitations to this
measurement arise from the intrinsic, quantum fluctuations in the particle
content of the hadronic shower that results from the fragmentation of the
nucleus hit by a neutrino. These effects have been extensively discussed for
ORCA-like detectors [45] and shown to limit the performances as long as
individual particles in the neutrino-induced cascades cannot be reconstructed.
This limitation is mitigated in instrumented water tanks à la HyperKamiokande,
thanks to the high density of sensors that allows a sufficient sampling of the
particles produced by the neutrino interaction.

In view of the above considerations, our favours go to the water-Cherenkov
technique to define the path towards a Next-Generation detector able to
reliably test the FeNiH hypothesis for the earth’s outer core composition.
Such a detector should combine a multi-Megaton target volume with a sufficient
density of sensors to achieve an efficient sampling of the neutrino event
topology. Because of the limited scalability of man-made water tanks, we
rather advocate for a hyperdense network of light sensors deployed in natural
reservoirs of water or ice. The detector topology could be either three-
dimensional (à la ORCA) or two-dimensional (à la Hyperkamiokande) as long as
it provides good angular coverage for neutrinos crossing the Earth’s outer
core. In this latter case, one possible configuration could consist in an
immersed “carpet" of densely packed photosensors looking downwards, monitoring
a large target volume of water or ice. The resulting loss in sensitivity to
near-horizontal neutrinos (i.e. $\theta\sim 90^{\circ}$) should have minimal
impact on the tomography measurement: all neutrinos crossing the core will
reach the detector at an angle larger than $57^{\circ}$ (corresponding to
$\theta_{Z}\gtrsim 147^{\circ}$, and also a significant fraction of those
traversing only the mantle will be detected and exploitable as a reference
sample. Existing proposals for neutrino detectors made in different contexts,
such as the hyper-dense 3D network Super-ORCA [46] or the multi-megaton
shallow-water network of steel tanks TITAND proposed for proton decay studies
[47] might also be worth reevaluating in terms of their capabilities for
neutrino tomography.

Finally, we emphasize that the timescale required for the measurements
described here is inversely proportional to the total detection volume
available. Deploying multiple detector units at different positions on the
surface of the globe - an effort currently being achieved with sparse wC
detectors focusing on high-energy neutrino astronomy - has the extra appeal to
allow for a fully 3-dimensional analysis with potentially higher scientific
reach. Such a network of detectors would not only be of great interest for
earth tomography purposes, but would also provide an unprecedented tool for
high-statistics and high-precision studies of atmospheric neutrinos.

## Acknowledgments

This research has benefited from the financial support of the CNRS
Interdisciplinary Program 80PRIME, the LabEx UnivEarthS at Université de Paris
(ANR-10-LABX-0023 and ANR-18-IDEX-0001) and the Institut Universitaire de
France.

## References

  * [1] F. Birch.  Composition of the Earth’s Mantle.  Geophysical Journal International, 4(Supplement_1):295–311, 1961\.
  * [2] J.-P. Poirier.  Light elements in the Earth’s outer core: A critical review.  Phys. Earth Planet. Inter., 85:319–337, 1994.
  * [3] S. Hirose, K. Labrosse et al.  Composition and State of the Core.  Annu. Rev. Earth Planet Sci., 41:657–691, 2013.
  * [4] W. F. McDonough.  Compositional Model for the Earth’s Core, volume 2, pages 547–566.  Elsevier, 2003.
  * [5] A. M. Dziewonski et al.  Preliminary reference Earth model.  Phys. Earth Planet. Int., 25 (4):297–356, 1981.
  * [6] J. Badro et al.  Core formation and core composition from coupled geochemical and geophysical constraints.  PNAS, 112 (40):12310–12314, 2015.
  * [7] E. Kaminski et al.  A two-stage scenario for the formation of the Earth’s mantle and core.  Earth and Planetary Science Letters, 365:97–107, 2013.
  * [8] K. Umemoto et al.  Liquid iron-hydrogen alloys at outer core conditions by first-principles calculations.  Geophys. Res. Lett., 42(18):7513–7520, 2015.
  * [9] T. Sakamaki et al.  Constraints on Earth’s inner core composition inferred from measurements of the sound velocity of hcp-iron in extreme conditions.  Science Advances, 2(2):e1500802, 2016.
  * [10] L. Yuan et al.  Strong Sequestration of Hydrogen Into the Earth’s Core During Planetary Differentiation.  Geophysical Research Letters, 47:e2020GL088303, 2020.
  * [11] Y. Fukuda et al. (Super-Kamiokande Collaboration).  Evidence for oscillation of atmospheric neutrinos. Phys. Rev. Lett., 81:1562–1567, 1998 [hep-ex/9807003].
  * [12] S. Mikheyev et al.  Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos.  Sov. J. Nucl. Phys., 42:913–917, 1985.
  * [13] L. Wolfenstein.  Neutrino Oscillations in Matter. Phys. Rev. D, 17:2369–2374, 1978.
  * [14] V. K. Ermilova et al.  Buildup of Neutrino Oscillations in the Earth.  JETP Lett., 43:453–456, 1986.
  * [15] A. Nicolaidis.  Neutrinos for Geophysics. Phys. Lett. B, 200:553–559, 1988.
  * [16] A. Nicolaidis et al.  Neutrino tomography of the earth. J. Geophys. Res., 96(B13):21811–21817, 1991.
  * [17] T. Ohlsson et al.  Reconstruction of the earth’s matter density profile using a single neutrino baseline. Phys. Lett. B, 512:357–364, 2001 [hep-ph/0105293].
  * [18] T. Ohlsson et al.  Could one find petroleum using neutrino oscillations in matter?. Europhys. Lett., 60:34–39, 2002 [hep-ph/0111247].
  * [19] M. Lindner et al.  Tomography of the earth’s core using supernova neutrinos. Astropart. Phys., 19:755–770, 2003 [hep-ph/0207238].
  * [20] W. Winter.  Neutrino tomography: Learning about the earth’s interior using the propagation of neutrinos. Earth Moon Planets, 99:285–307, 2006 [physics/0602049].  [,285(2006)].
  * [21] C. Rott et al.  Spectrometry of the Earth using Neutrino Oscillations. Sci. Rep., 5:15225, 2015 [1502.04930].
  * [22] W. Winter.  Atmospheric Neutrino Oscillations for Earth Tomography. Nucl. Phys., B908:250–267, 2016 [1511.05154].
  * [23] S. Bourret et al. (KM3NeT Collaboration).  Neutrino oscillation tomography of the Earth with KM3NeT-ORCA. PoS, ICRC2017:1020, 2018.
  * [24] S. Bourret et al. (KM3NeT Collaboration).  Earth tomography with neutrinos in KM3NeT-ORCA. EPJ Web Conf., 207:04008, 2019.
  * [25] J. C. D’Olivo et al.  Earth tomography with atmospheric neutrino oscillations. Eur. Phys. J. C, 80(10):1001, 2020.
  * [26] A. Kumar et al.  Validating the Earth’s core using atmospheric neutrinos with ICAL at INO. JHEP, 08:139, 2021 [2104.11740].
  * [27] K. J. Kelly et al.  DUNE atmospheric neutrinos: Earth tomography. JHEP, 05:187, 2022 [2110.00003].
  * [28] M. Honda et al.  Atmospheric neutrino flux calculation using the NRLMSISE-00 atmospheric model. Phys. Rev. D, 92(2):023004, 2015 [1502.03916].
  * [29] B. Abi et al. (DUNE Collaboration).  Volume I. Introduction to DUNE. JINST, 15(08):T08008, 2020 [2002.02967].
  * [30] C. Rubbia.  The Liquid Argon Time Projection Chamber: A New Concept for Neutrino Detectors.  1977\.
  * [31] K. Abe et al. (Hyper-Kamiokande Proto- Collaboration).  Physics potential of a long-baseline neutrino oscillation experiment using a J-PARC neutrino beam and Hyper-Kamiokande. PTEP, 2015:053C02, 2015 [1502.05199].
  * [32] K. Abe et al. (Hyper-Kamiokande Collaboration).  Hyper-Kamiokande Design Report. arXiv:1805.04163, 2018.
  * [33] S. Adrian-Martinez et al. (KM3Net Collaboration).  Letter of intent for KM3NeT 2.0. J. Phys. G, 43(8):084001, 2016 [1601.07459].
  * [34] S. Aiello et al. (KM3NeT Collaboration).  Determining the neutrino mass ordering and oscillation parameters with KM3NeT/ORCA. Eur. Phys. J. C, 82(1):26, 2022 [2103.09885].
  * [35] P. A. Čerenkov.  Visible Radiation Produced by Electrons Moving in a Medium with Velocities Exceeding that of Light. Physical Review, 52(4):378–379, 1937.
  * [36] I. E. Tamm et al.  Visible Radiation Produced by Electrons Moving in a Medium with Velocities Exceeding that of Light. Dokl. Acad. Sci. USSR, 14(3):107–112, 1937.
  * [37] C. Giganti et al.  Neutrino oscillations: the rise of the PMNS paradigm. Prog. Part. Nucl. Phys., 98:1–54, 2018 [1710.00715].
  * [38] Y. Abe et al.  Indication of reactor $\bar{\nu}_{e}$ disappearance in the double chooz experiment. Physical Review Letters, 108(13), 2012.
  * [39] F. P. An et al.  Observation of electron-antineutrino disappearance at daya bay. Physical Review Letters, 108(17), 2012.
  * [40] J. K. Ahn et al.  Observation of reactor electron antineutrinos disappearance in the reno experiment. Physical Review Letters, 108(19), 2012.
  * [41] I. Esteban et al.  The fate of hints: updated global analysis of three-flavor neutrino oscillations. JHEP, 09:178, 2020 [2007.14792].
  * [42] C. Andreopoulos et al.  The GENIE Neutrino Monte Carlo Generator. Nucl. Instrum. Meth. A, 614:87–104, 2010 [0905.2517].
  * [43] M. G. Aartsen et al. (IceCube Collaboration).  PINGU: A Vision for Neutrino and Particle Physics at the South Pole. J. Phys., G44(5):054006, 2017 [1607.02671].
  * [44] C. Rubbia et al.  Underground operation of the ICARUS T600 LAr-TPC: first results. JINST, 6:P07011, 2011 [1106.0975].
  * [45] S. Adrián-Martínez et al.  Intrinsic limits on resolutions in muon- and electron-neutrino charged-current events in the KM3NeT/ORCA detector. JHEP, 05:008, 2017 [1612.05621].
  * [46] J. Hofestädt et al.  Super-ORCA: Measuring the leptonic CP-phase with Atmospheric Neutrinos and Beam Neutrinos. PoS, ICRC2019:911, 2020 [1907.12983].
  * [47] Y. Suzuki et al. (TITAND Working Group Collaboration).  Multimegaton water Cherenkov detector for a proton decay search: TITAND (former name TITANIC).  In 11th International School on Particles and Cosmology, pages 288–296, 2001.
  * [48] S. Tagawa et al.  Compression of Fe–Si–H alloys to core pressures.  Geophysical Research Letters, 43(8):3686–3692, 2016.

## Tables

Label | FeNi |  FeNiSi2O4 |  FeNiSi7O2 | FeNiSiH | FeNiH
---|---|---|---|---|---
Composition | 95 wt% Fe | 94 wt% Fe | 91 wt% Fe | 93.2 wt% Fe | 94 wt% Fe
| 5 wt% Ni | 5 wt% Ni | 5 wt% Ni | 5 wt% Ni | 5 wt% Ni
| - | 2 wt% Si | 7 wt% Si | 6.5 wt% Si | 1 wt% H
| - | 4 wt% O | 2 wt% O | 0.3 wt% H | -
Z/A | 0.4661 | 0.4682 | 0.4691 | 0.4699 | 0.4714
Table 1: Models of outer core compositions considered in this study, showing the weight fraction of the different elements and corresponding average $Z/A$. The FeNi model is the benchmark alloy used in the present study to model the inner core composition, whereas other models introduce different combinations of light elements: FeNiSi2O4 [6], FeNiSi7O2 [7], FeNiSiH [48] and FeNiH [9]. In all models the Ni content is set to 5 wt% and Fe is the complement to 100% once light elements have been taken into account. All elements in the table have a $Z/A$ ranging between 0.46 and 0.50, except Hydrogen whose $Z/A=1$ is responsible for pulling the average $Z/A$ to higher values. Layer | Shells |  [$R_{-},R_{+}$] | Z/A
---|---|---|---
Inner core | 7 | 0 - 1221.5 | 0.466
Outer core | 12 | 1221.5 - 3480.0 | 0.466
Mantle + crust | 23 | 3480.0 - 6368.0 | 0.496
Table 2: Compositional layers in the benchmark Earth model used in the analysis when assuming a pure FeNi outer core. The columns indicate the number of constant density shells (each of them taking its value from PREM), the exact innermost and outermost radius (in km) and the assumed Z/A value. Detector |  $M$ (Mton) |  $E_{th}$ (GeV) |  $E_{pl}$ (GeV) | $\sigma(E)/E$ |  $\sigma_{\theta}$ (deg) |  $E_{th}^{class}$ (GeV) |  $E_{pl}^{class}$ (GeV) | $P_{max}^{class}$
---|---|---|---|---|---|---|---|---
ORCA-like | 8 | 2 | 10 | 25% | $30/\sqrt{E}$ | 2 | 10 | 85%
HyperKamiokande-like | 0.40 | 0.1 | 0.2 | 15% | $15/\sqrt{E}$ | 0.1 | 0.2 | 99%
DUNE-like | 0.04 | 0.1 | 0.2 | 5% | 5 | 0.1 | 0.2 | 99%
Next-Generation | 10 | 0.5 | 1.0 | $5\%+10\%/\sqrt{E}$ | $2+10/\sqrt{E}$ | 0.5 | 1 | 99%
Table 3: Inputs for the response functions of the detectors considered in this
study: total target mass; threshold and plateau energy for the detection
efficiency curve; energy and zenith resolutions; threshold and plateau energy
for the classification efficiency curve; maximal classification probability
achievable.

## Figures

Figure 1: Oscillations of Earth-crossing neutrinos. Left: A detector located
at the surface will receive atmospheric neutrinos having traversed the Earth
along different paths, as measured by their zenith angle $\theta_{z}$ with
respect to the vertical at the detector location. As these neutrinos cross the
Earth, the amplitude of their oscillations at different depths will be
affected by the local density of electrons. A modification in the electron
density of the outer core would then affect the flavor mix of neutrinos that
cross it. Right: To illustrate this effect, the number of interactions
expected for atmospheric muon neutrinos and antineutrinos
($\nu_{\mu}+\bar{\nu}_{\mu}$) is shown as a function of the neutrino energy in
a detector with 200 Mton-years exposure, assuming two extreme chemical
compositions chosen for the sake of the argument, corresponding respectively
to pure FeNi (a light-element-free outer core) and to pyrolite (same as the
mantle). Only neutrinos expected from the angular region
$\theta_{z}\in[160^{\circ},165^{\circ}]$ highlighted on Panel (a), with a
large intersection with the outer core, are shown. The observed oscillation
pattern reflects a combination of disappearing
($\nu_{\mu}\rightarrow\nu_{e},\nu_{\tau}$) and appearing
($\nu_{e}\rightarrow\nu_{\mu}$) muon-neutrinos. Differences related to the
outer core composition model are most visible for neutrino energies between 2
and 5 GeV.

Figure 2:  Predicted relative difference $\Delta\\!N/N$ in the number of
neutrino interactions as a function of the incoming (or “true") neutrino
energy, for the different outer core compositions described in Tab. 1, taking
FeNi as a reference. Upper panels are for $\nu_{\mu}+\bar{\nu}_{\mu}$, while
lower panels are for $\nu_{e}+\bar{\nu}_{e}$. Two incoming neutrino directions
have been considered: $\theta_{z}=165^{\circ}$ (left) and $160^{\circ}$
(right).

Figure 3: Upper panels: expected relative difference $\Delta\\!N/N$ in the
number of neutrino interactions, as a function of the incoming (or “true")
neutrino energy and zenith angle, for the composition models FeNiH vs. FeNi .
Left panel is for $\nu_{\mu}+\bar{\nu}_{\mu}$, while right panel is for
$\nu_{e}+\bar{\nu}_{e}$. The horizontal dashed lines correspond to the two
incoming directions considered for the plots in figure 2. Lower panels:
expected relative difference $\Delta\\!N/N$ in the number of detected events,
as a function of the reconstructed neutrino energy and zenith angle, for the
composition models FeNiH vs. FeNi , for the Next-Generation detector described
in section 3.3 and table 3). Left panel is for the track-like events (mostly
associated to $\nu_{\mu}$ interactions) and right panel is for cascade-like
events (mainly from $\nu_{e}$ and $\nu_{\tau}$ interactions), as defined in
section 3.3.

\begin{overpic}[width=212.47617pt]{figures/Figure4A.png}
\put(15.0,11.0){DUNE,\leavevmode\nobreak\ {\it tracks}} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure4B.png}
\put(15.0,11.0){DUNE,\leavevmode\nobreak\ {\it cascades}} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure4C.png}
\put(15.0,11.0){ORCA,\leavevmode\nobreak\ {\it tracks}} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure4D.png}
\put(15.0,11.0){ORCA,\leavevmode\nobreak\ {\it cascades}} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure4E.png}
\put(15.0,11.0){HyperK,\leavevmode\nobreak\ {\it tracks}} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure4F.png}
\put(15.0,11.0){HyperK,\leavevmode\nobreak\ {\it cascades}} \end{overpic}

Figure 4: $\Delta\chi^{2}$ sensitivity for discriminating between FeNi model
and FeNiH model in 20 years livetime of upcoming detectors. From top to
bottom, the panels show the signed $\Delta\chi^{2}$ maps (as defined in
Section 3.4) as a function of the reconstructed neutrino energy and zenith
angle, for DUNE, ORCA, and HyperKamiokande, for track-like events (left
column) and cascade-like events (right column) as defined in section 3.3. The
color scale is the same in all plots. The number indicated in each plot
corresponds to the total $\Delta\chi^{2}$ sensitivity summed (in absolute
value) over all bins of the $(E_{reco},\theta_{reco})$ plane.

Figure 5: Sensitivity of the individual upcoming detectors to the outer core
composition. Left: Sensitivity profile for the absolute precision in the Z/A
measurement, for 20 years operation of the HyperKamiokande, DUNE and ORCA
detectors. The vertical lines indicate the Z/A separation between specific
pairs of models of core composition. Right: Sensitivity bands as a function of
the detector livetime, for discriminating specific pairs of models. In the
right plot, the filled bands represent the discrimination power between
FeNiSi2O4 and FeNiH , while the hashed bands represent FeNi vs FeNiH . The
band limits in both plots correspond to the baseline and optimistic
assumptions on the detectors capability to filter out events which are
insensitive to neutrino flavor (as explained in section 3.3).

Figure 6: Sensitivity to the outer core composition for upcoming and next-
generation detectors. Left: $\Delta\chi^{2}$ profile for the absolute
precision in the Z/A measurement, for 20 years operation of a Next Generation
(NextGen) detector and for the combination of DUNE, ORCA and HyperKamiokande
(HK) over the same livetime. The vertical lines indicate the Z/A separation
between specific pairs of models of core composition. Right: Sensitivity bands
as a function of the detector livetime, for discriminating specific pairs of
models. In the right plot, the filled bands represent the discrimination power
between FeNiSi2O4 and FeNiH , while the hashed bands represent FeNi vs FeNiH
. The band limits correspond to the baseline and optimistic assumptions on
the detectors capability to filter out events which are insensitive to
neutrino flavor (as explained in section 3.3).

Figure 7:  $\Delta\chi^{2}$ sensitivity for discriminating FeNiSi2O4 vs. FeNiH
with a NextGen detector for various combinations of exposure, energy and
angular resolutions. The left plot is obtained for a fixed exposure of
$200\text{\,}\leavevmode\nobreak\ $Mton yr. The color scale indicates the
intervals between 1, 2 and 3 ${\sigma}$ ($68\%,95\%,99\%$ C.L.). Superimposed
on the maps are the lifetimes required to reach the 200 Mton yr exposure for
DUNE, ORCA and HK detectors, as well as their specific resolution in the
energy range of interest for the tomography measurement (3–7 GeV). The right
plot is obtained by assuming a linear relationship between the energy and
angular resolutions for the benchmark detectors, as suggested by the left
plot.

.

\begin{overpic}[width=212.47617pt]{figures/Figure8A.png} \end{overpic}

\begin{overpic}[width=212.47617pt]{figures/Figure8B.png} \end{overpic}

Figure 8: Sensitivity of the NextGen detector for discriminating between
FeNiSi2O4 and FeNiH outer core compositions after 20 years data taking. The
plots show the signed $\Delta\chi^{2}$ maps (as defined in equation (9)) as a
function of the reconstructed energy and zenith of the neutrinos, respectively
for track-like (left) and cascade-like (right) events. Figure 9: Precision of
the Z/A measurement achievable at 1$\sigma$ with the NextGen detector as a
function of running time. The crosses indicate the separation in Z/A between
pairs of models considered in this study. Not shown in the graph is the time
required to distinguish FeNiSi7O2 vs FeNiSi2O4 ($\simeq 90$ years) and
FeNiSi7O2 vs FeNiSiH ($\simeq 120$) years. If $n$ identical NextGen detectors
were running in parallel, the time scale would be approximately reduced by
that same factor $n$, i.e. FeNiSi7O2 vs FeNiSi2O4 could be distinguished in
about 20 years if 4 NextGen detectors would be taking data simultaneously.

Figure 10: Examples of response functions used for the modelling of the
neutrino detectors. Upper left: Effective detector mass ($M_{eff}$) as a
function of the neutrino true energy. Upper right: Classification efficiency
($\rho_{class}$) as a function of the true neutrino energy. Lower left:
Probability distribution function for the reconstructed zenith angle for a
neutrino with true energy $E_{\text{true}}=10$ GeV and true zenith angle
$\theta_{\text{true}}\approx 169^{\circ}$ (corresponding to a neutrino
trajectory grazing the inner-outer core boundary). Lower right: Probability
distribution function for the reconstructed energy for two specific values of
the true neutrino energy ($E_{\text{true}}=1\;\text{and}\;5$ GeV).


66institutetext: Center for Biomedical Image Computing and Analytics (CBICA),
University of Pennsylvania, Philadelphia, PA, USA 22institutetext: Department
of Pathology and Laboratory Medicine, Perelman School of Medicine, University
of Pennsylvania, Philadelphia, PA, USA 33institutetext: Department of
Radiology, Perelman School of Medicine, University of Pennsylvania,
Philadelphia, PA, USA 44institutetext: Department of Informatics, Technical
University of Munich, Munich, Germany 55institutetext: Department of
Biomedical Informatics, Stony Brook University, Stony Brook, NY, USA
66institutetext: Department of Computer Science, Stony Brook University, Stony
Brook, NY, USA
* Corresponding author: 66email<EMAIL_ADDRESS>
# Federated Learning for the Classification of Tumor Infiltrating Lymphocytes

Ujjwal Baid 112233 Sarthak Pati 11223344 Tahsin M. Kurc 55 Rajarsi Gupta 55
Erich Bremer 55 Shahira Abousamra 66 Siddhesh P. Thakur 112233 Joel H. Saltz
55 Spyridon Bakas 112233**

###### Abstract

We evaluate the performance of federated learning (FL) in developing deep
learning models for analysis of digitized tissue sections. A classification
application was considered as the example use case, on quantifiying the
distribution of tumor infiltrating lymphocytes within whole slide images
(WSIs). A deep learning classification model was trained using $50\times 50$
square micron patches extracted from the WSIs. We simulated a FL environment
in which a dataset, generated from WSIs of cancer from numerous anatomical
sites available by The Cancer Genome Atlas repository, is partitioned in 8
different nodes. Our results show that the model trained with the federated
training approach achieves similar performance, both quantitatively and
qualitatively, to that of a model trained with all the training data pooled at
a centralized location. Our study shows that FL has tremendous potential for
enabling development of more robust and accurate models for histopathology
image analysis without having to collect large and diverse training data at a
single location.

###### Keywords:

federated learning, classification, histopathology, digital pathology, tumor
infiltrating lymphocytes

##  1 Introduction

Advances in machine learning (ML), and particularly deep learning (DL), have
shown promise in healthcare. However, the availability of large amounts of
data with increased diversity is essential to produce accurate and
generalizable models [1]. This is currently addressed by pooling data to a
centralized location, typically facilitated via use-inspired consortia [2, 3,
4, 5, 6, 7]. However, this centralization is challenging (and at times
infeasible) due to numerous concerns relating to privacy, data-ownership,
intellectual property, as well as compliance with varying regulatory policies
(e.g., Health Insurance Portability and Accountability Act (HIPAA) of the
United States [8] and the General Data Protection Regulation (GDPR) of the
European Union [9]).

In contrast to this centralized paradigm, “federated learning” (FL) describes
an approach where ML/DL models are getting trained only by sharing model
updates, while all data are always retained locally within the acquiring
institution [10, 11, 12, 13, 14]. Several studies have demonstrated that the
performance of FL models is comparable to their equivalent CL-trained models
[15, 16, 11, 12, 17, 18, 19]. As ML/DL methods increasingly become the primary
means of analyzing large datasets, FL can offer a potential paradigm shift for
multi-institutional collaborations, alleviating the need for data sharing, and
hence increase access to geographically-distinct collaborators, thereby
increasing the size and diversity of data used to train ML/DL models.

Electronic capture (digitisation) and analyses of whole slide images (WSIs) of
tissue specimens are becoming ubiquitous. Digital Pathology interpretation is
becoming increasingly common, where many sites are actively scanning archived
glass tissue slides with commercially available high speed scanners to
generate high-resolution gigapixel WSIs. Alongside these efforts, a great
variety of AI algorithms have been developed to extract many salient tissue
and tumor characteristics from WSIs. Examples include segmentation of tumor
regions, histologic subtypes of tumors, microanatomic tissue compartments;
detection and classification of immune cells to identify tumor infiltrating
lymphocytes (TILs); and the detection and classification of cells and nuclei.
TILs are lymphoplasmacytic cells that are spatially located in tumor regions,
where their role as an important biomarker for the prediction of clinical
outcomes in cancer patients is becoming increasingly recognised [20, 21, 22].
Identification of the abundance and the patterns of spatial distribution of
TILs in WSI represent a quantitative approach to characterizing important
tumor immune interactions.

The Stony Brook Biomedical Informatics group has actively contributed to this
field of work for many years by characterizing the performance of DL pathology
algorithms. Specifically, we have previously developed a DL based method to
analyze WSIs to quantify distributions of TILs [23, 24]. This method
partitions a WSI into a regular mesh of image patches. Each image patch covers
an area of $50\times 50$ square microns, which is equivalent to $100\times
100$ and $200\times 200$ square pixels in a WSI captured at $20\times$ and
$40\times$ magnification level, respectively. The method trains a
classification model (based on a VGG16 pre-trained on ImageNet) to predict if
a given image patch is TIL-positive (i.e., the patch contains 2 or more
lymphocytes) or TIL-negative. This classification model has been trained and
evaluated with a set of TIL-positive and TIL-negative patches from multiple
cancer types with comprehensive analyses carried out to characterize
performance of this method [24].

In this present study, we study the performance of federated learning using DL
based detection of TILs as our model application. It is particularly useful to
leverage an algorithm and WSI datasets with well known performance
characteristics so that performance differences arising from FL can be easily
understood. In our study, we assume that each node is responsible for all
training required for slides from one or more cancer sites. This ensures that
we will see significant out-of-distribution impacts as training carried out on
WSIs from a given tumor type often imperfectly generalizes to WSIs from
different cancer sites.

##  2 Methods

###  2.1 Data

We created a training and validation dataset consisting of patches extracted
from WSIs of cancer from 12 anatomical sites, comprising breast (BRCA), cervix
(CESC), colon (COAD), lung (LUAD and LUSC), pancreas (PAAD), prostate (PRAD),
rectum (READ), skin (SKCM), stomach (STAD), uterus (UCEC), uvea of the eye
(UVM) cases, publicly available in The Cancer Genome Atlas (TCGA) [25]. Tables
1 and 2 show the explicit breakdown of the dataset into individual network
sites (Site 1-8).

Table 1: Data sharding to collaborative network sites according to anatomy. Node | Anatomy
---|---
Site1 | Cervical squamous cell carcinoma (CESC)
Site2 | Lung squamous cell carcinoma (LUSC)
Site3 | Colon adenocarcinoma (COAD)
Pancreatic adenocarcinoma (PAAD)
Uterine corpus endometrial carcinoma (UCEC)
Site4 | Rectum adenocarcinoma (READ)
Site5 | Stomach adenocarcinoma (STAD)
Site6 | Uveal melanoma (UVM)
Site7 | Lung adenocarcinoma (LUAD)
Site8 | Breast invasive carcinoma (BRCA)
Prostate adenocarcinoma (PRAD)
Skin cutaneous melanoma (SKCM)

###  2.2 Data Sharding

The complete dataset was split into vertical partitions, i.e., include all the
patients of a given type of cancer and assign all the associated training data
to a single site (as is shown in Table 2). With this strategy in mind, eight
virtual network sites were created with each virtual site assigned to a
separate computational node.

Table 2: Data split amongst different network sites for training and validation cohorts Cohort |  | Total
---
Patients
| Total
---
Patches
| Patients
---
in training
| Training
---
patches
| Patients
---
in validation
| Validation
---
patches
Site1 | 111 | 13144 | 89 | 10542 | 22 | 2602
Site2 | 217 | 25624 | 174 | 20517 | 43 | 5107
Site3 | 4 | 11276 | 3 | 11039 | 1 | 237
Site4 | 40 | 4743 | 32 | 3793 | 8 | 950
Site5 | 196 | 23258 | 156 | 18505 | 40 | 4753
Site6 | 24 | 2438 | 19 | 1938 | 5 | 500
Site7 | 10 | 23336 | 8 | 16265 | 2 | 7091
Site8 | 60 | 51522 | 48 | 39665 | 12 | 11857

###  2.3 Model Architecture

In this study we used the VGG network architecture [26, 27] (Figure 1), which
is well-known for performing classification and regression workloads, on the
ImageNet classification challenge [28].

Figure 1: Deep convolutional neural network based VGG type architecture. This
figure was created using PlotNeuralNetwork tool [29].

The VGG network used in the centralized configuration evaluated previously has
16 convolutional layers and 3 dense layers. We have modified the final
classifier layers to include a global average pooling layer followed by a
single dense layer, which allows greater flexibility for different types of
workloads [30]. VGG uses $3\times 3$ convolution filters and $2\times 2$ max-
pooling layers with a stride of 2 throughout the architecture. The original
architecture uses ReLU activation function [31] and categorical cross-entropy
loss function. The initial layers of the VGG perform feature extraction and
the last sigmoid layers act as the classifier.

The model configuration used in the federated learning setting, was
initialized using randomized weights and used as the starting point for all
sites. The VGG network is trained with patch size of $300\times 300$, are
normalised between $[0,1]$. Sigmoid layer is used as final layer for
classification. The model is trained with initial learning rate of $0.001$.
The cross entropy loss with adam optimizer is used to update the weights of
the model during training.

###  2.4 Model Training and Prediction

In a centralized training scenario, all of the training and validation patches
are used to train a classification model. In the inference (prediction) phase,
an input WSI is partitioned into patches of $50\times 50$ square microns. Each
patch is processed by the trained model, and the model assigns a probability
value between $[0.0,1.0]$ of the patch being TIL-positive. After all of the
patches are processed, a probability map is created to show the distribution
of TIL patches across the WSI.

Each “federated training round” is defined as each site training the model on
their data for a pre-defined number of epochs and providing model updates to
be aggregated. For this study, the model was trained for $500$ federated
rounds, where each site trained the model for a single epoch for each
federated round by leveraging an independent aggregation server [13].

###  2.5 Experimental Design

The complete study described here was facilitated by the model and training
infrastructure provided by the Generally Nuanced Deep Learning Framework
(GaNDLF) [32]. We federate [11, 12, 13] the model and the entire training
pipeline using the Open Federated Learning (OpenFL) library [33], which allows
the model to be trained across multiple sites across the collaborative network
without sharing any data. The quantitative performance evaluation was done
using balanced classification accuracy [34].

##  3 Results

###  3.1 Quantitative Evaluation

The performance evaluation in terms of balanced classification accuracy for
each individual site and the consensus model is shown in Table 3. In the
centralized training scenario we achieved accuracy of $0.75$.

Table 3: Performance evaluation in terms of balanced classification accuracy for individual sites and the consensus model on the validation dataset Model(row)\Data(col) | Site1 | Site2 | Site3 | Site4 | Site5 | Site6 | Site7 | Site8 | Average
---|---|---|---|---|---|---|---|---|---
Site 1 model | 0.90 | 0.88 | 0.78 | 0.80 | 0.86 | 0.97 | 0.86 | 0.82 | 0.86
Site 2 model | 0.90 | 0.90 | 0.72 | 0.84 | 0.88 | 0.97 | 0.84 | 0.79 | 0.85
Site 3 model | 0.85 | 0.84 | 0.87 | 0.82 | 0.82 | 0.98 | 0.96 | 0.84 | 0.87
Site 4 model | 0.78 | 0.80 | 0.56 | 0.86 | 0.81 | 0.93 | 0.63 | 0.72 | 0.76
Site 5 model | 0.88 | 0.88 | 0.70 | 0.86 | 0.88 | 0.97 | 0.80 | 0.79 | 0.85
Site 6 model | 0.78 | 0.74 | 0.59 | 0.76 | 0.72 | 1.00 | 0.97 | 0.73 | 0.79
Site 7 model | 0.82 | 0.79 | 0.75 | 0.76 | 0.75 | 0.97 | 0.97 | 0.80 | 0.83
Site 8 model | 0.83 | 0.81 | 0.71 | 0.81 | 0.79 | 0.99 | 0.94 | 0.88 | 0.84
Consensus model | 0.88 | 0.88 | 0.81 | 0.85 | 0.86 | 0.98 | 0.95 | 0.87 | 0.89

###  3.2 Qualitative Evaluation

Figure 2(a) shows a visualization of patch-level probabilities from the
original TIL classification model trained with all of the data in a
centralized location. The patch-level probabilities are stitched together to
generate a TIL heatmap overlaid on the source WSI. The RED-colored regions
indicate high probability of TILs predicted by the models. The TIL predictions
from the federated and site-specific models are shown in Figure 2(b) and (c),
respectively. The federated consensus model predicts a TIL map that is
qualitatively similar to the curated TIL map generated by the centralized
model. In this example, the site-specific model predictions also appear quite
similar even though this model is not trained with the entire training
dataset. Even though these TIL maps look similar overall, qualitative
differences are evident in terms of both the abundance and spatial
distribution of TILs when comparing the federated consensus model with the
site-specific model.

Figure 2: Predictions from different models displayed as heatmaps overlaid on
the source whole slide image. (a) Predictions from the original model trained
with all of the training data at a centralized location. (b) Predictions from
the consensus model trained with the federated learning approach using data
from all the sites. (c) Predictions from the model trained with data from a
single site only.

##  4 Discussion

In this study, we generated anatomic site-specific and a federated consensus
model to classify tissue regions that contain TILs. Overall, TIL maps
generated with the consensus model trained in a federated manner appear
comparable to the output of the original VGG16 TIL model, trained by centrally
collecting all of the data in a single location [23]. We have also shown that
there are differences in performance when the TIL models are trained with
anatomic site-specific data. We evaluated performance by using held-out local
validation data from each anatomic cancer site. We also performed qualitative
analysis on completely unseen held-out data.

The detection of TILs across WSIs of cancer was conceived based on the fairly
consistent appearance of lymphocytes in normal and cancer tissues across
anatomical sites. Lymphocytes are typically 8-12 um in size, round to ovoid in
shape with dark blue-purple nuclei, and contain scant cytoplasm. However,
detecting lymphocytes, which are called TILs within the spatial context of
cancer, is challenging based on the ability of the models to distinguish TILs
from cancer and other types of cells that may have overlapping features.
Therefore, the original TIL model was trained with a centralized approach to
provide the model with plenty of examples of non-lymphocytes and non-TILs
across the complex microscopic landscapes of both normal and cancer tissues to
develop a robust pan-cancer model.

Figure 3 expands upon Figure 2 to provide further insights into the nuances of
model performance, where site-specific models (from Sites 1, 2, 3, and 5)
generate TIL maps that appear quite similar. The TIL maps generated with
models using training data from sites 4 and 8 also appear similar to one
another with increased TILs in the left portion of the WSI at the left side of
the displayed figure, which is qualitatively different from the TIL maps from
site-specific models 1, 2, 3, and 5. The TIL map generated with the model from
site 7 appears more conservative in comparison to models from sites 1-5 and 8,
in terms of predicting less TILs globally. Notably, the model from site 6 does
not predict any TILs since uveal melanoma (UVM) of the eye is typically not
associated with TILs, which indicates the scarcity of TIL positive training
patches. Furthermore, UVM is actually used as a negative control to test the
specificity of the original TIL model in distinguishing TILs from melanoma
tumor cells that can have features which mimic lymphocytes.

Figure 3: TIL heatmaps from federated training corresponding to the
hematoxylin and eosin (H&E) stained WSI of breast cancer, TCGA-AO-A0JC-DX1.
The TIL map from the original VGG16 TIL model, trained with all of the
training data at a centralized location, is shown before displaying predicted
TIL maps from models trained at each anatomic site. The consensus model
trained with the federated training approach using data from all the sites is
shown afterwards. Visually representing the TIL maps from the various models
allows qualitative evaluation of model performance for each site and consensus
with respect to the original centralized VGG16 TIL model within the federated
training paradigm.

The performance of the models trained on site-specific data is dependent on
the distribution of training data, specifically whether all training labels
were present and if there were enough subjects in the site to warrant
meaningful training. For example, it is observed (Table 3) that Site 6 has a
perfect classification performance on its own validation dataset. After
further investigation we confirmed that the validation data of Site 6 includes
only non-TIL patches. This shows the value of training such models in a
federated setting, which allows sites that do not have enough data to be able
to learn from a global consensus, thus improving their own model performance.
Also, from Table 3 it can be concluded that the federated consensus model
performed better than the individual sites in terms of average balanced
classification accuracy.

We have evaluated the federated training paradigm using a single network
architecture, and this can be expanded to include multiple architecture types
to analyze which one is more suitable with regards to communication
efficiency, overall convergence stability, and data privacy. Additionally, a
more rigorous analysis with varying data splits is needed to provide a
holistic picture of the federated training paradigm in this vertical
situation. Finally, the entire study was based on simulating a federated
learning paradigm on a single institution’s network, but still using different
computational nodes, and needs to be expanded to multiple institutions to
showcase a more realistic deployment.

## Funding

Research reported in this publication was partly supported by the National
Institutes of Health (NIH) under award numbers NIH/NCI:U01CA242871,
NIH/NINDS:R01NS042645, NIH/NCI:U24CA189523, NIH/NCI:UG3CA225021,
NIH/NCI:4UH3CA225021, NIH/NCI:U24CA215109, and NCI/NIH:U24CA180924. The
content of this publication is solely the responsibility of the authors and
does not represent the official views of the NIH.

## References

  * [1] Z. Obermeyer and E. J. Emanuel, “Predicting the future—big data, machine learning, and clinical medicine,” The New England journal of medicine, vol. 375, no. 13, p. 1216, 2016.
  * [2] S. G. Armato III, G. McLennan, M. F. McNitt-Gray, C. R. Meyer, D. Yankelevitz, D. R. Aberle, C. I. Henschke, E. A. Hoffman, E. A. Kazerooni, H. MacMahon, et al., “Lung image database consortium: developing a resource for the medical imaging research community,” Radiology, vol. 232, no. 3, pp. 739–748, 2004.
  * [3] P. M. Thompson, J. L. Stein, S. E. Medland, D. P. Hibar, A. A. Vasquez, M. E. Renteria, R. Toro, N. Jahanshad, G. Schumann, B. Franke, et al., “The enigma consortium: large-scale collaborative analyses of neuroimaging and genetic data,” Brain imaging and behavior, vol. 8, no. 2, pp. 153–182, 2014.
  * [4] T. G. Consortium, “Glioma through the looking GLASS: molecular evolution of diffuse gliomas and the Glioma Longitudinal Analysis Consortium,” Neuro-Oncology, vol. 20, pp. 873–884, 02 2018.
  * [5] S. Bakas, D. R. Ormond, K. D. Alfaro-Munoz, M. Smits, L. A. D. Cooper, R. Verhaak, and L. M. Poisson, “iglass: imaging integration into the glioma longitudinal analysis consortium,” Neuro-oncology, vol. 22, no. 10, pp. 1545–1546, 2020.
  * [6] C. Davatzikos, J. S. Barnholtz-Sloan, S. Bakas, R. Colen, A. Mahajan, C. B. Quintero, J. Capellades Font, J. Puig, R. Jain, A. E. Sloan, et al., “Ai-based prognostic imaging biomarkers for precision neuro-oncology: the respond consortium,” Neuro-oncology, vol. 22, no. 6, pp. 886–888, 2020\.
  * [7] M. Habes, R. Pomponio, H. Shou, J. Doshi, E. Mamourian, G. Erus, I. Nasrallah, L. J. Launer, T. Rashid, M. Bilgel, et al., “The brain chart of aging: Machine-learning analytics reveals links between brain aging, white matter disease, amyloid burden, and cognition in the istaging consortium of 10,216 harmonized mr scans,” Alzheimer’s & Dementia, vol. 17, no. 1, pp. 89–102, 2021.
  * [8] G. J. Annas et al., “Hipaa regulations-a new era of medical-record privacy?,” New England Journal of Medicine, vol. 348, no. 15, pp. 1486–1490, 2003.
  * [9] P. Voigt and A. Von dem Bussche, “The eu general data protection regulation (gdpr),” A Practical Guide, 1st Ed., Cham: Springer International Publishing, vol. 10, p. 3152676, 2017.
  * [10] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in Artificial intelligence and statistics, pp. 1273–1282, PMLR, 2017.
  * [11] M. J. Sheller, G. A. Reina, B. Edwards, J. Martin, and S. Bakas, “Multi-institutional deep learning modeling without sharing patient data: A feasibility study on brain tumor segmentation,” in International MICCAI Brainlesion Workshop, pp. 92–104, Springer, 2018.
  * [12] M. J. Sheller, B. Edwards, G. A. Reina, J. Martin, S. Pati, A. Kotrotsou, M. Milchenko, W. Xu, D. Marcus, R. R. Colen, et al., “Federated learning in medicine: facilitating multi-institutional collaborations without sharing patient data,” Scientific reports, vol. 10, no. 1, pp. 1–12, 2020\.
  * [13] N. Rieke, J. Hancox, W. Li, F. Milletari, H. R. Roth, S. Albarqouni, S. Bakas, M. N. Galtier, B. A. Landman, K. Maier-Hein, et al., “The future of digital health with federated learning,” NPJ digital medicine, vol. 3, no. 1, pp. 1–7, 2020.
  * [14] I. Dayan, H. R. Roth, A. Zhong, A. Harouni, A. Gentili, A. Z. Abidin, A. Liu, A. B. Costa, B. J. Wood, C.-S. Tsai, et al., “Federated learning for predicting clinical outcomes in patients with covid-19,” Nature medicine, vol. 27, no. 10, pp. 1735–1743, 2021.
  * [15] K. Chang, N. Balachandar, C. Lam, D. Yi, J. Brown, A. Beers, B. Rosen, D. L. Rubin, and J. Kalpathy-Cramer, “Distributed deep learning networks among institutions for medical imaging,” Journal of the American Medical Informatics Association, vol. 25, no. 8, pp. 945–954, 2018.
  * [16] A. Nilsson, S. Smith, G. Ulm, E. Gustavsson, and M. Jirstrand, “A performance evaluation of federated learning algorithms,” in Proceedings of the Second Workshop on Distributed Infrastructures for Deep Learning, pp. 1–8, 2018\.
  * [17] K. V. Sarma, S. Harmon, T. Sanford, H. R. Roth, Z. Xu, J. Tetreault, D. Xu, M. G. Flores, A. G. Raman, R. Kulkarni, et al., “Federated learning improves site performance in multicenter deep learning without data sharing,” Journal of the American Medical Informatics Association, vol. 28, no. 6, pp. 1259–1264, 2021.
  * [18] C. Shen, P. Wang, H. R. Roth, D. Yang, D. Xu, M. Oda, W. Wang, C.-S. Fuh, P.-T. Chen, K.-L. Liu, et al., “Multi-task federated learning for heterogeneous pancreas segmentation,” in Clinical Image-Based Procedures, Distributed and Collaborative Learning, Artificial Intelligence for Combating COVID-19 and Secure and Privacy-Preserving Machine Learning, pp. 101–110, Springer, 2021.
  * [19] D. Yang, Z. Xu, W. Li, A. Myronenko, H. R. Roth, S. Harmon, S. Xu, B. Turkbey, E. Turkbey, X. Wang, et al., “Federated semi-supervised learning for covid region segmentation in chest ct using multi-national data from china, italy, japan,” Medical image analysis, vol. 70, p. 101992, 2021.
  * [20] B. Mlecnik, M. Tosolini, A. Kirilovsky, A. Berger, G. Bindea, T. Meatchi, P. Bruneval, Z. Trajanoski, W.-H. Fridman, F. Pagès, et al., “Histopathologic-based prognostic factors of colorectal cancers are associated with the state of the local immune reaction,” Journal of clinical oncology, vol. 29, no. 6, pp. 610–618, 2011.
  * [21] G. Badalamenti, D. Fanale, L. Incorvaia, N. Barraco, A. Listi, R. Maragliano, B. Vincenzi, V. Calo, J. L. Iovanna, V. Bazan, et al., “Role of tumor-infiltrating lymphocytes in patients with solid tumors: Can a drop dig a stone?,” Cellular immunology, vol. 343, p. 103753, 2019.
  * [22] G. E. Idos, J. Kwok, N. Bonthala, L. Kysh, S. B. Gruber, and C. Qu, “The prognostic implications of tumor infiltrating lymphocytes in colorectal cancer: a systematic review and meta-analysis,” Scientific reports, vol. 10, no. 1, pp. 1–14, 2020.
  * [23] J. Saltz, R. Gupta, L. Hou, T. Kurc, P. Singh, V. Nguyen, D. Samaras, K. R. Shroyer, T. Zhao, R. Batiste, et al., “Spatial organization and molecular correlation of tumor-infiltrating lymphocytes using deep learning on pathology images,” Cell reports, vol. 23, no. 1, pp. 181–193, 2018\.
  * [24] S. Abousamra, M. Gupta, L. Hou, R. Batiste, T. Zhao, A. Shankar, A. Rao, C. Chen, D. Samaras, T. Kurc, et al., “Deep learning-based mapping of tumor infiltrating lymphocytes in whole slide images of 23 types of cancer,” Frontiers in oncology, p. 5971, 2021.
  * [25] S. Abousamra, L. Hou, R. Gupta, C. Chen, D. Samaras, T. Kurc, R. Batiste, T. Zhao, S. Kenneth, and J. Saltz, “Learning from thresholds: fully automated classification of tumor infiltrating lymphocytes for multiple cancer types,” arXiv preprint arXiv:1907.03960, 2019.
  * [26] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” arXiv preprint arXiv:1409.1556, 2014.
  * [27] A. Ben-Cohen, I. Diamant, E. Klang, M. Amitai, and H. Greenspan, “Fully convolutional network for liver segmentation and lesions detection,” in Deep learning and data labeling for medical applications, pp. 77–85, Springer, 2016.
  * [28] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in 2009 IEEE conference on computer vision and pattern recognition, pp. 248–255, Ieee, 2009.
  * [29] H. Iqbal, “Harisiqbal88/plotneuralnet v1.0.0,” Dec. 2018.
  * [30] T.-Y. Hsiao, Y.-C. Chang, H.-H. Chou, and C.-T. Chiu, “Filter-based deep-compression with global average pooling for convolutional networks,” Journal of Systems Architecture, vol. 95, pp. 9–18, 2019.
  * [31] A. F. Agarap, “Deep learning using rectified linear units (relu),” arXiv preprint arXiv:1803.08375, 2018.
  * [32] S. Pati, S. P. Thakur, M. Bhalerao, U. Baid, C. Grenko, B. Edwards, M. Sheller, J. Agraz, B. Baheti, V. Bashyam, et al., “Gandlf: A generally nuanced deep learning framework for scalable end-to-end clinical workflows in medical imaging,” arXiv preprint arXiv:2103.01006, 2021.
  * [33] G. A. Reina, A. Gruzdev, P. Foley, O. Perepelkina, M. Sharma, I. Davidyuk, I. Trushkin, M. Radionov, A. Mokrov, D. Agapov, et al., “Openfl: An open-source framework for federated learning,” arXiv preprint arXiv:2105.06413, 2021.
  * [34] K. H. Brodersen, C. S. Ong, K. E. Stephan, and J. M. Buhmann, “The balanced accuracy and its posterior distribution,” in 2010 20th international conference on pattern recognition, pp. 3121–3124, IEEE, 2010.