# Theoretical Antineutrino Detection, Direction and Ranging at Long Distances

Glenn R. Jocher<EMAIL_ADDRESS>Integrity Applications Incorporated,
15020 Conference Center Drive, Chantilly, VA, 20151 USA  Daniel A. Bondy
<EMAIL_ADDRESS>Integrity Applications Incorporated, 15020
Conference Center Drive, Chantilly, VA, 20151 USA  Brian M. Dobbs
<EMAIL_ADDRESS>Integrity Applications Incorporated, 15020
Conference Center Drive, Chantilly, VA, 20151 USA  Stephen T. Dye
<EMAIL_ADDRESS>College of Natural Sciences
Hawaii Pacific University, Kaneohe, HI 96744 USA  Department of Physics and
Astronomy
University of Hawaii, Honolulu, HI, 96822 USA  James A. Georges III
<EMAIL_ADDRESS>Integrity Applications Incorporated, 15020
Conference Center Drive, Chantilly, VA, 20151 USA  John G. Learned
<EMAIL_ADDRESS>Department of Physics and Astronomy
University of Hawaii, Honolulu, HI, 96822 USA  Christopher L. Mulliss
<EMAIL_ADDRESS>Integrity Applications Incorporated, 15020
Conference Center Drive, Chantilly, VA, 20151 USA  Shawn Usman
<EMAIL_ADDRESS>InnoVision Basic and Applied Research Office, Sensor
Geopositioning Center
National Geospatial-Intelligence Agency, 7500 GEOINT Dr., Springfield, VA,
22150 USA

###### Abstract

In this paper we introduce the concept of what we call “NUDAR” (NeUtrino
Direction and Ranging), making the point that measurements of the observed
energy and direction vectors can be employed to passively deduce the exact
three-dimensional location and thermal power of geophysical and anthropogenic
neutrino sources from even a single detector. Earlier studies have presented
the challenges of long-range detection, dominated by the unavoidable inverse-
square falloff in neutrinos, which force the use of kiloton scale detectors
beyond a few kilometers. Earlier work has also presented the case for multiple
detectors, and has reviewed the background challenges. We present the most
precise background estimates to date, all handled in full three dimensions, as
functions of depth and geographical location. For the present calculations, we
consider a hypothetical 138 kiloton detector which can be transported to an
ocean site and deployed to an operational depth. We present a Bayesian
estimation framework to incorporate any a priori knowledge of the reactor that
we are trying to detect, as well as the estimated uncertainty in the
background and the oscillation parameters. Most importantly, we fully employ
the knowledge of the reactor spectrum and the distance-dependent effects of
neutrino oscillations on such spectra. The latter, in particular, makes
possible determination of range from one location, given adequate signal
statistics. Further, we explore the rich potential of improving detection with
even modest improvements in individual neutrino direction determination. We
conclude that a 300 MWth reactor can indeed be geo-located, and its operating
power estimated with one or two detectors in the hundred kiloton class at
ranges out to a few hundred kilometers. We note that such detectors would have
natural and non-interfering utility for scientific studies of geo-neutrinos,
neutrino oscillations, and astrophysical neutrinos. This motivates the
development of cost effective methods of constructing and deploying such next
generation detectors.

###### keywords:

antineutrino , neutrino , geo-neutrino , reactor , geo-reactor , oscillation

††journal: Physics Reports

###### Contents

  1. 1 Introduction
    1. 1.1 Types of neutrinos
    2. 1.2 Detection mechanism: inverse beta decay
    3. 1.3 Electron antineutrino production in reactors
    4. 1.4 Detectors
    5. 1.5 Background Sources
    6. 1.6 Other physics and astrophysics
  2. 2 Antineutrino Detector Model
    1. 2.1 Inverse beta decay model
    2. 2.2 Photon propagation model
    3. 2.3 Detector hardware model
    4. 2.4 Validation results
      1. 2.4.1 KamLAND energy resolution validation
      2. 2.4.2 CHOOZ direction resolution validation
      3. 2.4.3 TREND energy and direction resolution MC results
  3. 3 Geospatial Model
    1. 3.1 Reactor Neutrinos
    2. 3.2 Geo-neutrinos
    3. 3.3 Non-neutrino background
    4. 3.4 Uncertainty model
  4. 4 Estimation theory
    1. 4.1 Optimal MAP estimator
    2. 4.2 Applied suboptimal MAP estimator
      1. 4.2.1 Reactor geolocation
      2. 4.2.2 Antineutrino oscillation parameter estimation
  5. 5 Oscillation parameter estimation results
    1. 5.1 Estimator Cramer-Rao lower bound
    2. 5.2 Cramer-Rao lower bound verification results
    3. 5.3 Optimal detector placement via CRLB results
    4. 5.4 Oscillation parameter search MC results
  6. 6 Reactor geolocation results
    1. 6.1 No unknown reactor (Reactor Exclusion)
    2. 6.2 High reactor background (Europe-Atlantic)
      1. 6.2.1 One detector
      2. 6.2.2 Two detectors
    3. 6.3 Low reactor background with nearby reactor (Southern Africa)
      1. 6.3.1 One detector, detector #1
      2. 6.3.2 One detector, detector #2
      3. 6.3.3 Two detectors
    4. 6.4 High reactor background with nearby reactors (Europe-Mediterranean)
      1. 6.4.1 Three detectors
  7. 7 Conclusions and Recommendations
  8. 8 Acronyms
  9. 9 Acknowledgements
  10. A Geo-reactor results

##  1 Introduction

In this paper we present a careful discussion of scientific and applied
experiments based on the measurement of electron antineutrinos from distant
sources (tens to hundreds of km) in a $10^{34}$ proton class detector. Because
of the well-known ghostly nature of the neutrinos, we focus on large detectors
and strive to employ all information potentially available, including our best
present knowledge of all background processes. One cannot escape the small
neutrino cross-sections or the inverse-square falloff in neutrino flux with
range, which requires the detector size to grow with the square of the range.
One also cannot escape the increasing importance of the background at long
distances as the source signal becomes a relatively small component of the
total signal. Since some (dominant) background sources depend on cosmic rays
which come into the atmosphere, one is driven to place the larger detectors
deeper underground (or under water) at longer distances.

###  1.1 Types of neutrinos

Since the original observation of neutrinos by Cowan and Reines in the mid
1950’s [1], operating in close proximity to nuclear reactors, people have
speculated on the possible detection of reactors from longer ranges. Note that
we focus here on the detection of electron antineutrinos in the energy range
of about 1 to 11 MeV. The flux of neutrinos from reactors is relatively small
compared to antineutrinos ($<1\%$). Additionally, flux from the sun is
composed almost entirely of neutrinos and not antineutrinos. Experiments have
set strong limits on the flux of antineutrinos from the sun [4]. The neutrino
flux from radioactive decays throughout the Earth is composed almost entirely
of electron antineutrinos [93]. Hence we can speak of the electron
antineutrinos from reactors and the Earth (geo-neutrinos), without ambiguity,
as neutrinos.

A useful complication (as we shall see) is that the flux of electron
antineutrinos consists of varying fractions of all three types of
antineutrinos (electron, muon and tauon) when observed from distances beyond a
few kilometers. For the energies and the detection mechanism under
consideration here, however, only the electron type is detectable due to the
masses of the muon and tauon. All three types of neutrinos may be sensed with
a coherent neutrino scattering detector [94], but such devices will only be
useful at small distances ($<$1km) from a reactor and will be dominated by
solar neutrinos at larger distances.

###  1.2 Detection mechanism: inverse beta decay

Remarkably, the detection mechanism for reactor neutrinos remains exactly the
same as in the original Reines-Cowan discovery studies: inverse beta decay
(IBD), detected in a scintillating medium with surrounding optical sensors. In
the charge exchange process, the electron antineutrino becomes a positron and
the target free proton becomes a neutron. The positron receives energy
proportional to the incoming neutrino energy, less the threshold of 1.8 MeV
for the reaction, and the positron scatters most often nearly perpendicular to
the incoming direction. The neutron, being much heavier than the positron,
acquires little velocity and rolls forward with tens of keV of kinetic energy
while the positron may have several MeV. One may think of the positron as
keeping the kinetic energy of the neutrino and the neutron as inheriting the
momentum.

The neutron slows to thermal velocities and gets captured within microseconds,
as compared to the positron stopping and annihilating with an electron
(usually after forming positronium) within nanoseconds. The net signature of
the IBD reaction is two spatially correlated depositions of energy, resulting
from positron energy loss and annihilation and from nuclear de-excitation
after neutron capture. These are detected by optical sensors as flashes of
light in a scintillating liquid medium, two pulses occurring close in space
(less than 1 m) and time (few microseconds) with the second pulse having a
known intensity. This distinctive double hit signature distinguishes this
reaction from the many processes which produce only one flash (such as from
solar neutrinos and many types of radioactivity). Large detectors such as the
1000 ton Kamioka Liquid-scintillator Anti-Neutrino Detector (KamLAND) [74] and
Borexino [56] have demonstrated clean extraction of a few IBD events per month
from a single pulse background rate more than a million times greater. We
discuss this process in more detail later in the discussion of neutrino fluxes
and cross sections in Section 4.

###  1.3 Electron antineutrino production in reactors

The source of the neutrinos in the nuclear fire of the reactor originates in
beta decay, which is simply neutron decay to a less massive proton, an
electron (the beta particle) and the generally unseen electron antineutrino.
Nuclear reactors operate by fission chain reactions, wherein the fission of a
heavy nucleus such as U-235 becomes two lighter but unstable elements, freeing
typically 7 neutrons, and liberating approximately 200 MeV in total energy
[15]. Most of the resulting nuclei are well away from the nuclear stable
valley and are neutron rich, resulting in many beta decays. The calculation of
the flux of neutrinos involves, as one may imagine, substantial nuclear
physics and we use the results of existing calculations herein. It should not
be a surprise that the spectrum may be approximated by an exponentially
falling flux with neutrino energy.

In reality though, both the neutrino flux and energy spectrum from a reactor
are not constant in time. They depend upon the design of the reactor, its
operation, and the reactor fuel mix. (The measured flux at the detector is
calculable to within a few percent given the source information.) There is
also a well -known predicted difference in the neutrino spectra between
uranium and plutonium, with roughly a 5% slope difference above 4 MeV. This
“burn-up” effect has been observed by the SONGS1 collaboration, operating
cubic meter scale detectors about 20 m from a power reactor [25]. This effect
has been the source of much interest as to the possibility of directly
detecting the fuel mix and observing the enrichment of Pu over time as the U
is consumed (see Nieto et al. [90]). However, since the fuel is always a
mixture of U and Pu, the statistical requirements for spectral determination
are on the order of tens of thousands of events, and something not as yet
achieved by detectors more than 20m from a fission source [25]. Thus, for
longer ranges, it is reasonable to assume a reactor energy spectrum which is
constant in time.

Thermal power output from a reactor (typically three times the electrical
power) tracks the neutrino flux to within about one percent. Since this is
also about the precision to which power output can be measured at reactors,
and of the order of the accuracy of the flux calculations, we can calculate
the expected neutrino flux from a reactor facility, which usually reports
power output (typically daily or more often) to about 3% including all sources
of error. See Lasserre, et al. [2] for more discussion of this topic.
Generally, the reactor can be treated as a point source of neutrino emission
for ranges beyond tens of meters since most of the emission occurs within
about a one meter radius of the center.

The inverse beta decay cross section for $\bar{\nu}_{e}$ upon free protons is
basic to weak interactions, and is well studied. It is proportional to the
square of the visible neutrino energy in this energy regime, and on the scale
of $10^{-42}\mathrm{cm}^{2}$ (see Equation 49). The target and detection
medium are usually the same for this purpose, typically a long organic
molecule with roughly twice as many hydrogen nuclei as carbon nuclei.
Interestingly, for many materials the density of hydrogen in chemical form is
not far from that of liquid hydrogen. Hence one cannot do much to increase the
free proton density. The protons bound in heavier nuclei do not help, since it
requires many MeV of energy to get them out of the nucleus, hence only free
protons count.

###  1.4 Detectors

With IBD energy depositions in the range of a few MeV, the optical signal in
large detectors (tons) will not be terribly large. For economic reasons the
configuration of such instruments has generally been of a large volume of
liquid surrounded by optical sensors, generally photomultipliers. Smaller
detectors, such as employed in laboratory-scale double beta decay and direct
dark matter searches, utilize a variety of techniques such as crystals and
solid state detectors. The price of these options, however, limits their
practical use to detectors on the scale of tons. At the kiloton scale, the
only viable options at present are organic fluids or water. Fortunately, these
liquids may be “doped” with a variety of materials in small concentration to
enhance light output and spectral match to the light sensors.

We do not perform a detailed investigation of possible liquid scintillators
for this study, but rather we assume an organic liquid such as those employed
in the KamLAND and Borexino detectors, and that proposed for the 50 kiloton
LENA detector [64]. Of course, purified water costs far less but it produces
roughly 30 times less light (the organic liquid producing light from
scintillation and the water only from the Cherenkov radiation). For
comparison, the KamLAND scintillation detector produces signals of about 250
photoelectrons per million electron volts of ionization (PE/MeV) whereas the
water Cherenkov Super-Kamiokande [52] produces about 8-10 PE/MeV. Thus water
should have additives to increase light production, and such is under active
study at present. Therefore, we conservatively assume organic liquids in this
study.

Neutron absorption depends upon the choice of medium as well. In water and
organic liquids the IBD-produced neutron generally wanders for a hundred or
more microseconds before capture by a hydrogen nucleus. The deuterium formed
de-excites with emission of a 2.2 MeV gamma ray. The addition of an absorber
with a larger neutron capture cross section and a larger de-excitation energy
improves the measurement resolution by causing the neutron capture to be
closer in time and space to the neutrino interaction. A favorite material
considered for addition to water at present is gadolinium [92], which is
inexpensive, has water soluble salts and an emission of order 8 MeV.

###  1.5 Background Sources

There are several types of background sources which can cause false neutrino
signatures. Fortunately, the main discriminant for IBD is the spatial
coincidence (within a meter or so) and the temporal coincidence (within a few
microseconds) between the prompt (positron) and delayed (neutron capture)
signals, along with the prompt signal providing a known energy release.
Background sources which may fake neutrino signatures fall into two
categories, one which simulates both prompt and delayed signals, and one where
the prompt and delayed signals are of separate uncorrelated origin but
randomly occur in near coincidence. These latter events, the so called
“accidentals”, are generally restricted to the lower end of the energy range
since the frequency of single noise hits increases at lower energies. The
sources of background can be from internal residual radioactivity, cosmic rays
(only muons are important at depths greater than a few meters) passing through
the detector, cosmic rays outside the detector which make products such as
fast neutrons which enter the detector [51], and finally from radioactivity
outside the detector. The external radioactivity of concern is largely gamma
rays since alphas and electrons do not travel far.

Of course, there are neutrinos which come from other reactors around the
world, there are neutrinos from radioactive decays within the Earth, there are
cosmic ray neutrinos, and there are potentially neutrinos from beyond the
Earth. We deal with these sources one at a time, here qualitatively and later
numerically:

  1. 1.

Accidental Coincidences. It is difficult to make general statements about the
rate of random coincidences, since they depend entirely upon environmental and
detector construction details. However, based upon the experience of the two
extant large neutrino instruments KamLAND and Borexino, we know that this rate
can be kept to nearly ignorable levels. Much of the source of single pulses
will originate in the detector structure and environs, so there is a clear
benefit to larger detectors with smaller surface-to-volume ratios.

  2. 2.

Detector Internal Radioactivity. With care in detector preparation, the
neutrino community has learned to fabricate detectors with negligible levels
of internal radioactivity. The U and Th decays chains present a continuing
concern, and radon presents a dangerous and subtle pollutant. However, the
technology to reach extremely low levels for internal radioactivity has been
well demonstrated in the last several decades.

  3. 3.

Detector External Radioactivity. External radioactivity, such as from 40K only
causes problems near the detector surface. Two meters (water equivalent)
thickness of shielding is generally enough to reduce it significantly. Such a
background is easily revealed by the clustering of reconstructed event
locations that occurs near the surface.

  4. 4.

Penetrating Cosmic Ray Muons. Muons originating from cosmic ray interactions
high in the atmosphere (typically 20 km above) penetrate to the deepest mines,
having been observed two miles underground and under water. Since the incoming
cosmic ray spectrum falls steeply with energy ($\sim E^{-2.7}$) the energy of
muons falls as well, and the flux falls steeply with depth. A strong
constraint on detector depth arises simply from having a rate of muons so high
that the electronics become saturated if the depth is too shallow. We deal
with this quantitatively in this study. To give some sense of scale, current
large detectors need to be under more than about 2000 meters of water
equivalent (MWE) to have low enough rates to reduce electronics dead time to
negligible levels. These events are due to long lived 9Li and 8He isotopes,
and are generally referred to as “cosmogenic” background (see Section 3.3).

  5. 5.

External Cosmic Ray Interaction Products. Fast neutrons from muon interactions
in the detector surroundings pose a similar concern to entering external
radioactivity products. These neutrons can scatter elastically from protons
and then be absorbed, mimicking the prompt and delayed signals of a true IBD
event. Protection against such comes from possible shielding with muon
detectors in the outer layers. But in the case of large detectors, again,
these events are confined to the near surface and are discernible.

  6. 6.

Geo-neutrinos. The study of geo-neutrinos and their origins is an interesting
goal in itself, which we do not discuss at any length in this paper. We
carefully consider geo-neutrinos because they represent the largest background
for applied experiments (such as those related to reactors). The uranium and
thorium decay chains have energies extending up to 3.27 MeV. Hence, one common
way to reduce geo-neutrino background is simply to reject events with visible
energies above some threshold, but this approach can discard nearly half of
the desired (and never too strong) reactor signal. The rates of the geo-
neutrino events are low in current detectors, only recently being reported for
the first time by KamLAND [4] and Borexino [56].

If direction resolution improves in the future, then event reconstruction
(discussed below) can potentially identify geo-neutrinos events by isolating
those events arriving from below the horizon. This would dramatically improve
the signal-to-noise ratio for applied experiments. In such experiments, geo-
neutrino data can be harvested on a non-interfering basis to provide excellent
parallel science. In this paper we utilize all the data, employing all
possible information and a Bayesian approach to include the geo-neutrino data
in applied experiments.

  7. 7.

Atmospheric Cosmic Ray Produced Neutrinos. Fortunately for this subject, most
of the cosmic rays are of much higher energy, in the hundreds of MeV range
with average around a GeV. These events have a low enough rate not to be a
problem and are well separated in detected energy. With large enough
detectors, the study of these events themselves will become a useful and non-
interfering byproduct. See Honda et al. [106] for more information on
atmospheric neutrinos.

  8. 8.

Extraterrestrial Neutrinos. The search for extraterrestrial neutrinos, such as
those from ancient supernovae and from gamma ray bursts, has been a grand goal
for many in the neutrino field for decades. The only occurrence of such a
detection was during the Supernova 1987A, observing a total of 19 interactions
between the Irvine-Michigan-Brookhaven (IMB) detector [95] and Kamioka [96]
detectors. Nearby supernovae in our galaxy are rare (a few per century) and
thus pose no problem for neutrino experiments, but provide a potentially
wonderful opportunity for parallel science.

  9. 9.

Other Reactors. In Figure 11 we present a map of the calculated flux of
neutrinos from power and research reactors around the world. One cannot
currently distinguish a neutrino’s reactor-of-origin on an event by event
basis since all nuclear reactors share similar spectra (to a few percent) and
current detectors have nearly isotropic direction resolution. Statistically
however, neutrino oscillation modifies the reactor spectra and reduces the
flux (measured over a wide range of energy) at great distances. If one can
achieve sufficient angular measurement resolution, then one may reject events
based on arrival directions.

There have been a number of suggestions of a natural nuclear reactor, at the
Earth’s core [101] or perhaps at the inner core boundary [102] or at the core-
mantle boundary [103]. While geologists generally reject such models, they
have not yet been ruled out at the level below 3 TWth [56]. Such would provide
a significant, and detectable, background for a very large detector located
far from any nearby reactors.

###  1.6 Other physics and astrophysics

The list of possible scientific experiments that can be accomplished with the
sort of detector we are discussing is quite impressive, both in the fields of
elementary particle physics and in astrophysics. Study of these topics will
generate considerable interest in the scientific community and will result in
many theses, publications, and in the training of a cadre of professionals in
these fields.

The physics and astrophysics study list includes the following:

  1. 1.

Long baseline studies using a neutrino beam from an accelerator

  2. 2.

Proton decay searches, particularly in the Kaon decay modes, not well studied
by SuperKamiokande

  3. 3.

Search for supernova burst neutrinos, most likely from our galaxy

  4. 4.

Search for relic neutrinos summed from all past type II supernovae

  5. 5.

Atmospheric neutrinos, moving beyond SuperKamiokande

  6. 6.

Solar neutrinos, at least careful monitoring of solar output versus time

  7. 7.

Reactor neutrino properties, depending upon the proximity and power of
reactors

  8. 8.

Geo-neutrinos both as object of study and background, depending upon location
and as discussed below

  9. 9.

Pion decay at rest neutrinos from a possible portable low energy accelerator
(the DAEDALUS idea [98])

  10. 10.

Possible experiments employing radioactive sources for short range searches
for or studies of sterile neutrinos

These are already discussed in the specific framework of the LENA “whitepaper”
for a 50 kiloton liquid scintillation detector in Europe [64], and we shall
not elaborate further upon them here. First, it should be noted that all of
these are low duty factor studies which carry on without interference and in
parallel with applied experiments… they are essentially “free”. Some, such as
those involving a neutrino beam from an accelerator will depend upon
placement. Others, such as employing a radioactive source may require access
to the detector. For the most part, the important searches for proton decay,
studies of atmospheric neutrinos, and for rare astrophysical phenomena will go
along unimpeded wherever the detector is placed. The presence of such
detectors would serve to engage a large scientific community.

In the following we make explicit and conservative models of background
sources, in three dimensions, and we test the capabilities of the simulated
detectors in a variety of situations via simulation. We construct a
hypothetical detector in Section 2, which we call TREND (Test REactor
antiNeutrino Detector), and explicitly define its properties, resolution and
the expected background. In Section 3, we construct a thorough geospatial
model including the variation of all background sources with location and
depth. In Section 4, we present our optimal estimator theory and apply it. In
Section 5, we discuss the synergistic determination of neutrino properties in
terms of neutrino oscillation parameters. Finally in Section 6, we explore the
applied experiment of geolocating a reactor in several different detector
configurations111The choice of locations was designed to match those studied
previously by a French group [2], to explore different geographical
constraints on detector placements, and to sample different background
environments.. We finish with a summary of our findings and we highlight the
areas deserving further research.

##  2 Antineutrino Detector Model

The Detector Model determines an arbitrary detector’s neutrino measurement
resolution (both angle and energy) via Monte Carlo (MC) simulations. It also
determines the energy-dependent efficiencies associated with various candidate
selection criteria. The Detector Model is composed of three distinct parts: a
simulation of particles (and their energy depositions) produced by inverse
beta decay in liquid scintillator, a model for photon production and
propagation, and a detector hardware model for the photon detection process.
The Detector Model has been validated two ways, against CHOOZ experimental
angular resolution and against KamLAND experimental energy resolution.

Figure 1: $E_{\bar{\nu}_{e}}=5$MeV randomly placed inverse beta decay event
within TREND, about 100 ns after the event. The event vertex occurred about 4m
from the detector’s lower cylinder wall. $\bar{\nu}_{e}$ trajectory in green,
photons in red, PMT outlines in grey. Cylinder is 96m long and 46m in
diameter. The 1 m thick outer region of purified water for Cherenkov detection
of cosmic muons is not shown for purposes of visual clarity.

Determination of a detector’s measurement resolution is a precursor to the
placement of the detector within the Geospatial Model discussed in Section 3.
The resulting energy measurement resolution consists of an energy-dependent
table of standard deviations which correspond to the energy measurement
uncertainty probability density function (pdf). This uncertainty pdf is
assumed Gaussian. The resulting angle measurement resolution consists of an
energy-dependent table of angular signal-to-noise ratio (SNR) values. The
angular SNR is defined as the mean displacement between the prompt and delayed
signal divided by the direction vector uncertainty per Cartesian axis. A
detector’s measurement resolution is employed in two ways by the Geospatial
Model. First, it is used to add the correct amount of energy and angle
measurement noise to the measurements of detectors placed in the Geospatial
Model. Second, the detector resolution is passed to the maximum a posteriori
(MAP) estimator discussed in Section 4 where it is used to account for the
impact of measurement noise on parameter estimation.

###  2.1 Inverse beta decay model

The Detector Model begins with output from a GEANT [88] (Version 4.9.1 patch
3) simulation. The GEANT simulation provides the locations, timestamps,
kinetic energy and energy deposition of all of the particles and daughter
particles of the positron and neutron produced by the inverse beta decay
event. A large, homogeneous volume of linear alkyl benzene (LAB) scintillation
fluid was chosen as the target medium. Positrons and neutrons were directly
introduced into GEANT simulations in a manner consistent with the statistical
properties described by Vogel [30] while enforcing the appropriate
conservation laws (e.g. positron-neutron momentum along all three dimensions).

The GEANT dataset used in this paper consists of 10,000 sample IBD events at
10 different discrete neutrino energies from 2MeV through 11MeV. These 100,000
sample events were used in the Detector Model to characterize the measurement
resolutions of TREND, CHOOZ and KamLAND, assuming random placement of events
within each detector. It is important to note that the Geospatial Model does
not directly utilize these 100,000 sample events. Instead, it uses the energy-
dependent measurement resolution that was derived from these 100,000 sample
events. This allows statistical outliers (a real possibility) to occur in the
Geospatial Model even if such outliers were not included in the 100,000 sample
events.

###  2.2 Photon propagation model

The Photon Model (written in MATLAB [89]) accepts GEANT output and generate MC
photons via scintillation and Cherenkov radiation [91]. The model assumes an
amalgamation of parameters derived from various sources to model photon
scintillation, Cherenkov radiation, attenuation and re-emission. The
scintillation spectrum was defined by Eljen Technology’s EJ-254 1% boron-doped
plastic scintillator datasheet [75]. The scintillation spectrum has a peak at
421 nm, a scintillation yield of 9200 photons per MeV (deposited), and a
scintillation decay constant of 2.2 ns (fast component only). Photo-multiplier
tube (PMT) quantum efficiency (QE) curves were modeled upon a notional
absorption spectra which is peaked at 425 nm and extends down to 200 nm in the
ultraviolet region.

Oleg Perevozchikov’s thesis [77] was used to define an energy-dependent re-
emission efficiency curve, scaled to about 20% efficiency at 425 nm. The
wavelength-dependent attenuation length curve was taken from Lightfoot [76]
and scaled to give a scintillation-weighted mean of slightly over 20 m, per
Lasserre’s example [2]. The refractive index was defined by the Sellmeier
equation shown in Equation 1, scaled to provide a scintillation-weighted mean
index of refraction of about $n=1.58$.

$n(\lambda)=1+\frac{B_{1}\lambda^{2}}{\lambda^{2}-C_{1}}+\frac{B_{2}\lambda^{2}}{\lambda^{2}-C_{2}}+\frac{B_{3}\lambda^{2}}{\lambda^{2}-C_{3}}$
(1)

Here $B_{1}=1.44$, $C_{1}=81.76\mathrm{nm}^{2}$, $B_{2}=B_{3}=0$, and
$C_{2}=C_{3}=0\mathrm{nm}^{2}$.

All of the optical properties referenced above are displayed in Figure 2.
Figure 4 shows the photons created in one MC run compared to the optical
properties described above (and shown in Figure 2) for a neutrino inverse beta
decay event centered within TREND.

Figure 2: Photon model optical properties. For our emission spectrum we used
Eljen Technology’s EJ-254 1% boron-doped LS datasheet [75]. Oleg
Perevozchikov’s thesis [77] defines our energy-dependent re-emission
efficiency curve shape, scaled to 20% at 425nm. Attenuation length curve
derived from Lightfoot [76] and scaled according to Lasserre et al. [2].
Sellmeier equation fitted to refractive index curve to achieve a scintillation
wavelength-weighted mean of about $n=1.58$. Notional PMT QE curve set to peak
at 425nm and survive down to 200nm.

Poisson statistics were applied according to Equation 2 to create discrete
integer photon counts from GEANT’s continuous energy deposition values. In
Equation 2, $N$ is the number of photons produced at a single energy
deposition point.

$N=\mathcal{P}\left(\bar{N}\right)$ (2)

$N$ is a Poisson random number, sampled from the mean photon count, $\bar{N}$,
defined in Equation 3

$\bar{N}=y(kB)dE$ (3)

where $dE$ is the energy deposition given by GEANT (in MeV), $y$ is the
scintillation yield (in photons/MeV) and $kB$ is the “Birk’s quenching
factor222True quenching factor values vary by scintillator and are a function
of particle ionization state, kinetic energy, and mass energy. These values
are typically experimentally determined, but notional values were assumed for
simplicity, independent of energy.”, a ratio of visible energy released to
total energy released. In the case of light ionizing particles such as
electrons and positrons a value of $kB=1.00$ was assumed. In the case of heavy
ionizing particles such as protons a value of $kB=0.10$ was assumed.

Charged particles exceeding the local speed of light in the medium produce
Cherenkov radiation. Equation 4 defines the scaled Cherenkov emission
spectrum.

$f(\lambda)=\left\\{\begin{array}[]{l l}\beta\geq\frac{1}{n(\lambda)}&\quad
2\pi\alpha q^{2}\mu\left(\frac{1-(\beta
n(\lambda))^{-2}}{\lambda^{2}}\right)\times 10^{6}\\\
\beta<\frac{1}{n(\lambda)}&\quad 0\\\ \end{array}\right.$ (4)

Here, $\beta=\frac{|\vec{v}|}{c}$, is the ratio of the particle speed through
the fluid to the speed of light in a vacuum, $n(\lambda)$ is the unitless
fluid index of refraction at wavelength $\lambda$ (in nanometers), $\mu$ is
the relative permeability of the fluid (assumed to be 0.999992), $q$ is the
particle charge in units of elementary charge, and the fine structure constant
$\alpha=1/137.036$.

The units of $f(\lambda)$ are the mean number of Cherenkov photons produced
per mm of particle travel distance per nanometer of photon wavelength. Figure
3 plots this value across a typical $\lambda$-$\beta$ domain.

Determination of the (Poisson mean) number of Cherenkov photons produced per
mm from $\lambda_{1}$ to $\lambda_{2}$ requires an integration of Equation 4
across the $\lambda_{1}$ to $\lambda_{2}$ spectrum:

$\bar{N}=\int\limits_{\lambda_{1}}^{\lambda_{2}}f(\lambda)d\lambda$ (5)

Cherenkov radiation “cone” half-angles $\theta_{c}$ were assumed to be
continuously variable over the portion of each charged particle’s trajectory
where $\beta\geq\frac{1}{n(\lambda)}$. These angles, $\theta_{c}$, defined in
Equation 6 and shown in Figure 3, are a function of a Cherenkov photon’s
wavelength $\lambda$ (in nm) and parent particle’s speed $\beta$.

$\theta_{c}(\lambda)=\arccos\left(\frac{1}{\beta n(\lambda)}\right)$ (6)

The parent particle speed $\beta$ may be defined as a function of its kinetic
energy KE and rest mass $m_{0}$ by Equation 7.

$\beta=\frac{\sqrt{\mathrm{KE}^{2}+2m_{0}c^{2}\mathrm{KE}}}{\mathrm{KE}+m_{0}c^{2}}$
(7)

Figure 3: On the right panel, Equation 4 is evaluated: Cherenkov photon
production per mm of particle path-length as a function of photon wavelength
$\lambda$ and parent particle speed $\beta$. On the left panel, Equation 6 is
evaluated: Cherenkov angles $\theta_{c}$ as a function of wavelength $\lambda$
and parent particle speed $\beta$.

A single-exponential time delay (2.2 ns) was applied to all scintillation
photon emissions, and photon “start” times were randomly sampled from this
exponential distribution. We expect effects of a second slower exponential to
be small and ignore them in the present simulation. Photon attenuation and re-
emission was also modeled. Wavelength dependencies including attenuation
length, re-emission fraction, refractive index and QE were defined for each
individual photon as a function of each photon’s wavelength333No two photons
share the same exact wavelength (or the same exact wavelength-dependent
optical properties) as the scintillation spectrum is treated as continuous;
i.e. it is randomly sampled using a linear interpolant between defined table
lookup points..

Photons which are absorbed and then re-emitted are probabilistically assigned
a new wavelength using a re-emission spectrum which is different than the
original scintillation spectrum. For each such photon, the re-emission
spectrum is scaled to ensure that energy conservation is not violated (i.e.
the new wavelength is longer than its predecessor). Photon re-emission
directions are assumed isotropic for all re-emitted photons regardless of
their properties or source of origin. This causes photons to gradually lose
directional “information”444The term “information” here denotes any
information which may contribute to a mathematical reconstruction of an
event’s location and $dE$. over time as they repeated absorptions and re-
emission. This is especially true of high energy Cherenkov photons (see Figure
2) because of the shorter attenuation lengths in the ultraviolet/blue region
of the spectrum.

Figure 4: Photon histograms from one TREND MC IBD event, using LS optical
properties from Figure 2. Blue histograms represent MC photons. Red lines are
estimator-assumed emission spectra. Green lines are the estimator-assumed
observed spectra, whose $\lambda$-weighted mean values are used (since PMTs
are “colorblind”) to estimate the prompt and delayed “bookends” (positions and
energies). Deviation between the emitted and observed spectra are caused by
operation of the QE curve on all observed photons.

###  2.3 Detector hardware model

Any arbitrary detector design can be selected for use with the Inverse Beta
Decay Model (described in Section 2.1) and the Photon Model (described in
Section 2.2) to quickly test different detector designs of varying volume on
the same GEANT dataset. The Detector Hardware Model consists primarily of a
physical description of the size and shape of the scintillating volume and its
enclosure(s), the size and shape of the photo-detectors, the optical coverage
due to photo-detectors, and the wavelength-dependent QE curve for the photo-
detectors. The physical structure of TREND (the detector selected for this
study) was largely adopted from the Secret Neutrino Interaction Finder (SNIF)
detector proposed by Lasserre et al. [2]. We are, however, proposing to
operate TREND quite differently. We propose a different set of selection
criteria, to operate without a low energy cut, to utilize energy and angle
measurements, to use a different estimation technique, and to operate at
deeper depths. In future work, we hope to further optimize TREND using the
flexibility of the Detector Hardware Model.

In addition to TREND, Regions I and II of CHOOZ were also modeled for angular
measurement resolution validation, as well as KamLAND for energy measurement
resolution validation. Our Detector Model incarnations of these two real world
detectors are shown in Figure 5. As evidenced in Figure 5, PMTs were modeled
as flat squares lying tangent to the detector wall. PMT temporal response
functions and wavelength-dependent QE curves were imported from manufacturer’s
datasheets. Rayleigh scattering, and photon reflection at the LS/PMT boundary
are not modeled555Concerns over the omission of Rayleigh scattering and photon
reflection at the LS/PMT boundary may be allayed by the positive validation
results for CHOOZ and KamLAND presented in Section 2.4..

For simplicity both CHOOZ and KamLAND were modeled as cube-shaped in our
Detector Model, though we did model TREND as a cylinder. We notice almost no
adverse effects in our validations due to simplifying the detector spherical
geometries into a cube shapes.

Figure 5: CHOOZ and KamLAND models as created by the Detector Model. Each
detector (modeled as a cube rather than a sphere for simplicity) shows the
same GEANT IBD event at detector center. Purple PMTs (square “pixels” lying on
the outer region walls) indicate active voltage jumps from photon detection
(the images were taken about 500 ns after $\bar{\nu}_{e}$ annihilation). Note
separate inner regions (thick grey lines) and outer regions (thin grey lines).

###  2.4 Validation results

To establish confidence in the measurement resolution values that the Detector
Model predicts for TREND, two separate validation tests were carried out. The
first against experimental KamLAND energy resolution [74] (Section 2.4.1), and
the second against experimental CHOOZ direction vector resolution [3] (Section
2.4.2).

####  2.4.1 KamLAND energy resolution validation

The Kamioka Liquid-scintillator Anti-Neutrino Detector (KamLAND) is located in
the Kamioka mine, which is approximately 1 km ($\sim$2700 MWE) below the
summit of Mt. Ikenoyama, Gifu prefecture, Japan (36.43∘ N, 137.31∘ E). The
presently operating detector measures electron antineutrinos from nuclear
reactors (Eguchi et al. [41]) and the earth (Araki et al. [5]), using
$\sim$1000 tons of ultra-pure scintillating liquid monitored by $\sim$1900
PMTs providing $\sim$34% photo-coverage. The detector realizes an energy
resolution of $6.5\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$ (Abe et al. [51]).

KamLAND was simulated with a refractive index of 1.44, a (scintillation
wavelength-weighted) mean attenuation length of 8.46m, a LS exponential decay
constant of 6ns, a yield of 9200 photons/MeV, (scintillation wavelength-
weighted) mean quantum efficiency of about 29%, a (scintillation wavelength-
weighted) mean re-emission efficiency of about 13%, a 1150m3 inner region and
a 3050m3 outer region with 324 PMT’s per face (1944 total). A fixed quenching
factor of 0.1 was assumed for heavy ionization particles, and the
scintillation spectrum was lifted from our EJ-254 datasheet [75]. Each PMT was
assumed capable of 1 ns timing resolution, with no flat-fielding effects. Each
PMT was modeled as a square, 0.47m per side, separated from its neighbor by
0.81 m (center-to-center), providing 34% detector surface area coverage. An
earlier alternate configuration is also modeled with 22% PMT coverage. This
earlier configuration uses identical fluid properties as above but has 225
PMT’s per face (1350 total). Each PMT in this configuration was modeled as a
square, 0.45m per side, separated from its neighbors by 0.97 m (center-to-
center), providing 22% total surface area coverage.

KamLAND energy validation results are shown in Table 1. Note that these
validation results pertain specifically to center-detector point-sources of
energy, not IBD sources. A full comparison between IBD energy resolution and
point-source energy resolution, as well as center vs off-center event
resolution, is presented in Section 2.4.3.

The full comparison over energy can be seen in Figure 6. Our Detector Model
produced very similar resolution values across the energy spectrum when
compared to actual KamLAND experimental data, both before and after the
KamLAND PMT additions of 2003. This provides us with high confidence in the
$\left.8.9\%\right|_{E_{\mathrm{vis}}=1\mathrm{MeV}}$ energy resolution that
our Detector Model produced for TREND in Table 4.

Comparison Metric | Cited [74] | Our MC
---|---|---
Energy Resolution $1\sigma$ | $6.2\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$ | $6.24\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$
(34% PMT coverage) |  | (82.4% candidates)
Energy Resolution $1\sigma$ | $7.3\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$ | $7.69\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$
(22% PMT coverage) |  | (82.0% candidates)
Table 1: KamLAND $\bar{\nu}_{e}$ visible-energy resolution validation results
for center-detector point sources of energy. Values in right column represent
nonlinear least squares fits to our KamLAND MC energy estimate errors. We
compared the KamLAND detector before and after the 2003 additional PMT
installation cited by [74]: “The central detector PMT array was upgraded on
February 27, 2003 by commissioning 554 20-inch PMTs, increasing the photo-
cathode coverage from 22% to 34% and improving the energy resolution from
$7.3\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$ to
$6.2\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$.” Our KamLAND “candidates” met
all the criteria specified for CHOOZ candidates in Table 2. Figure 6: KamLAND
validation results for center-detector point sources of energy. Our Detector
Model MC results (red and blue points) are shown vs KamLAND [74] cited values
(red and blue solid lines). Blue results show 34% PMT coverage validation, red
results show 22% PMT coverage validation.

####  2.4.2 CHOOZ direction resolution validation

The CHOOZ detector was located $\sim$1 km from the 4.25 GWth nuclear reactor
near the village Chooz, along the River Meuse in northern France. To reduce
the cosmic ray muon flux the detector was placed underground with an
overburden of $\sim$300 MWE. The active target of the detector consisted of 5
tons of Gd-loaded scintillating liquid. Scintillation light was collected by
$\sim$200 PMTs, achieving $\sim$15% photo-coverage and an energy resolution of
$9\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$ (Apollonio et al. [3]). This
detector measured the direction of reactor antineutrinos with a 1$\sigma$
resolution of 18∘.

CHOOZ was simulated with a refractive index of 1.47, a (scintillation
wavelength-weighted) mean attenuation length of 3.4m, a LS exponential decay
constant of 7ns, a yield of 9200 photons/MeV, (scintillation wavelength-
weighted) mean quantum efficiency of about 29%, a (scintillation wavelength-
weighted) mean re-emission efficiency of about 13%, a 5.6m3 inner region and a
22m3 outer region with 6 PMT’s per face (24 total). A fixed quenching factor
of 0.1 was assumed for heavy ionization particles, and the scintillation
spectrum was lifted from our EJ-254 datasheet [75]. Each PMT was assumed
capable of 8 ns timing resolution, with no flat-fielding effects. Each PMT was
modeled as a square, 0.54m per side, separated from its neighbor by 1.40 m
(center-to-center), providing 15% detector surface area coverage.

With CHOOZ, as with TREND, the direction of the incoming neutrino is
calculated by forming a vector between the estimated location of neutron
capture (the delayed signal) and the estimated location of the positron
annihilation (the prompt signal). Because the neutron can bounce many times
before its capture, the neutrino directional information recorded for each IBD
event is typically very poor using this method. Because of this CHOOZ (and
others) often quote their direction resolution as a spherical 1$\sigma$
uncertainty after a fixed large number of IBD events. In this work, we
introduce a metric called “angular SNR” which quantifies the directional
information contained per single IBD event. The angular SNR (SNR$\approx$.07
in CHOOZ), defined in Equations 53 \- 54, is the mean positron-neutron
displacement ($\sim$15mm in CHOOZ) divided by the norm ($\sim$210mm) of the
combined neutron ($\sim$200mm) and positron ($\sim$60mm) location
uncertainties.

CHOOZ direction validation results are shown in Table 2. The Detector Model
values compared well to the CHOOZ experimental values across all the different
comparison metrics, including candidate selection cuts, mean neutron
displacement, neutron position resolution, and the primary metric of interest,
direction vector SNR. This successful validation establishes grounds for
confidence in the 0.045 SNR direction vector resolution our Detector Model
produced for TREND in Table 4

Comparison Metric | Cited [3] | Our MC
---|---|---
positron energy cut | 97.8% | 99.9%
positron-geode distance cut | 99.9% | 99.7%
neutron capture cut | 84.6% | 82.5%
capture energy containment cut | 94.6% | 94.7%
neutron-geode distance cut | 99.5% | 97.1%
neutron delay cut | 93.7% | 99.8%
positron-neutron distance cut | 98.4% | 99.9%
neutron multiplicity cut | 97.4% | 100%*
combined cut | 69.8% | 75.1%
Mean Neutron Displacement | 19mm (17mm sim) | 16.6mm
Neutron Position Resolution $1\sigma$ | 190.0mm | 181.9mm
Angular $1\sigma$, $n$=2500 events | 18.0o (19.0o sim) | 18.9o
SNR | 0.100 (.091 sim) | 0.091
Table 2: CHOOZ direction resolution validation results. Energy dependent
values weighted to a 1km range reactor spectrum. *Neutron multiplicity cut not
implemented because our GEANT dataset was comprised of single-neutron inverse
beta decay events only. CHOOZ own simulation results [3] shown in parenthesis
next to CHOOZ experimental results [3].

####  2.4.3 TREND energy and direction resolution MC results

Once validated against CHOOZ and KamLAND experimental results, the Detector
Model was tasked with determining the energy and direction resolution of the
proposed TREND detector. TREND “candidate” events were selected to meet all
the CHOOZ candidate criteria specified in Table 2 with the exception of three
TREND-specific modifications:

  1. 1.

Neutron delay cut increased from 100$\mu$s in CHOOZ to 200$\mu$s in TREND.

  2. 2.

Positron energy cut increased from $E_{\mathrm{vis}}^{e^{+}}<$8MeV in CHOOZ to
$E_{\mathrm{vis}}^{e^{+}}<$12MeV in TREND.

  3. 3.

Positron-neutron distance cut increased from 1m in CHOOZ to 2m in TREND.

We believe the time and distance increases above were warranted by the poor
position resolution and vastly larger size of TREND compared to CHOOZ. The
positron energy cut increase was implemented to capture reactor antineutrinos
to the very end of their detectable energy spectrum (near 11MeV). Figure 7
shows how increasing the positron-neutron distance cut from 1m to 2m lets us
accept a much larger percentage of neutrinos as candidates, while increasing
the neutron delay cut from 100$\mu$s to 200$\mu$s has a similar, but smaller,
effect.

For TREND, the scintillator optical properties assumed by the Photon Model
were defined by the various sources cited in Section 2.2. The modeled TREND
detector (shown intercepting a $\bar{\nu}_{e}$ in Figure 1) is a cylinder 96m
high with a radius of 23m, comprising 16906 PMTs of 0.2m2 area each to provide
almost 20% optical coverage. The fiducial mass of 138,000 metric tons brings
the free proton count to $\sim 10^{34}p^{+}$. Table 3 provides a complete list
of parameter inputs used by the Detector Model to construct TREND, including
the candidate cut criteria used to produce the statistical results presented
in Table 4.

Figure 8 shows all 16906 square-shaped PMTs in the TREND detector model, as
seen from the perspective of an observer placed at the $\bar{\nu}_{e}$ vertex
shown in the Figure 1 IBD event. Summing over the angular area within the
pixels shown in Figure 8 can provide a useful measure of the solid angle
coverage from the IBD vertex. The solid angle coverage is naturally a function
of the simpler TREND surface area666Note that even though the TREND PMT
surface area coverage is about 19.9%, the PMT solid angle coverage will vary
slightly depending on the observer’s location within the TREND detector. As an
example, solid angle coverage from the $\bar{\nu}_{e}$ vertex near the
detector wall shown in Figure 1 is about 20.4%, while the solid angle coverage
as seen by an observer exactly in the center of the detector is 19.9%, more
closely matching the surface area coverage. coverage (19.9%), but is a more
accurate measure of how much energy the detector can expect to capture from a
single IBD event.

Figure 7: $E_{\bar{\nu}_{e}}=5$MeV MC statistics, 10,000 $\bar{\nu}_{e}$
events. Prompt “bookends” in red, delayed “bookends” in blue. Select candidate
criteria are also shown in this figure: the 200$\mu$s time cut (99.8%
candidates), 2m range cut (96.5% candidates), and
4MeV$<E^{n^{0}}_{\mathrm{vis}}<$12MeV neutron energy cut (83.3% candidates).

A third very useful metric, the energy capture efficiency, or ECE, is also
shown in Figure 8. The ECE includes attenuation effects from the fluid, and
can predict the total fraction of energy that a neutrino will deposit on the
detector PMTs. In this respect it is a much more useful metric than surface
area coverage or even solid angle coverage. The ECE will naturally vary as a
function of an event’s location within the detector. Comparing the ECEs of the
two vertices shown in Figure 8, we can see that a detector-center event in
TREND deposits only about 5% of its energy at the actual PMT surfaces, while
the example event near the detector wall deposits much more energy at the PMT
surfaces, about 12.5%. This indicates to us that off-center events are
actually preferred to detector-center events, as they will yield greater
photon collection per their greater ECE, and they should provide us with
consequentially better energy and position resolution.

100,000 events were used for the TREND MCs, 10,000 each at $\bar{\nu}_{e}$
energies of 2meV through 11MeV. Some example histograms produced from the
10,000 5MeV MCs can be seen in Figure 7. Application of select candidate
“cuts” are also shown in Figure 7, including the neutron delay cut, neutron
energy cut and positron-neutron distance cut.

Figure 8: 2D PMT solid angle projections for center vertex and off-center vertex cases. The off-center vertex case in the upper plot corresponds to the IBD event shown in Figure 1. Different colors indicate six different PMT regions: front (green, 3456 PMTs), top (blue, 3456 PMTs), bottom (purple, 3456 PMTs), back (pink, 3456 PMTs), right circular cap (aqua, 1541 PMTs) and left circular cap (yellow, 1541 PMTs). The effective energy capture efficiency values shown include attenuation effects inside the detector. The mean attenuation length modeled for our LS was 21.1m. TREND has a fiducial-volume radius of 23m and a height of 96.5m, so photon attenuation plays a significant role, harming detector-center events more than detector-edge events. TREND Parameter | Value
---|---
Cylinder Radius | 23.0m
Cylinder Height | 96.5m
Volume |  $\sim$160,000m3
Mass | 138kT
Protons | $10^{34}\mathrm{p}^{+}$
Attempted PMTs | 17000
Achieved PMTs | 16906
Optical Coverage | 19.9%
PMT Time Resolution |  $<1$ns
PMT Area | .203m2
PMT Shape | flat square
Mean PMT QE | 23.1%
Peak Absorption Wavelength | 425nm
GEANT Scintillator | 0.5% Gd doped LS
Peak Emission Wavelength | 421nm
Yield | 9200photons/MeV
Scintillator Decay Constant | 2.20ns
Mean Refractive Index | 1.58
Mean Attenuation Length | 21.1m
Mean Re-Emission Fraction | 13.2%
Measured Cherenkov Fraction | $\sim 2\%$
Positron energy cut |  $E^{e^{+}}_{\mathrm{vis}}<12$MeV
Positron-wall distance cut |  $||x_{e^{+}}-x_{\mathrm{wall}}||>30$cm
Neutron capture cut |  $12$MeV$>E^{n^{0}}_{\mathrm{vis}}>4$MeV
Capture energy containment cut |  $E^{n^{0}}_{\mathrm{vis}}<6$MeV
Neutron-wall distance cut |  $||x_{n^{0}}-x_{\mathrm{wall}}||>30$cm
Neutron delay cut |  $t_{n^{0}}-t_{e^{+}}<200$ns
Positron-neutron distance cut |  $||x_{e^{+}}-x_{n^{0}}||<2$m
Neutron multiplicity cut | num$\left(n^{0}\right)=1$
Table 3: Assumed TREND characteristics used in TREND detector MCs. All candidate cut percentages shown are energy weighted mean values for a 1km range reactor spectra. Metric | TREND MC Results
---|---
Positron energy cut | 99.9%
Positron-geode distance cut | 97.1%
Neutron capture cut | 83.3%
Capture energy containment cut | 98.8%
Neutron-geode distance cut | 97.8%
Neutron delay cut | 99.8%
Positron-neutron distance cut | 96.5%
Neutron multiplicity cut | 100.0%
Combined candidate cut | 78.5%
Mean Neutron Displacement |  $20.7$mm
Positron Annihilation Point Resolution $1\sigma$ |  $374.5$mm
Neutron Capture Point Resolution $1\sigma$ |  $243.7$mm
Direction Vector Resolution $1\sigma$ |  $461.3$mm
Angular Direction $1\sigma$, n=2500 events | 37.0o
SNR | 0.045
Energy Resolution $1\sigma$ | $\left.8.9\%\right|_{E_{\mathrm{vis}}=1\mathrm{MeV}}$
Table 4: TREND MC results, energy dependent values weighted to a 1km range
reactor spectrum. Candidate events used in these statistics meet all the CHOOZ
candidate criteria specified in Table 2 with three TREND specific
modifications: neutron delay cut (200$\mu$s from 100$\mu$s), positron energy
cut ($E_{\mathrm{vis}}^{e^{+}}<12$MeV from $E_{\mathrm{vis}}^{e^{+}}<8$MeV),
and positron-neutron distance cut (2m from 1m).

The mean TREND energy resolution based on 100,000 randomly placed inverse beta
decay events (shown in Figure 9) was found to be about
$\left.8.9\%\right|_{E_{\mathrm{vis}}=1\mathrm{MeV}}$ at Evis=1MeV. This
random-vertex IBD energy resolution is not directly comparable to the KamLAND
center-detector point-source energy resolution, though we can directly compare
center-detector point-source energy resolutions. In TREND we find a center-
detector point-source energy resolution of 10.1%, about 50% worse than the
cited KamLAND energy resolution of
$7.3\%/\sqrt{E_{\mathrm{vis}}[\mathrm{MeV}]}$. TREND surface area coverage is
similar to KamLAND, however in TREND fewer photons reach the detector walls
due to the inevitably larger travel distances mandated by a 138kT detector,
which the increased attenuation length in TREND apparently does not completely
make up for.

Figure 9: TREND Detector energy resolution results. Three energy resolutions
are shown: 1. Center-detector point-source resolution (red), 2. Random-
location point-source resolution (green), and 3. Random-location
$\bar{\nu}_{e}$ energy resolution (blue). $\bar{\nu}_{e}$ statistics were
harvested from candidate events out of 10,000 possible GEANT neutrino IBD
events at integer energy levels of 2MeV through 11MeV (100,000 total events).
CHOOZ maximum likelihood estimator (MLE) [3] was used to estimate energies and
locations within TREND Detector Model for all three types of sources. Energy
estimate bias can be seen in left plot, useful for detector calibration.

Figure 9 shows the complete TREND energy resolution picture, including
measurement bias on the left panel. Estimation biases in this context usually
arise from a mismatch between the truth model and the estimator assumptions;
in our case this may be caused by several real-world effects occurring in the
Detector Model which are absent from the CHOOZ MLE, such as photon re-
emission. Such biases would scale almost linearly with $\bar{\nu}_{e}$ energy;
and indeed in Figure 9 we do observe bias to be a linear function of
$\bar{\nu}_{e}$ energy.

Determination of measurement bias in neutrino physics is useful for detector
calibration to remove such biases from future energy measurements.
Accordingly, in this paper we assume that all such biases have been removed
from energy measurements and we concern ourselves solely with measurement
noise. Three energy resolutions are shown in Figure 9: 1. Center-detector
point-source resolution (red), 2. Random-location point-source resolution
(green), and 3. Random-location $\bar{\nu}_{e}$ energy resolution (blue). This
figure substantiates our claim that detector-centered events have worse
resolution (at least for low-energy events) in very large detectors.

Figure 9 also allows us to directly compare IBD energy resolution with point-
source energy resolution. Interestingly, we see that IBD energy resolution and
point-source energy resolutions are comparable up until about Evis=3MeV, at
which point they begin to diverge as the IBD energy resolution suffers to a
greater degree than the point-source energy resolution. We believe this
deviation in resolution at the higher $\bar{\nu}_{e}$ energies to be a product
of several factors:

  1. 1.

Cherenkov radiation. Cherenkov radiation is directional in nature (produced in
significant quantities by the positron streak and any “knockoff” electrons),
and not suitably compensated for by the simple CHOOZ MLE estimator, which
assumes isotropic radiation.

  2. 2.

Extended source vs point source. IBD events are, naturally, extended sources
of radiation, a conglomeration of many smaller point-sources of radiation
produced by the IBD positron streak, gammas, and neutron-proton collisions,
among others. This extended-source nature violates the point-source assumption
of the simple CHOOZ MLE estimator.

  3. 3.

Prompt gamma energy floor. A portion of an IBD event’s energy is not a
function of $\bar{\nu}_{e}$ energy, but rather a constant value (an energy
“floor”) of energy release produced by two prompt gammas created upon positron
annihilation. Each gamma releases a fixed 511keV of energy, and this radiation
is not a function of the original antineutrino energy, as the simple
$kE^{1/2}$ fit equation seeks to provide. Thus we would not expect IBD events
to conform to the standard energy resolution fit equation either. This is the
reason that we have fit our IBD MC events with some additional constants
rather than trying to apply the usual $kE^{1/2}$ fit.

Direction vector resolution for TREND was found to be about SNR=0.045
(weighted for the neutrino energy spectrum of a reactor at a distance of 1km),
or 21mm/461mm, from a mean neutron capture point resolution of 243mm and a
mean positron annihilation point resolution of 374mm, which combine to form a
reconstruction vector resolution of 461mm per cartesian axis.

Note that the reconstruction vector resolution 461mm $1\sigma$ noise value
includes additive noise incurred by neutron random-walk natural to inverse
beta decay, and is not simply the norm of the 374mm positron CE $1\sigma$ and
243mm neutron CE $1\sigma$ values777A simple norm of the 374mm positron CE
$1\sigma$ and 243mm neutron CE $1\sigma$ values would produce a false smaller
reconstruction vector resolution $1\sigma$ of only 446m..

##  3 Geospatial Model

The Geospatial Model calculates the expected rate of detections within the
scintillator volume ($\bar{n}$) from all modeled sources at a specified
detector location. This rate includes neutrinos from reactors and the Earth,
as well as non-neutrino background which survives all selection criteria
applied, and accounts for the modeled energy-dependent detector efficiency and
veto-related duty cycle losses. The model also computes the elevation-azimuth-
energy measurement probability density function (pdf) for each detector,
smeared by the modeled measurement resolution. For simplicity, the measurement
resolution found for neutrino events were applied to all detections, including
those arising from non-neutrino background. The total number of detections in
each detector, $n_{z}$, was randomly determined according to Poisson
statistics. The resulting noisy measurements were randomly selected from the
appropriate measurement pdf and fed to the estimator detailed in Section 4 for
parameter estimation.

This process is repeated one thousand times to construct an a posteriori
parameter space pdf. The pdf defines the location and thermal power
observability of an unknown neutrino source or, in the case of oscillation
parameter estimation, the oscillation parameter observability. MC results for
an unknown neutrino source are presented in Section 6 and oscillation
parameter estimation results are presented in Section 5.

The Geospatial Model operates on a three dimensional representation of the
Earth. It uses the National Oceanic and Atmospheric Administration (NOAA)
Earth TOPOgraphical 1 (ETOPO1) “ice” data [83] to represent the land/ice
surface and the ocean bathymetric data, relative to the World Geodetic System
84 (WGS84) ellipsoid. The ETOPO1 dataset consists of 233 million tiles in a
10800 by 21600 matrix, providing global elevation resolution of 1 arc-minute.
For the oceans, the “surface” was taken to be the zero-tide elevation
determined via the National Geospatial-Intelligence Agency (NGA) Earth Gravity
Model 2008 [84] (EGM2008).

Several different coordinate systems are employed by the Geospatial Model.
Detections are measured in the local North-East-Down (NED) reference frame.
This reference frame is different for each detector in a given scenario. For
estimating the location of an unknown neutrino source, the parameter space is
expressed in the WGS84 Earth centered, Earth-fixed reference frame. This
reference frame can be expressed in (Latitude, Longitude, Height) or in a
Cartesian coordinate system that rotates with the Earth. When the unknown
source is constrained to the local terrain height, the parameter space reduces
to just two local dimensions, WGS84 Latitude and Longitude. The standard
Direction Cosine Matrix (DCM) defined in Equation 11 is used to rotate between
the NED and ECEF reference frames. Note that this rotation depends on the
detector’s geodetic latitude $\phi$ and longitude $\lambda$.

$\displaystyle
C_{NED_{d}}^{ECEF}=\left[\begin{array}[]{ccc}-\cos\lambda\sin\phi&-\sin\phi\sin\lambda&\cos\phi\\\
-\sin\lambda&\cos\lambda&0\\\
-\cos\phi\cos\lambda&-\cos\phi\sin\lambda&-\sin\phi\end{array}\right]$ (11)

###  3.1 Reactor Neutrinos

Neutrinos are detected in TREND by the inverse beta decay (IBD) reaction,
where a $\bar{\nu}_{e}$ interacts with a free proton in the liquid
scintillator (12).

$\bar{\nu}_{e}+p\rightarrow e^{+}+n$ (12)

In this study, the reactor energy spectrum is approximated as an exponential
fall-off with respect to a 2nd order polynomial in neutrino energy (see
Equation 48 in Section 4). The scaling of the total neutrino flux assumes
approximately 200 MeV and 6 neutrinos per fission, with an average of two
neutrinos per fission created above the inverse beta decay energy threshold of
$E_{\bar{\nu}_{e}}\approx 1.8$MeV. These assumptions yield $1.872\times
10^{20}$ detectable neutrinos per GWth per second emitted isotropically from a
reactor. This rate is comparable to that produced by a typical pressurized
water reactor at the beginning of its fuel cycle [68]. The TREND yearly
predicted $\bar{\nu}_{e}$ inverse beta decay observation rate for a 300MWth
source is shown in Figure 10 as a function of source range.

Figure 10: Number of predicted $\bar{\nu}_{e}$ inverse beta decay events
observed over 1 year (100% duty cycle) by a $10^{34}p^{+}$ detector, as a
function of distance from a 300 MWth reactor. This assumes 100% efficiency for
both the detector and the reactor, perfect measurements of neutrino energy, no
energy cut, and no background contribution.

For the longer ranges and smaller reactors in this study, it would be
difficult to extract much information about fuel mix from the spectrum. Thus
long range detection of reactors should focus upon detection, location, and
estimation of the time-averaged power output. If there is sufficient signal,
it might be possible to detect interruptions in the reactor operation over
time scales of months and crudely monitor long term power output.

The 1999 International Atomic Energy Agency (IAEA) data on nuclear power
reactors [85] was used to model the location and power of known reactors. The
data include 436 active reactor cores distributed among 206 locations (total
1063GWth) and 35 reactor cores (total 93GWth) under construction. The reported
electrical power was converted to thermal power assuming an efficiency of 33%,
regardless of the reactor design. A reactor duty cycle of 80% was assumed and
applied in the Geospatial Model. The neutrino event rate per 1032 protons per
year (assuming 100% efficiency) due solely to these active reactor cores can
be seen in Figure 11.

Figure 11: IAEA known reactor background for a $10^{32}\mathrm{p}^{+}$
detector, saturated at 100 events per year.

###  3.2 Geo-neutrinos

The underlying interior structure and composition of the Earth is, in some
regards, still poorly understood. The concentration and distribution of radio-
isotopes, whose decay chains are capable of producing significant neutrino
flux, does not escape this uncertainty. Therefore modeling of the
distribution, energy spectra, and total flux of geo-neutrinos remains a
challenging task on its own. The Geospatial Model does not attempt to be a hi-
fidelity model for geo-neutrino research, but rather it provides for a
practical representation of this complex and significant background neutrino
source.

The modeling of this flux was broken down into the mantle and the crust. The
Earth’s core was assumed to have no significant contribution to neutrino flux
(i.e. no geo-reactor) in this analysis. For the mantle (radii from 3480km to
6291km), the spherically symmetric density profile in the Preliminary
Reference Earth Model [86] (PREM) was used. For elemental abundances, the two-
layer stratified model suggested by Fiorentini et al. [69] was used. For the
mantle above a depth of 670 km, the elemental abundances were 6.5 ppb, 17.3
ppb, and 78 ppm for 238U, 232Th, and 40K respectively. Below this, the
abundances were 13.2 ppb, 52 ppb, and 160 ppm, respectively. The mantle
isotope abundances were assumed to be constant with values of approximately
99.3%, 100%, and 0.01% for 238U, 232Th, and 40K respectively.

The Geospatial Model uses CRUST 2.0 [70] to describe the Earth’s crust. This
model consists of seven layers (to which we add an 8th) of $2^{\circ}\times
2^{\circ}$ tiles beginning at the ETOPO1 surface and typically descending to
depths of 30km-70km. The lowest points of the 7th layer tiles range from
6302km Earth radius at the poles to 6366km Earth radius near the equator. A
spherical mantle is assumed in the Geospatial Model, and seamlessly joined to
the CRUST 2.0 by introducing an 8th “Mantle Adjoining” layer with properties
identical to the upper mantle. Figure 12 shows how the thickness of the CRUST
2.0 and Mantle Adjoining Layer tiles varies across the Earth. Typical CRUST
2.0 tiles are a few km thick, while some layers (like 1. Ice) are primarily
made up of zero-thickness tiles across wide swaths of the Earth’s surface. Our
Mantle Adjoining Layer tiles are considerably thicker than the CRUST 2.0
tiles. The thinnest Mantle Adjoining Layer tile, at about 11km thick, can be
found beneath the Himalayas in Nepal, while the thickest, at 76km thick, is
found right at the equator. These Mantle Adjoining tiles create a seamless
Earth geo-neutrino source model with no air gaps.

Figure 12: CRUST 2.0 + Mantle Adjoining Layer thickness, in km. The Mantle
Adjoining Layer connects the ellipsoidal 7th CRUST 2.0 layer, the “lower crust
(Mohorovicic discontinuity)”, to the spherical upper mantle. Each 8th layer
tile’s thickness is defined by the gap between the spherical mantle below it
and the floor of the 7th layer tile above it. Zero-thickness tiles are omitted
from this figure.

The eight layers we assume are:

  1. 1.

Ice

  2. 2.

Water

  3. 3.

Soft Sediment

  4. 4.

Hard Sediment

  5. 5.

Upper Crust

  6. 6.

Middle Crust

  7. 7.

Lower Crust (Mohorovicic Discontinuity)

  8. 8.

Mantle Adjoining Layer (not part of CRUST 2.0)

Figure 13: Crust background for a $10^{32}\mathrm{p}^{+}$ detector, saturated
at 100 events per year. Figure 14: Crust+Mantle background for a
$10^{32}\mathrm{p}^{+}$ detector, saturated at 100 events per year. Figure 15:
Random Crust and Mantle geo-neutrino flux vs range from a TREND detector off
the coast of Europe-Atlantic. Refer to Section 6.2 for more information
concerning this detector’s placement and its environment.

Each crust tile has a defined location, thickness (shown in Figure 12), and
density (shown in Figure 16). Elemental abundances for each layer were taken
from Mantovani et al. [59]. Using the CRUST 2.0 densities and volumes, the
Mantovani et al. elemental abundances for the appropriate layer, and the same
constant isotope abundances assumed for the mantle, the total mass of each
isotope of interest can be calculated for each tile. The relationship between
total abundance and neutrino flux is determined from isotope half-life and
multiplicity (the number of neutrinos emitted per decay). The complete
neutrino emissions calculated for each tile in each crust layer is shown in
Figure 17 for 238U and 232Th.

The neutrino energy spectra for the isotopes considered were taken from
Enomoto [78]. The flux and energy spectrum, together, fully define the
neutrino source signal (assumed to be isotropic) for each tile and for each
isotope. In Section 4.1 these source spectra are combined with the detector
cross section, defined in Equation 49 and shown in Figure 29, to create the
observed geo-neutrino spectrum for each tile and isotope.

Figure 16: CRUST 2.0 + Mantle Adjoining Layer density, in g/cm3. Density is
seen to generally increase as depth increases. The top ice layer naturally
shows the lowest density, slightly under g/cm3, while the bottom Mantle
Adjoining Layer density is about 3.4g/cm3. Zero-thickness tiles are omitted
from this figure. Figure 17: CRUST 2.0 + Mantle Adjoining Layer emissions
(238U in upper, 232Th in lower), in $\bar{\nu}_{e}$ per second per tile. This
figure combines density and thickness information from CRUST 2.0 [70] and PREM
[86] with elemental abundances from Mantovani et al. [59], (along with 232Th
half life of $20.3\times 10^{9}$ years and multiplicity of 4$\bar{\nu}_{e}$
per 232Th decay, and 238U half life of $6.45\times 10^{9}$ years and
multiplicity of 6$\bar{\nu}_{e}$ per 238U decay) to compute the total
$\bar{\nu}_{e}$ emissions from a tile per second isotropically. Layers 5 and
6, the Upper Crust and Middle Crust emit the most neutrinos due to their high
elemental abundance [59], despite having lower density and thickness (as seen
in Figures 16 and 12) than many other layers.

The computation burden of modeling the crust geo-neutrino signal, even at this
level of fidelity, is quite significant. For each detector in a scenario,
there are 129,600 (90x180x8) observed energy spectra for each isotope. Each
spectra must be distorted by neutrino oscillations with the appropriate range.
Additionally, the crust tiles near a given detector are too large (about 200km
across at the equator) to be adequately modeled as point sources. Because of
this fact, a “smart” discrete integration process was developed to recursively
subdivide the nearest tiles into progressively smaller sub-tiles until the
contribution of each is below a threshold of 0.1 events per detector-year.
With the large detector sizes in this study, the creation of hundreds of
thousands of additional sub-tiles to replace the nearest high-flux crust tiles
is required. Figure 23 is a Google Earth overlay showing the uppermost layer
of these sub-tiles following smart integration after TREND placement within
the Geospatial Model.

The process described for the crust must be repeated for the mantle in order
to compute the total geo-neutrino signal at a given detector. The mantle is
much easier since spherical symmetry is assumed here and because the “smart”
discrete integration process is rarely needed. The Geospatial Model can
predict the geo-neutrino flux for detector placements at locations around the
Earth. The mantle component will vary neutrino flux only with detector depth,
however the crust component will vary neutrino flux with detector latitude,
longitude and depth. Figure 13 shows the observed crust-only measurement rate
calculated for a $10^{32}\mathrm{p}^{+}$ detector located anywhere in the
world and Figure 14 shows the combined crust+mantle geo-neutrino flux on the
same scale. Figure 15 shows MC geo-neutrino flux contributions separately for
the crust and mantle as a function of range away from a detector placed off
the southern coast of Europe-Atlantic888The Europe-Atlantic scenario is one of
four used for reactor searches in this paper, and Detector #1 shown in Figure
15 is discussed in detail in Section 6.2..

###  3.3 Non-neutrino background

Non-neutrino background can be categorized into two types. The first type is
one where a single complex event mimics both the prompt and delayed signals.
The second type is one where the prompt and delayed signals are of separate
uncorrelated origin, the so called “accidental” background. Accidental
background may be caused by internal residual radioactivity, cosmic rays (only
muons are important at depths greater than a few meters) passing through the
detector, cosmic rays outside the detector which make products such as fast
neutrons which enter the detector, and finally radioactivity outside the
detector, largely gamma rays being the concern (alphas and electrons do not
travel far). The three specific components making up the total non-neutrino
background considered in this study are:

  1. 1.

Accidentals. Accidental prompt and delayed signals are caused by two separate
uncorrelated sources, and occur close enough in time and space to fool a
typical inverse beta decay filter.

  2. 2.

Fast Neutrons. A fast neutron may mimic an inverse beta decay prompt signal
during its random walk process, recoiling off free protons in the liquid
scintillator, and then when finally captured on a dopant (such as Gd) it may
produce a signal indistinguishable from a normal inverse beta decay delayed
signal.

  3. 3.

Cosmogenic 9Li/8He. Cosmogenic isotopes produced by showering muons fall into
this category, and are most problematic due to the long lifetimes of 8He and
9Li isotopes in the detector, which may trick a time- or space-based
coincidence filter designed to detect inverse beta decay events.

Fortunately, inverse beta decay is a very distinctive process that involves a
spatial coincidence on the order of a meter and the temporal coincidence on
the order of 10 $\mu$s between the prompt (positron) and delayed (neutron
capture) signals, followed by a neutron capture providing a known energy
release. Thus, a “delayed coincidence filter” can be used in the data
processing to dramatically reduce the number of background events mistaken for
neutrinos. This and other processing filters reduce the non-neutrino
background, usually at the expense of reducing the effective fiducial volume,
reducing the duty cycle, and/or decreasing the detection efficiency for true
neutrino sources.

In order to estimate the count rate of the non-neutrino background in the
detector, a series of scaling relationships were established based upon
KamLAND and Borexino background measurements. Table 5 shows a summary of these
estimated rates, along with the predicted reactor neutrino and geo-neutrino
background rates. The cosmogenic 9Li/8He and fast neutron count rates and
uncertainties were carefully scaled from KamLAND [72] and the accidental count
rate and uncertainty from Borexino [56]. The geo-neutrino and reactor neutrino
background count rates were computed as described in Section 3 for one of the
scenarios. The background count rates are given for the entire fiducial volume
(no fiducial cuts) and for one year of live-time. The actual background counts
expected over one year must include reductions caused by the muon-related and
cosmogenic-related downtimes, which are very sensitive to detector depth. Note
that the geo-neutrino signal does not disappear entirely with the
$E_{\mathrm{vis}}\geq 2.6$MeV energy cut due to blurring from the energy
measurement resolution.

| KamLAND | Borexino | TREND extrapolated from
---|---|---|---
|  |  |  KamLAND & Borexino
Flat eq. depth | 2,050 m | 3,050 m | 3,500 m
Scintillator | C11.4H21.6 | C9H12 | C16H30
H/m3 | $6.60\;10^{28}$ | $5.30\;10^{28}$ | $6.24\;10^{28}$
C/m3 | $3.35\;10^{28}$ | $3.97\;10^{28}$ | $3.79\;10^{28}$
density | 0.78 | 0.88 | 0.86
Mass (tons) | 912 | 278 | 138,000
Volume (m3) | 1170 | 316 | 160,000
Radius (m) | 6.5 | 4.25 | 23
Cyl. Length (m) | — | — | 96.5
$\mu-$Section (cm2) | $1.3\;10^{6}$ | $0.57\;10^{6}$ | $4.4\;10^{7}$
$\mu-$Flux (cm-2s-1) | $1.6\;10^{-7}$ | $0.3\;10^{-7}$ | $1.4\;10^{-8}$
$\mu-$Energy (MeV) | 219 | 276 | 295
$\mu-$Rate (s-1) | $2.13\;10^{-1}$ | $1.6\;10^{-2}$ | $6.2\;10^{-1}$
$\mu-$DT (200 $\mu$s) | $4\;10^{-5}$ | $0.3\;10^{-5}$ | $12\;10^{-5}$
Co-DT (1500 ms) | $1.0\;10^{-1}$ | $7.5\;10^{-3}$ | $1.5\;10^{-1}$
Exposure (H.y) | $2.44\;10^{32}$ | $6.02\;10^{30}$ | $10^{34}$
$E_{\mathrm{vis}}$ threshold [MeV] | 0.9 | 1 | 0.9 | 1 | 2.6
Accidental Rate | 80.5$\pm$0.1 | 0.080$\pm$0.001 | 277$\pm$3 | 183$\pm$2 | 0.3$\pm$0.003
9Li/8He Rate | 13.6$\pm$1.0 | 0.03$\pm$0.02 | 20$\pm$1 | 20$\pm$1 | 18$\pm$1
Fast n Rate |  $<$9 (4.5$\pm$2.6) |  $<$0.05 (0.025$\pm$0.014) | 31$\pm$18 | 31$\pm$18 | 26$\pm$15
Geo-$\nu$ | 69.7 | 2.5$\pm$0.2 | 2,111 | 2,088 | 18
Known Reactor-$\nu$ | — | — | 499 | 499 | 352
Table 5: Breakdown of the expected background count rates for TREND. The
$\mu-$DT and Co-DT are the estimates of muon and cosmogenic induced downtime
(DT). The flat equivalent depths and muon fluxes are taken from [54].

The accidental rate was modeled as having two components: one that is uniform
within the fiducial volume (dictated by the radiopurity of the LS) and another
that is an exponential function of the (closest) distance from the inner-most
structure. The scale length of the exponential component (0.3 m) and its
magnitude relative to the uniform component at the edge of the structure
(1,400) were inferred from KamLAND measurements of the accidental rate using
various fiducial cuts [71]. The shape of the exponential component was assumed
to hold in both Borexino and TREND, with the absolute scaling being determined
by a fit to the Borexino measurements. For the TREND detector, there is not a
significant fiducial volume cut, so the predicted accidental rate includes
regions of increased rate (per unit volume) near the edges of the cylinder.

The fast neutron rate was modeled as an exponential function of the (closest)
distance from the outer structure that might produce fast neutrons without
having them be tagged to a passing muon. The exponential curve was defined
such that 90% of the rate (per unit volume) comes from within 3 m of the
“edge”. This is somewhat similar to the assumption made by Lasserre et al [2]
that “all” of the fast neutron rate comes uniformly from within 3 m of the
“edge”. This corresponds to a scale length of 1.3 m for the exponential. The
shape of the exponential was assumed to hold in both KamLAND and TREND, with
the absolute scaling being determined by a fit to the KamLAND measurements.
For the TREND detector, this “edge” was considered to be the cylinder itself.

In constructing Table 5, several significant discrepancies with the background
calculations of Lasserre et al. [2] were found. The first discrepancy involved
the calculation of the expected cosmogenic 9Li/8He count rate. Lasserre et al.
scaled the expected cosmogenic 9Li/8He count rate from KamLAND without
accounting for the fact that KamLAND uses a different veto time (2s) than the
$\sim$600ms time proposed for SNIF. This lead Lasserre et al. to significantly
underestimate the expected cosmogenic count rate. This, in turn, contributed
to their conclusion that SNIF can be deployed at depths as shallow as 1500m.
In order to reduced the count rate back to an acceptable level, the veto time
was increased from 600ms in SNIF to 1500ms in TREND.

The second discrepancy involved the calculation of downtime related to
cosmogenic 9Li/8He veto. In their calculations, Lasserre et al. vetoed all
muons with a 3m radius cylinder through the detector along the track; they
never veto the entire fiducial volume for their proposed 600ms veto time. This
did not impose a penalty, in terms of reduced duty cycle, for the large
unsegmented design. Based on KamLAND, approximately 1.5% of muons passing
through the detector were showering[72] for an average muon path length of
approximately 8.7m. Scaling this to the average muon path length computed for
TREND ($\sim$37m) yields a prediction that 6.4% of the muons are showering.
This prediction is highly uncertain and the scaling ignores any dependence on
the average muon energy with depth. To be conservative, we assume 15% of the
muons passing through TREND are showering (regardless of depth) and subject to
veto of the entire fiducial volume. These assumptions caused our estimate of
the cosmogenic downtime to be significantly larger. Figure 18 shows the
computed cosmogenic downtime percentage as a function of depth (for a veto
time of 1500ms). This analysis indicates that a depth of approximately 3500m
is required to reduce the cosmogenic downtime to values less than 15%. A
deeper ocean site increases the stand-off distances involved, making the TREND
concept even more challenging in practice. This study will show, however, that
TREND can still be effective using the estimator presented in this paper.

Figure 18: Downtime due to Cosmogenic 9Li/8He veto as a function of detector
depth in meters of water equivalent and detector volume. The dashed white
lines represent downtimes of 15%, 50%, and 90%. The solid white line indicates
the 160,000 m3 fiducial volume of TREND. The discontinuity below 113 m3 occurs
because the detector becomes too small to use the r=3m veto cylinder; the
entire volume is vetoed even for well-reconstructed muon tracks. Note: this
figure assumes spherical detector geometry, that 0.17% of muons are showering
per meter of average muon path length (regardless of depth), and that the
entire volume is vetoed for showering muons. The calculations for TREND use
the actual cylindrical geometry and a conservative estimate of 15% showering
muons; this yields a larger downtime than implied by this figure (e.g. 15%
downtime at a depth of $\sim$3,500m).

The third discrepancy involved the calculation of the expected accidental
count rate. Lasserre et al. scaled the accidental count rate from Borexino
without accounting for the fact that Borexino uses a different delayed
coincidence filter ($\tau$=1280 $\mu$s, $\Delta$R=1 m) than the filter they
proposed ($\tau$=30 $\mu$s, $\Delta$V=1 m3). This lead Lasserre et al. to
significantly overestimate the expected accidental count rate. While this is
good news for the technical feasibility of TREND, Lasserre et al. correctly
point out that a concept as futuristic as TREND can safely expect the higher
level of Borexino’s radiopurity due to technological improvements over the
next 30 years. Based on detailed modeling of the TREND detector design, a
modified delayed coincidence filter is proposed here ($\tau$=200 $\mu$s,
$\Delta$R=2 m) to balance the competing desires to minimize the expected
accidental count rate while preserving a relatively high neutrino detection
efficiency.

Figure 19 shows the non-neutrino background as a function depth. It is
interesting to note that the muon-related non-neutrino background no longer
dominates for depths greater than 2000m. What appears to be preventing TREND
from operating at a shallow depth (2000m - 2500m) is the high cosmogenic
downtimes (100% - 70%) incurred primarily by the vetoing of showering muons in
such a large unsegmented detector.

Figure 19:  Non-neutrino background rate as a function of the detector depth
in meters of water equivalent for an exposure of $10^{34}\mathrm{p}^{+}$-year.
The limitation line for accidentals depends upon the detailed detector design,
as it will be largely due to radiation from outside the detector or the
detector structure and components. This limitation may be considerably lower
in practice for TREND, and in any event is largely confined to the peripheral
region of the detector. The same is true of the entering fast neutron
background.

There is a well-known seasonal variation in the muon flux on the surface and
at depth due, primarily, to seasonal variations in the atmospheric temperature
profile [97]. Ambrosio et al. found a 2% seasonal variation with the MACRO
detector located at the San Grasso underground laboratory in Italy.
Furthermore, the relative magnitude of this variation is expected to increase
with the average muon energy and, thus, with increasing depth. For simplicity,
this paper assumes that the (Poisson mean) muon flux, and its contribution to
the non-neutrino background, is constant in time regardless of the detector
depth. In actual operation, the expected temporal variation in muon flux
during an experiment can be accounted for if it is deemed necessary.

For simplicity an isotropic angular distribution is assumed for each non-
neutrino background. Thus, unlike reactor and geo-neutrino background, the
non-neutrino background have the same energy spectrum in all directions seen
by a detector. Even if this assumption is not strictly true, the non-neutrino
background very likely have a slowly and smoothly varying angular
distribution. The nearly isotropic directional measurement resolution of the
detector further reduces the importance of modeling the angular distribution
of non-neutrino background with high fidelity.

The non-neutrino spectra were all taken from empirical observations. For fast
neutron background and cosmogenic background the energy spectra defined by
Mention et al. [55] were used. For accidentals the energy spectrum was
obtained by fitting an exponential distribution (see Equation 13) to KamLAND
data [79]. Figure 20 shows the fit. Being based on observations from other
detectors, these spectra likely have their own energy measurement resolution
and energy-dependent efficiencies applied to them. For this study, the
observed spectra were simply assumed to represent the source spectra already
smeared by our modeled detector energy resolution. Thus the non-neutrino
spectra are somewhat notional for this detector. This approach is not perfect,
but it was deemed more than sufficient for this study. The three non-neutrino
energy pdfs can be seen in normalized form in Figure 21 alongside four
neutrino source categories: known reactors, unknown reactor, crust and mantle.

$p(E_{\mathrm{vis}})=\left(2.25\times 10^{5}\right)e^{-6.185E_{\mathrm{vis}}}$
(13)

Figure 20: Exponential least squares fit to experimental KamLAND data from
[79]. Fit was performed using KamLAND data from $E_{\mathrm{vis}}=0.9$MeV to
$E_{\mathrm{vis}}=1.7$MeV, then applied over the wider energy spectrum seen in
this figure. Figure 21: Smeared energy spectrum of each measurement source
type in $\theta_{k}$ as seen by TREND detector #3 in our Europe-Mediterranean
unknown-reactor scenario (refer to Section 6.4 for more information concerning
this detector’s placement and its environment). The known-reactor spectrum is
a mixture distribution ($Z=w_{n}\sum\limits_{n}H(\theta_{n},d)+v(d)$) of all
the individual known-reactor spectra from around the world, which tend to
“wash out” the effect of oscillations from each source when summed together.
See Abe et. al. [51] for more information concerning cosmogenic background
energy spectra.

###  3.4 Uncertainty model

The estimator presented in Section 4 naturally incorporates the expected count
rate for background sources, similar to the approach taken by Lasserre et al.
[2]. This work additionally allows the mean count rate for each background
source, the mean free proton count of each detector, and the neutrino
oscillation parameters to have realistic Gaussian-distributed systematic
uncertainty. The inclusion of these real-world uncertainties addresses many
questions raised by simpler assumptions and increases MC fidelity. While this
widens the uncertainty contours about any maximum a posteriori estimate, it
acknowledges the real uncertainties involved in accurate parameter estimation
and leads to a more accurate a posteriori uncertainty interpretation about the
maximum a posteriori estimates.

The a priori systematic uncertainties used by the Geospatial Model are shown
in Table 6. The uncertainty model assumes that the energy spectra and angular
distributions for the various flux sources have no uncertainty to their
intrinsic shapes. All of the uncertainty for a particular point source, or for
a class of point sources, is assumed to be located in the mean flux rate
itself. For flux originating from discrete sources (crust with CRUST 2.0
tiles, known reactors), there may be a systematic uncertainty that affects the
entire class ($\sigma_{\mathrm{syst_{class}}}$) as well as a systematic
uncertainty which is independent for each source in the class
($\sigma_{\mathrm{syst_{source}}}$). For continuous background sources (non-
neutrinos), where $\sigma_{\mathrm{syst_{source}}}$ = 0%, the “number of
sources” column indicates the number of discrete sample points used to
approximate the continuous background in the Geospatial Model. It is important
to note that the background uncertainties described here are distinct from the
statistical uncertainty that naturally arises from Poisson noise.

The assumption that the uncertainty can be restricted to the mean flux rate
seems quite reasonable for several reasons. We assume that careful modeling
and extensive calibration runs would be performed on TREND during its design
and construction to quantify the spectral character of non-neutrino background
surviving the various processing filters. For known-reactor flux and geo-
neutrino flux, the spectral shape of the intrinsic signal is relatively well
understood (to within a few percent). The noted sterile neutrino uncertainty
in Table 6 refers to the present controversy referred to as the Reactor
Neutrino Anomaly [108]. Reanalysis of old data and calculation presently point
to a deficit in near reactor counting rates as compared to calculated reactor
neutrino fluxes, of the order of $6.3+/-2.7\%$. This plus other hints may
point towards a small level oscillation coupling from normal neutrinos into
sterile neutrinos. The mass difference determined range is implied to be a few
meters, such that the oscillations are in equilibrium for almost all near-
reactor experiments. Detectors located further away, as is the case for all
situation considered in this paper, the result of the putative sterile
neutrinos would be a few percent reduction in the expected flux, and some
rather small shift in oscillation parameters. At present there are also some
peculiarities in the spectra of the experiments (Daya Bay, Double Chooz and
RENO) measuring $\theta_{13}$ from about 2 km away from reactors and reported
at conferences. This may be pointing to problems with the calculations of
reactor borne neutrino spectra. Of course it is reasonable to expect that this
anomaly will be resolved by the time the long range detectors discussed herein
will be built, and will require in any event only minor adjustment of mixing
parameters, but will change no conclusions about our results.

Finally, small uncertainties in the spectral shapes are smoothed out by the
energy measurement noise and the poor directional resolution of TREND999“Poor”
directional resolution of TREND defined as $\sim$0.045 SNR per Section 2.4.3..

Source | Uncertainty | Number of Sources | Comment
---|---|---|---
| $\pm 1\sigma_{\mathrm{syst_{class}}}\pm 1\sigma_{\mathrm{syst_{source}}}$ |  |
Known reactor flux | $\pm 2.0\%\pm 3.4\%$ | 471 reactor cores | does not include sterile $\nu$ uncertainty
Crust flux | $\pm 20.0\%\pm 8.0\%$ |  $\sim 129600$ tiles | 100% correlated across 238U and 232Th
Mantle flux | $\pm 50.0\%\pm 0\%$ |  $\sim 50396$ tiles | 100% correlated across 238U and 232Th
Cosmogenic 9Li/8He | $\pm 3.3\%\pm 0\%$ | 4584/detector | KamLAND [50] extrapolation
Fast neutron flux | $\pm 10.0\%\pm 0\%$ | 4584/detector | KamLAND [50] (capped at 10.0%)
Accidental | $\pm 1.3\%\pm 0\%$ | 4584/detector | Borexino [72, 56] isopurity extrapolation
Detector proton count | $\pm 0.0\%\pm 2.4\%$ | 1/detector | i.e. fiducial volume uncertainty
$\Delta m^{2}_{12}$ | ${}^{+2.4\%}_{-2.8\%}\pm 0\%$ | $\sim 180467$ | per Fogli et al. [63]
$\Delta m^{2}_{13}$ | ${}^{+3.4\%}_{-2.6\%}\pm 0\%$ | $\sim 180467$ | per Fogli et al. [63]
$\sin^{2}\theta_{12}$ | ${}^{+5.2\%}_{-5.7\%}\pm 0\%$ | $\sim 180467$ | per Fogli et al. [63]
$\sin^{2}\theta_{13}$ | ${}^{+10.0\%}_{-9.8\%}\pm 0\%$ | $\sim 180467$ | per Fogli et al. [63]
Table 6: Breakdown of the a priori systematic uncertainties assumed by the
Geospatial Model. These uncertainties are assumed zero-mean normally
distributed, and are each randomly sampled once per Geospatial Model MC.

Significant effort was spent on creating a realistic uncertainty model,
populated with reasonable values. The uncertainty model also reflects the long
distance, multi-detector nature of the scenarios in this study by allowing for
the fact that the detected signal is not dominated by one detector or by one
source. The systematic uncertainties for the non-neutrino background were
generally taken from KamLAND and Borexino measurements. The one exception was
fast neutrons, whose highly uncertain rate was based upon a measured upper
limit. It was assumed that, for a futuristic detector such as TREND,
measurement methods would exist to quantify this to 10%. The 3.4% systematic
uncertainty of each known reactor core corresponds to the combined reactor-
related system uncertainty from S. Abe et al.[50]. The 2.4% systematic
uncertainty of each detector corresponds to the combined detector-related
system uncertainty from S. Abe et al.[50]. Finally, oscillation parameter
uncertainty was defined per Fogli et al. [63] and applied to all neutrino
sources present in the Geospatial Model, including 471 reactor cores, all
crust tiles and the mantle.

##  4 Estimation theory

A Bayesian maximum a posteriori (MAP) estimator is presented here, suitable
for a wide variety of neutrino experiments. The estimator utilizes all the
available measurements and exploits all available prior information such as
the measurement noise (detector hardware), uncertainty in the background
signal, and constraints on the parameters to be estimated. The
optimal101010Optimal is defined here as using all readily available
measurements (count rate, energy, angle, no energy cut) and marginalizing over
(or co-estimating) all parameters with significant uncertainty. estimator is
presented, as well as two tractable suboptimal approximations appropriate for
real world implementation with finite computing resources. A few example
applications include the refinement of neutrino oscillation parameters, the
characterization of the geo-neutrino background, or placing an upper limit on
the power of a geo-reactor.

In this paper we test the estimator on two separate applications. First, it is
used to refine $\bar{\nu}_{e}$ oscillation parameters and to reduce systematic
uncertainty related to our a priori knowledge of backgrounds (e.g. geo-
neutrinos and reactor neutrinos). Second, it is used to geolocate an unknown
reactor. Both of these applications assume that four TREND detectors are
operating simultaneously and that they combine their detections for the
estimation process.

###  4.1 Optimal MAP estimator

In the Bayesian approach, one starts with the a priori probability density
function of the parameter space, $p(\theta)$, from which one can obtain its a
posteriori pdf:

$p(\theta|Z)=\frac{p(Z|\theta)p(\theta)}{p(Z)}=\frac{1}{c}p(Z|\theta)p(\theta)$
(14)

where $c=p(Z)$ acts as a normalization constant111111In the Bayesian approach,
$p(Z)$ is a constant scalar value (sometimes referred to as the ‘marginal
likelihood’) equal to an integral across the entire measurement space $Z$.It
is deterministic in nature since it is not evaluated at any specific
measurement. Since it is not a function of the parameter space $\theta$ and
does not vary as a function of the random measurements, it follows that the a
posteriori probability may be expressed as proportional to the numerator only
of Bayes’ equation., which, since it does not depend on $\theta$ can be
considered irrelevant for maximization [82]. $Z$ is a stacked vector of
$n_{z}$ measurements. $p(\theta)$ is the Bayesian prior which describes our a
priori knowledge of the parameter space $\theta$. $p(Z|\theta)$ is the
conditional probability of $Z$ given $\theta$, commonly referred to as the
likelihood function or, more simply, as the likelihood.

The measurement vector $Z$, corrupted by a detector-dependent zero-mean white
noise $v(d)$, is a function of the parameter vector $\theta$ and the detector
$d$. In a linear system, $H$ would be a Jacobian of partial derivatives
expressing the relationship between $Z$ and $\theta$: $Z=H\theta+v$. In our
nonlinear system, $H$ is a function which translates the parameter vector
$\theta$ into the smeared measurement space $Z$:

$Z=H(\theta,d)+v(d)$ (15)

The parameter vector $\theta$ is composed of a stacked vector of $n_{u}$
unknown sources (e.g. reactors) to be estimated ($\theta_{u}$), $n_{k}$ known
sources ($\theta_{k}$), and various physical constants contained in
$\theta_{c}$. $\theta_{k}$ contains the known background sources: geo-
neutrinos, known-reactor neutrinos, accidentals, fast neutrons, and
cosmogenics. $\theta_{c}$ contains all the physical constants employed by the
estimator, such as the solar mixing angle and the constants contained within
PREM used for geo-neutrino flux calculations.

$\displaystyle\theta$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\theta_{u}\\\ \theta_{k}\\\
\theta_{c}\end{array}\right]$ (19) $\displaystyle\theta_{u}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{c}\theta_{u_{1}}\\\ \vdots\\\
\theta_{u_{n_{u}}}\end{array}\right]$ (23) $\displaystyle\theta_{k}$
$\displaystyle=$ $\displaystyle\left[\begin{array}[]{c}\theta_{k_{1}}\\\
\vdots\\\ \theta_{k_{n_{k}}}\end{array}\right]$ (27)

To simplify the problem, we assume each unknown reactor in $\theta_{u}$ can be
accurately modeled as a point source which emits neutrinos isotropically, and
that all the unknown reactors are constrained to lie on the Earth’s surface.
Following these assumptions each unknown reactor in $\theta_{u}$ can now be
defined by just four elements: geodetic latitude $\phi$, geodetic longitude
$\lambda$, thermal power output $p$ (i.e. flux), and zero-range energy
spectrum $s$121212Note that the zero-range spectrum $s$ is a function of
energy, $s(E)$, but for simplicity it is assumed that all sources have one of
a finite number of well-defined energy spectra. Thus $s$ is treated as a
discrete parameter.. For $n_{u}$ unknown reactors:

$\theta_{u}=\left[\begin{array}[]{cccc}\phi_{1}&\lambda_{1}&p_{1}&s_{1}\\\
\vdots&\vdots&\vdots&\vdots\\\
\phi_{n_{u}}&\lambda_{n_{u}}&p_{n_{u}}&s_{n_{u}}\end{array}\right]$ (28)

We further simplify the problem by assuming that only one unknown source
exists in the vicinity of our detector. Thus $n_{u}=1$ and $\theta_{u}$
reduces to one row. For this 1x4 row vector $\theta_{u}$, we apply a four-
dimensional uniform (but bounded) a priori probability $U(\theta_{u})$. This
effectively constrains the unknown source to exist within a certain latitude-
longitude-flux-spectrum domain with equal a priori probability at every point.
$p(\theta_{u})$ need not be uniform, however, or even linear. Any informative
a priori probability may be used for $p(\theta_{u})$, but it must be properly
normalized to reflect a 100% confidence that the unknown source exists
somewhere within the full bounds of the four-dimensional parameter space under
consideration. Note that behavior of the estimator when the unknown source is
not, in reality, contained within $p(\theta_{u})$ is examined later on in
Section 6.1.

$\theta_{k}$ represents the known background sources expected at the detector
location, independent of the unknown sources in $\theta_{u}$. 232Th and 238U
isotope geo-neutrino sources in the mantle and crust are considered here,
though potassium 40K isotopes are ignored as their neutrino energies are below
the inverse beta decay energy measurement threshold of about
$E_{\bar{\nu}_{e}}\geq 1.8$MeV. The combined flux from all neutrino sources
for a $10^{32}p^{+}$ detector are shown in Figure 22.

We have eight known categories of background in $\theta_{k}$:

$\theta_{k}=\left[\begin{array}[]{c}\theta_{\mathrm{IAEA}}\\\
\theta_{\mathrm{crust_{{}^{238}U}}}\\\ \theta_{\mathrm{crust_{{}^{232}Th}}}\\\
\theta_{\mathrm{mantle_{{}^{238}U}}}\\\
\theta_{\mathrm{mantle_{{}^{232}Th}}}\\\ \theta_{\mathrm{accidental}}\\\
\theta_{\mathrm{fast\ }n^{0}}\\\
\theta_{\mathrm{cosmogenic}}\end{array}\right]$ (29)

Figure 22: Background $\bar{\nu}_{e}$ flux expected at varying locations
worldwide from the $\bar{\nu}_{e}$ sources in $\theta_{k}$. Mantle flux
computed from PREM [86] and Mantovani et al. [73]. Crust data computed from
CRUST 2.0 model and Mantovani et al. [73]. Reactor flux computed from publicly
available IAEA data [85].

Nearly all the neutrino source locations in $\theta_{k}$ are most naturally
expressed in geodetic, or ECEF reference frames, though the non-neutrino
background may be more naturally expressed in a local spherical coordinate
system with respect to the detector center-point. Some coordinate system
transformations may need to be defined if sources are defined with respect to
different origins, or expressed in different reference frames131313Note that
these low level coordinate system transformations are implicitly assumed in
all the equations in this paper..

Some sources in $\theta_{k}$ are assumed to be point sources, such as the
known reactors in $\theta_{\mathrm{IAEA}}$. Others are clearly volume sources,
such as the mantle and crust 238U and 232Th isotopes. These volume sources can
be treated as integrations of point sources, where each point in the discrete
integral would occupy a row in $\theta_{k}$. Discrete volume integral step
sizes must be carefully chosen to ensure accurate approximation to a
continuous volume integral, especially when dealing with volumes in the
immediate vicinity of the detector such as the Earth’s oceans and crust.

To avoid excessive computational burden caused by small integral step-sizes
over large volumes, a recursive “smart” integration technique may be
implemented, varying local integral step-size efficiently while retaining high
model fidelity. Large volume sources are recursively broken into groups of
smaller volume sources until a predetermined flux threshold is no longer
exceeded by each resultant discrete integral point. In our tests we found a
threshold of 0.1 events/detector/year per discrete point source to offer an
excellent compromise between discrete integral accuracy and computational
burden. Figure 23 shows one such variable step-size crust integration for a
detector submersed in the waters of the North Atlantic.

Figure 23: Example “smart” spherical integration of the Earth’s upper crust, a
volume source. Gray squares are original CRUST 2.0 Upper Crust tiles. Each
grey tile is $2^{\circ}\times 2^{\circ}$. Red tiles are $1^{\circ}\times
1^{\circ}$, green tiles are $0.5^{\circ}\times 0.5^{\circ}$ and so on, up to a
30th generation tile, at $\left(4\times
0.5^{-30}\right)^{\circ}\times\left(4\times 0.5^{-30}\right)^{\circ}$. Note
that smart discrete integration is performed in three dimensions
simultaneously, though the height axis is not shown in this figure.

It is natural for the expected flux of known sources in $\theta_{k}$ to
contain some degree of systematic uncertainty. An example of this would be the
large systematic uncertainty in the expected flux of geo-neutrinos produced
from the mantle. We assume here that any systematic uncertainty in
$\theta_{k}$ can be fully (and accurately) defined by an a priori probability,
$p(\theta_{k})$. Uncertainty in many elements of $\theta_{c}$, such as the
solar mixing angle uncertainty, is also expected as well. We, therefore, also
assume that a well-defined prior probability $p(\theta_{c})$ can describe
these uncertainties.

There are two mutually exclusive options141414Approximation is a third option,
used in our suboptimal approximations in Section 4.2. for optimally dealing
with such parameter space uncertainty: marginalization or estimation. In
marginalization, the parameter is treated as a nuisance parameter, and a
weighted integration through the parameter dimension is performed. This
results in a marginal probability which is no longer conditional on the
uncertain parameter:

$p(x)=\int\limits_{-\infty}^{\infty}p(x,i)p(i)di$ (30)

where $i$ is the uncertain parameter, $p(i)$ is the Bayesian prior of $i$,
$p(x,i)$ is the conditional likelihood and $p(x)$ is the marginal likelihood.
Note that $p(x)$ is independent of $i$ following marginalization. This method
can be extrapolated to correlated multidimensional uncertainties through the
use of higher dimensional integrals, albeit at a much high computational cost.
It may be more easily applied to independent (zero-covariance)
multidimensional uncertainties via the product of several single dimensional
integrals.

The second option when dealing with uncertain nuisance parameters in
$\theta_{k}$ or $\theta_{c}$ is to co-estimate more accurate values for them
alongside the original parameters of interest, $\theta_{u}$ in this case.
Assuming at least partial observability of the uncertain nuisance parameters
from the measurements, $Z$, this option may deliver superior results to
marginalization.

There is no correct answer for deciding when to marginalize or optimize, the
choice depends largely on the problem at hand and the computational resources
available. Marginalization usually expends more computational resources as it
requires more objective function evaluations. In general, information is lost
during marginalization as a higher dimensional a priori pdf collapses into a
lower dimensional a priori pdf. This results in flatter a posteriori pdf with
smaller gradients and larger confidence bounds about the maximum a posteriori
estimate. For these reasons this method should only be applied when there is
quite high uncertainty in the “known” parameter value, and poor prospects for
estimating a better value (i.e. low parameter observability). Marginalization
may also be appealing when faced with higher degrees of freedom, nonlinear or
multimodal distributions, or deteriorated a posteriori parameter space
topology inhibiting reliable optimizer performance.

So, recognizing the need to marginalize some of the uncertainty dimensions
while optimizing others, we update Bayes’ equation with the requisite
multidimensional integral through the uncertain nuisance dimensions we have
chosen to marginalize:

$p(\theta|Z)=\frac{1}{c}\int\limits_{-\infty}^{\infty}\cdots\int\limits_{-\infty}^{\infty}p(Z|\theta)p(\theta)d\theta_{1}\cdots\theta_{n}$
(31)

Now, the estimation problem can be defined as the maximization of Bayes’
posterior probability density $p(\theta|Z)$ over the parameters in $\theta$
that remain following marginalization of the nuisance parameters:

$\mathrm{arg}_{\theta}\mathrm{max}~{}p(\theta|Z)$ (32)

In other words, we want to create a maximum a posteriori (MAP) estimator to
find a value for $\theta$ which maximizes $p(\theta|Z)$ conditional on
$p(\theta)$, post-marginalization of the nuisance elements of
$\theta$.151515Note that while this paper primarily covers estimation of
oscillation parameters and unknown nuclear reactors that may exist in
$\theta_{u}$, the estimator is not limited to these applications. Any
uncertain element in $\theta$, including the neutrino physics constants in
$\theta_{c}$, such as possible sterile neutrino mixing angles and mass
differences, may be estimated by allowing an optimizer to travel along those
additional dimensions in $\theta_{c}$ rather than marginalizing them.

Prior to optimization (maximization of the likelihood function, or
alternatively, minimization of the negative log-likelihood function), each
element of the Equation 32 must be defined. Bayes’ normalization constant,
$p(Z)$, or $c$ is addressed first. $c$ is intended to normalize Bayes’ a
posteriori pdf into a proper161616A “proper” pdf integrates to unity, (i.e.
cdf=1). pdf, however since $c$ is irrelevant during optimization, it does not
require definition and is ignored from here onward.

$p(\theta_{u})$ has been defined to be a 4D uniform distribution,
$p(\theta_{u})=U(\theta_{u})$, across our latitude-longitude-power-spectrum
parameter space. A priori pdfs for our other parameters have been
assumed171717The exact form that these uncertainties take is not important for
description of the estimator, though every effort should be make to accurately
represent them when these equations are implemented in an experiment. for
$p(\theta_{k})$ and $p(\theta_{c})$ as well. The only remaining term in Bayes’
equation we have not yet explicitly defined is the likelihood function,
$p(Z|\theta)$.

We begin by examining the measurement vector $Z$. $Z$ is composed of energy
measurements ($Z_{E}$) in the first column and direction vector measurements
($Z_{V}$) in the last three columns.

$Z=\left[\begin{array}[]{cc}Z_{E}&Z_{V}\end{array}\right],\;\;\;\;\;Z_{E}=\left[\begin{array}[]{c}E_{1}\\\
\vdots\\\
E_{n_{Z}}\end{array}\right],\;\;\;\;\;Z_{V}=\left[\begin{array}[]{ccc}V_{1}^{x}&V_{1}^{y}&V_{1}^{z}\\\
\vdots&\vdots&\vdots\\\
V_{n_{Z}}^{x}&V_{n_{Z}}^{y}&V_{n_{Z}}^{z}\end{array}\right]$ (33)

$Z_{V}$ is simply the Cartesian vector which connects the estimated IBD
delayed center-of-energy (CE) to the prompt CE, following a neutrino IBD event
within a liquid scintillator detector. In our Detector Model (Section 2),
prompt and delayed CEs are separately estimated by two variants of the CHOOZ
[3] maximum likelihood estimator (MLE) for each event181818$Z_{V}$ is
expressed here in the ECEF reference frame, though $Z_{V}$ is naturally first
recorded in the local detector reference frame before rotation into the ECEF
frame. This is a natural prerequisite to combining direction vector
measurements from multiple detectors with different orientations in the MAP
estimator. Units of $Z_{V}$ assumed in this paper are mm..

$Z_{V}=\widehat{CE}_{\mathrm{delayed}}-\widehat{CE}_{\mathrm{prompt}}=\left[\begin{array}[]{ccc}\hat{x}^{CE}_{1}&\hat{y}^{CE}_{1}&\hat{z}^{CE}_{1}\\\
\vdots&\vdots&\vdots\\\
\hat{x}^{CE}_{n_{Z}}&\hat{y}^{CE}_{n_{Z}}&\hat{z}^{CE}_{n_{Z}}\end{array}\right]_{\mathrm{delayed}}-\left[\begin{array}[]{ccc}\hat{x}^{CE}_{1}&\hat{y}^{CE}_{1}&\hat{z}^{CE}_{1}\\\
\vdots&\vdots&\vdots\\\
\hat{x}^{CE}_{n_{Z}}&\hat{y}^{CE}_{n_{Z}}&\hat{z}^{CE}_{n_{Z}}\end{array}\right]_{\mathrm{prompt}}$
(34)

The two CHOOZ MLEs co-estimate the visible energy as well as position of both
CEs. The prompt visible-energy estimate, $\widehat{E}_{\mathrm{prompt}}$ is
highly correlated with the actual neutrino energy and is converted to an
estimate of neutrino energy via Equation 35. The MAP estimator has been
designed to exploit these neutrino energy estimates because the detected
neutrino energy spectrum encodes significant information about the range and
spectrum of the contributing sources. Thus we retain the corrected CHOOZ MLE
prompt energy estimate, $Z_{E}$ in units of $\mathrm{MeV}$, as an estimate of
the neutrino energy. Note that the delayed visible-energy estimate is used in
the candidate selection cut, but not directly in the MAP estimator.

$Z_{E}=\left[\begin{array}[]{c}\widehat{E}_{1}\\\ \vdots\\\
\widehat{E}_{n_{Z}}\end{array}\right]_{\mathrm{prompt}}+(E_{n^{o}})-(E_{p^{+}})-(E_{e^{+}})\approx\left[\begin{array}[]{c}\widehat{E}_{1}\\\
\vdots\\\
\widehat{E}_{n_{Z}}\end{array}\right]_{\mathrm{prompt}}+0.782\mathrm{MeV}$
(35)

The measurements $Z$ natively reside within the
$\left[\begin{array}[]{cccc}E&V_{x}&V_{y}&V_{z}\end{array}\right]$ measurement
space, and each individual measurement may be thought of as a random sample of
this continuous 4D measurement space (effectively a 4D pdf). The MAP estimator
exploits only 3 degrees of freedom (energy, azimuth of direction, and
elevation of direction) from each measurement; it does not utilize the CE-to-
CE distance. One may thus normalize the direction vectors and convert them to
spherical elevation and azimuth angles, “throwing away” the unneeded CE-to-CE
distance information in the process (assuming this event has already been
passed candidate selection cut).

Examples of such 3D measurement spaces are shown in Figures 24 to 27 for each
source category (unknown reactor, known reactors, crust, mantle and non-
neutrinos) as well as for the cumulative signal. Figures 24 and 26 compare the
full 3D measurement space at two different (highly exaggerated) levels of
detector angular resolution. Figures 25 and 27 show representations of Figures
24 and 26 respectively where the energy dimensional has been fully integrated
(or marginalized), resulting in a marginal 2D elevation-azimuth measurement
space. The unique directional characteristics of each background source can be
seen in these figures, including the isotropic nature of the non-neutrino
background signal we have chosen to assume.

Figure 24: Europe-Mediterranean scenario TREND detector #1 measurement space
(showing 2.5 angular SNR), a mixture distribution
($Z=w_{n}\sum\limits_{n}H(\theta_{n},d)+v(d)$) of measurement spaces defined
by $\theta_{u}$ and $\theta_{k}$. 90% probability confidence isosurface is
shown in silver; red line represents true direction of unknown reactor. More
information concerning this detector within the Europe-Mediterranean unknown-
reactor scenario can be found in Section 6.4. Figure 25: Europe-Mediterranean
scenario TREND detector #1 2D elevation-azimuth marginal measurement space
following integration (across energy) of Figure 24 3D measurement spaces.
Continent outlines are shown in white. More information concerning this
detector within the Europe-Mediterranean unknown-reactor scenario can be found
in Section 6.4. Figure 26: Europe-Mediterranean scenario TREND detector #1
measurement space (showing 8.4 angular SNR), a mixture distribution
($Z=w_{n}\sum\limits_{n}H(\theta_{n},d)+v(d)$) of measurement spaces defined
by $\theta_{u}$ and $\theta_{k}$. 90% probability confidence isosurface is
shown in silver; red line represents true direction of unknown reactor. More
information concerning this detector within the Europe-Mediterranean unknown-
reactor scenario can be found in Section 6.4. Figure 27: Europe-Mediterranean
scenario TREND detector #1 2D elevation-azimuth marginal measurement space
following integration (across energy) of Figure 26 3D measurement spaces.
Continent outlines are shown in white. More information concerning this
detector within the Europe-Mediterranean unknown-reactor scenario can be found
in Section 6.4.

Before defining $p(Z|\theta)$ we address a somewhat unique aspect of this
problem: information is contained not just within the measured energies and
direction vectors, but also in the number of measurements: $n_{Z}$, defined by
the row count of measurement vector $Z$. This aspect of the problem demands
application of Poisson statistics within the likelihood $p(Z|\theta)$. Thus a
new measurement count likelihood $p(n_{Z}|\theta)$ is combined with the energy
and directional likelihood $p(Z_{E},Z_{V}|\theta)$ to get a more explicit
definition of the overall likelihood $p(Z|\theta)$:

$p(Z|\theta)=p(n_{Z}|\theta)p\left(Z_{E},Z_{V}|\theta\right)$ (36)

where the measurement count likelihood $p(n_{Z}|\theta$) is simply a Poisson
pdf ($\mathcal{P}$):

$p\left(n_{Z}|\theta\right)=\mathcal{P}\left(n_{Z},\bar{n}\right)$ (37)
$\bar{n}=\sum\limits_{j=1}^{n_{u}}\bar{n}\left(\theta_{u}^{j}\right)+\sum\limits_{l=1}^{n_{k}}\bar{n}\left(\theta_{k}^{l}\right)$
(38)

The quantity $n_{Z}$, a scalar integer, is the observed number of candidate
measurements in the detector, $\bar{n}(\theta_{k}^{l})$ is the expected number
of candidate measurements at the detector from known source $\theta_{k}^{l}$,
$n_{k}$ is the number of known sources, $\bar{n}(\theta_{u}^{j})$ is the
expected number of candidate measurements from unknown source
$\theta_{u}^{j}$, and $n_{u}$ is the number of unknown sources.

The second component of the likelihood, $p(Z_{E},Z_{V}|\theta)$, is the
conditional likelihood of all the energy and direction measurements given
$\theta$. Considering all the measurements to be independent of each other,
the likelihood of all $n_{Z}$ measurements equals the product of the
likelihoods of each of the measurements:

$p\left(Z_{E},Z_{V}|\theta\right)=\prod\limits_{i=1}^{n_{Z}}p\left(Z_{E}^{i},Z_{V}^{i}|\theta\right)$
(39)

The energy and direction likelihood for one single measurement $i$,
$p(Z_{E}^{i},Z_{V}^{i}|\theta)$, corrupted by true energy-dependent white
measurement noise $v_{E}(E^{i},d)$, but assumed (in the estimator) corrupted
by measured energy-dependent white measurement noise $\sigma(Z_{E}^{i},d)$,
equals the sum of a convex191919A convex combination is a weighted sum of a
finite number of points where all coefficients are non-negative and sum up to
unity. combination of likelihoods representing a finite set of possible
sources, commonly known as a mixture distribution. Stated in words, the
mixture distribution simply states that each measurement must come from some
source included in the model (be it a known source or an unknown source that
we are trying to estimate).

$p\left(Z_{E}^{i},Z_{V}^{i}|\theta\right)=\sum\limits_{j=1}^{n_{u}}w_{j}p\left(Z_{E}^{i},Z_{V}^{i}|\theta_{u}^{j},\theta_{c}\right)+\sum\limits_{l=1}^{n_{k}}w_{l}p\left(Z_{E}^{i},Z_{V}^{i}|\theta_{k}^{l},\theta_{c}\right)$
(40)

The mixture weight $w_{j}$ for unknown source $j$, $\theta_{u}^{j}$, is
defined as the number of expected measurements from source $\theta_{u}^{j}$
divided by the total number of measurements expected from all sources in
$\theta$, known and unknown:

$w_{j}=\frac{\bar{n}\left(\theta_{u}^{j}\right)}{\sum\limits_{j}\bar{n}\left(\theta_{u}^{j}\right)+\sum\limits_{l}\bar{n}\left(\theta_{k}^{l}\right)}$
(41)

The mixture weight $w_{l}$ for the known source $l$, $\theta_{k}^{l}$, is
correspondingly defined:

$w_{l}=\frac{\bar{n}\left(\theta_{k}^{l}\right)}{\sum\limits_{j}\bar{n}\left(\theta_{u}^{j}\right)+\sum\limits_{l}\bar{n}\left(\theta_{k}^{l}\right)}$
(42)

where $\bar{n}(\theta_{u}^{j})$ is just the cumulative distribution function
(cdf) of the numerator inside the integral in Equation 44,
$\mathrm{cdf}_{d}(\theta_{u}^{j})$, defined in Equation 45 below for detector
d. Thus $\bar{n}(\theta_{u}^{j})=\mathrm{cdf}_{d}(\theta_{u}^{j})$ and
$\bar{n}(\theta_{k}^{l})=\mathrm{cdf}_{d}(\theta_{k}^{l})$.

The likelihood of a single measurement $i$ given a single unknown source
$\theta_{u}^{j}$202020Note that the process is nearly identical when defining
the likelihood of a single measurement $i$ given a single known source
$\theta_{k}^{l}$ can now be defined. The likelihood of energy and direction
measurement $i$ given unknown source $\theta_{u}^{j}$ and constants
$\theta_{c},$ prior to weighting by $w_{j}$ can finally be split into
independent energy and direction likelihoods:

$p\left(Z_{E}^{i},Z_{V}^{i}|\theta_{u}^{j},\theta_{c}\right)=p\left(Z_{E}^{i}|\theta_{u}^{j},\theta_{c}\right)p\left(Z_{V}^{i}|\theta_{u}^{j},\theta_{c}\right)$
(43)

allowing each of these independent likelihoods to be defined in turn. First
the energy likelihood is defined: the likelihood of energy measurement $i$
given unknown source $\theta_{u}^{j}$ and constants $\theta_{c}$, prior to
weighting by $w_{j}:$

$p\left(Z_{E}^{i}|\theta_{u}^{j},\theta_{c}\right)=\left.\int\limits_{1.8\mathrm{MeV}}^{11\mathrm{MeV}}\frac{s_{\theta_{u}^{j}}{(E)s_{d}(E)f_{\theta_{c}}\left(E,r_{d}^{j}\right)}}{\mathrm{cdf}_{d}\left(\theta_{u}^{j}\right)}\mathcal{N}\left(Z_{E}^{i},\sigma_{d}^{i}\right)dE\right|_{E=E_{\bar{\nu}}}$
(44)

where the $\mathrm{cdf}_{d}(\theta_{u}^{j})$ merely normalizes the unsmeared
energy spectrum in the numerator into a proper pdf:

$\mathrm{cdf}_{d}\left(\theta_{u}^{j}\right)=\left.\int\limits_{1.8\mathrm{MeV}}^{11\mathrm{MeV}}s_{\theta_{u}^{j}}{(E)s_{d}(E)f_{\theta_{c}}\left(E,r_{d}^{j}\right)}dE\right|_{E=E_{\bar{\nu}}}$
(45)

Matter oscillation effects are assumed to be negligible and, as such, we
ignore the Mikheyev–Smirnov–Wolfenstein (MSW) effect [87]. Following this
simplifying assumption, the $\bar{\nu}_{e}$ survival fraction
$f_{\theta_{c}}\left(E,r_{d}^{j}\right)$ used in Equations 44 and 45 can now
be defined for neutrino sources. At a given range $r$ (in meters) and neutrino
energy $E$ (in MeV), $f_{\theta_{c}}\left(E,r_{d}^{j}\right)$ is defined by:

$\displaystyle
f\left(E,r\right)=1-\cos^{4}\theta_{13}\sin^{2}2\theta_{12}\sin^{2}\frac{1.27\Delta
m^{2}_{12}r}{E}$
$\displaystyle-\sin^{2}2\theta_{13}\cos^{2}\theta_{12}\sin^{2}\frac{1.27\Delta
m^{2}_{23}r}{E}$
$\displaystyle-\sin^{2}2\theta_{13}\sin^{2}\theta_{12}\sin^{2}\frac{1.27\Delta
m^{2}_{31}r}{E}$ (46)

where $\Delta m_{23}^{2}=\Delta m^{2}_{13}-\Delta m^{2}_{12}$. We’ve assumed
the following nominal values for the neutrino mixing angles and mass constants
contained in $\theta_{c}$ per Fogli et al. [63]:

$\displaystyle\theta_{12}=0.566$ $\displaystyle\theta_{13}=0.232$
$\displaystyle\sin^{2}2\theta_{12}=0.856$
$\displaystyle\sin^{2}2\theta_{13}=0.156$ $\displaystyle\Delta
m^{2}_{12}=8.00\times 10^{-5}eV^{2}$ $\displaystyle\Delta
m^{2}_{13}=2.50\times 10^{-3}eV^{2}$ $\displaystyle\Delta
m^{2}_{23}=2.42\times 10^{-3}eV^{2}$

After applying the above constants, Equation 4.1 can be approximated by
Equation 47 which is evaluated in Figure 28:

$f\left(E,r\right)\approx 0.532+0.029\cos\left(\frac{6.147r}{E}\right)+\
0.071\cos\left(\frac{6.350r}{E}\right)+\
0.368\cos\left(\frac{0.203r}{E}\right)$ (47)

Figure 28: $\bar{\nu}_{e}$ survival fraction versus range, produced by the
evaluation of Equation 47. Note that, as reflected in this figure, Equation 47
does not evaluate the solid angle as a function of range. This task is
performed by the mixture weights in Equation 40.

In Equations 44 and 45, $s_{\theta_{u}^{j}}$ is the assumed zero-range
neutrino energy spectrum from unknown source $j$ in the unknown parameter
vector $\theta_{u}$. For a reactor, $s_{\theta_{u}^{j}}$ may be approximated
by:

$s_{\theta_{u}^{j}}\left(E\right)\left[\frac{1}{\mathrm{sec}\times\mathrm{GW_{th}}}\right]\approx
1.407\times\left(1.872\times 10^{20}\right)\left(\frac{1}{3}\right)\times
e^{-(0.3125E+0.25)^{2}}$ (48)

Equation 48, evaluated in Figure 29, uses the 2nd order exponential flux model
first proposed by Vogel and Engel [107] with updated model coefficients. It
assumes 200MeV per fission and 6 antineutrinos of all energies per fission to
get $1.872\times 10^{20}$ antineutrinos per GWth per second, per Bernstein et
al. [68]. The fraction $\frac{1}{3}$ represents the proportion of neutrinos
produced above $E_{\bar{\nu}}\geq 1.8$MeV, per Gil-Botella [80].

The detector IBD cross section $s_{d}$ is shown in Figure 29; it is defined by
the constants in $\theta_{c},$ and approximated by:

$s_{d}(E)\left[\frac{\mathrm{cm}^{2}}{10^{32}\mathrm{p}^{+}}\right]\approx(0.095E-0.123)\times\sqrt{(E-1.293)^{2}-0.261}\times
10^{-42}$ (49)

Figure 29: $\bar{\nu}_{e}$ reactor emission and detection spectra. The
measured spectrum is the product of the emitted and the detected spectra, and
represents the inverse beta decay $\bar{\nu}_{e}$ spectrum that the detector
should expect to measure from a nearby nuclear power plant. The detection
spectra should be independent of the $\bar{\nu}_{e}$ source, and applies to
geo-neutrinos as well as reactor neutrinos.

In Equation 44, $\mathcal{N}(Z_{E}^{i},\sigma_{d}^{i})$ describes the assumed