1307/1307.2832.md

Theoretical Antineutrino Detection, Direction and Ranging at Long Distances

Glenn R. Jocher

Integrity Applications Incorporated, 15020 Conference Center Drive, Chantilly, VA, 20151 USA

Daniel A. Bondy

Integrity Applications Incorporated, 15020 Conference Center Drive, Chantilly, VA, 20151 USA

Brian M. Dobbs

Integrity Applications Incorporated, 15020 Conference Center Drive, Chantilly, VA, 20151 USA

Stephen T. Dye

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

Integrity Applications Incorporated, 15020 Conference Center Drive, Chantilly, VA, 20151 USA

John G. Learned*

Department of Physics and Astronomy
University of Hawaii, Honolulu, HI, 96822 USA

Christopher L. Mulliss

Integrity Applications Incorporated, 15020 Conference Center Drive, Chantilly, VA, 20151 USA

Shawn Usman

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 MW_th 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

---

---

*Corresponding author

Email addresses: glenn.jocher@gmail.com (Glenn R. Jocher), dbondy@integrity-apps.com (Daniel A. Bondy),## Contents

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

## 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 1950s [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 ( $< 1\text{km}$ ) 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}\text{cm}^2$ (see Equation 37). 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.1. 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.
2. 3. Detector External Radioactivity. External radioactivity, such as from $^{40}\text{K}$ 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.
3. 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 $^9\text{Li}$ and $^8\text{He}$ isotopes, and are generally referred to as “cosmogenic” background (see Section 3.3).
4. 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.
5. 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.
6. 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.
7. 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.
8. 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 TW_th [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 accelerator1. 2. Proton decay searches, particularly in the Kaon decay modes, not well studied by SuperKamiokande
2. 3. Search for supernova burst neutrinos, most likely from our galaxy
3. 4. Search for relic neutrinos summed from all past type II supernovae
4. 5. Atmospheric neutrinos, moving beyond SuperKamiokande
5. 6. Solar neutrinos, at least careful monitoring of solar output versus time
6. 7. Reactor neutrino properties, depending upon the proximity and power of reactors
7. 8. Geo-neutrinos both as object of study and background, depending upon location and as discussed below
8. 9. Pion decay at rest neutrinos from a possible portable low energy accelerator (the DAEDALUS idea [98])
9. 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 configurations¹. 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\text{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).

¹The 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.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)n(\backslash lambda)\; =\; 1\; +\; \backslash frac\{B\_1\backslash lambda^2\}\{\backslash lambda^2\; -\; C\_1\}\; +\; \backslash frac\{B\_2\backslash lambda^2\}\{\backslash lambda^2\; -\; C\_2\}\; +\; \backslash frac\{B\_3\backslash lambda^2\}\{\backslash lambda^2\; -\; C\_3\}\; \backslash quad\; (1)$$

Here $B_1 = 1.44$ , $C_1 = 81.76\text{nm}^2$ , $B_2 = B_3 = 0$ , and $C_2 = C_3 = 0\text{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.

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}(\bar{N})(2)N\; =\; \backslash mathcal\{P\}(\backslash bar\{N\})\; \backslash quad\; (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)\backslash bar\{N\}\; =\; y(kB)dE\; \backslash quad\; (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 factor²", 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.

---

²True 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.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.

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 }_{{\lambda }_1}^{{\lambda }_2}f(\lambda )d\lambda (5)\backslash bar\{N\}\; =\; \backslash int\_\{\backslash lambda\_1\}^\{\backslash lambda\_2\}\; f(\backslash lambda)\; d\backslash lambda\; \backslash quad\; (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)\backslash theta\_c(\backslash lambda)\; =\; \backslash arccos\; \backslash left(\; \backslash frac\{1\}\{\backslash beta\; n(\backslash lambda)\}\; \backslash right)\; \backslash quad\; (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{K{E}^2+2m_0{c}^{2}KE}}{KE+m_0{c}^2}(7)\backslash beta\; =\; \backslash frac\{\backslash sqrt\{KE^2\; +\; 2m\_0c^2KE\}\}\{KE\; +\; m\_0c^2\}\; \backslash quad\; (7)$$

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 wavelength³.

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

³No 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.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$ .

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”⁴ 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.

⁴The term “information” here denotes any information which may contribute to a mathematical reconstruction of an event’s location and $dE$ .### 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 modeled⁵.

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^\circ$ N, $137.31^\circ$ 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_{\text{vis}}[\text{MeV}]}$ (Abe et al. [51]).

⁵Concerns 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.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 1150m³ inner region and a 3050m³ 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 $8.9%|{E{\text{vis}}=1\text{MeV}}$ energy resolution that our Detector Model produced for TREND in Table 4.

Comparison Metric	Cited [74]	Our MC
Energy Resolution $1\backslash sigma$ (34% PMT coverage)	$6.2\backslash \%/\backslash sqrt\{E\_\{\backslash text\{vis\}\}[\backslash text\{MeV\}]\}$	$6.24\backslash \%/\backslash sqrt\{E\_\{\backslash text\{vis\}\}[\backslash text\{MeV\}]\}$ (82.4% candidates)
Energy Resolution $1\backslash sigma$ (22% PMT coverage)	$7.3\backslash \%/\backslash sqrt\{E\_\{\backslash text\{vis\}\}[\backslash text\{MeV\}]\}$	$7.69\backslash \%/\backslash sqrt\{E\_\{\backslash text\{vis\}\}[\backslash text\{MeV\}]\}$ (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{\text{vis}}[\text{MeV}]}$ to $6.2%/\sqrt{E_{\text{vis}}[\text{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 GW_th 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_{\text{vis}}[\text{MeV}]}$ (Apollonio et al. [3]). This detector measured the direction of reactor antineutrinos with a $1\sigma$ resolution of $18^\circ$ .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.6m³ inner region and a 22m³ 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 41 - 42, 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$\backslash sigma$	190.0mm	181.9mm
Angular 1$\backslash sigma$, $n=2500$ events	18.0$^\backslash circ$ (19.0$^\backslash circ$ sim)	18.9$^\backslash circ$
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_{\text{vis}}^{e^+} < 8\text{MeV}$ in CHOOZ to $E_{\text{vis}}^{e^+} < 12\text{MeV}$ 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.2m² 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 coverageis naturally a function of the simpler TREND surface area⁶ 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\text{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\text{s}$ time cut (99.8% candidates), 2m range cut (96.5% candidates), and $4\text{MeV} < E{\text{vis}}^{n0} < 12\text{MeV}$ 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.

The mean TREND energy resolution based on 100,000 randomly placed inverse beta decay events (shown in Figure 9) was found to be about $8.9%|{E{\text{vis}}=1\text{MeV}}$ at $E_{\text{vis}}=1\text{MeV}$ . 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_{\text{vis}}[\text{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 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.

⁶Note 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.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	$\backslash sim\; 160,000\backslash text\{m\}^3$
Mass	138kT
Protons	$10^\{34\}\backslash text\{p\}^+$
Attempted PMTs	17000
Achieved PMTs	16906
Optical Coverage	19.9%
PMT Time Resolution
PMT Area	$.203\backslash text\{m\}^2$
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	$\backslash sim\; 2\backslash \%$
Positron energy cut
Positron-wall distance cut	$\backslash |x\_\{e+\}\; -\; x\_\{\backslash text\{wall\}\}\backslash |\; >\; 30\backslash text\{cm\}$
Neutron capture cut	$12\backslash text\{MeV\}\; >\; E\_\{\backslash text\{vis\}\}^\{n0\}\; >\; 4\backslash text\{MeV\}$
Capture energy containment cut
Neutron-wall distance cut	$\backslash |x\_\{n0\}\; -\; x\_\{\backslash text\{wall\}\}\backslash |\; >\; 30\backslash text\{cm\}$
Neutron delay cut
Positron-neutron distance cut
Neutron multiplicity cut	$\backslash text\{num\}(n^0)\; =\; 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.7mm
Positron Annihilation Point Resolution $1\backslash sigma$	374.5mm
Neutron Capture Point Resolution $1\backslash sigma$	243.7mm
Direction Vector Resolution $1\backslash sigma$	461.3mm
Angular Direction $1\backslash sigma$, n=2500 events	37.0°
SNR	0.045
Energy Resolution $1\backslash sigma$	$8.9\backslash \%|\_\{E\_\{\backslash text\{vis\}\}=1\backslash text\{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\text{s}$ from $100\mu\text{s}$ ), positron energy cut ( $E_{\text{vis}}^{e+} < 12\text{MeV}$ from $E_{\text{vis}}^{e+} < 8\text{MeV}$ ), and positron-neutron distance cut (2m from 1m).

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.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 $E_{\text{vis}}=3\text{MeV}$ , 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 $\text{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$ values⁷.

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 8 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$ .

---

⁷A 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.$$C_{NED_d}^{ECEF} = \begin{bmatrix} -\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{bmatrix} \quad (8)$$

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 (9).

$${\bar{\nu }}_e+p\to e^{+}+n(9)\backslash bar\{\backslash nu\}\_e\; +\; p\; \backslash rightarrow\; e^+\; +\; n\; \backslash quad\; (9)$$

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 36 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\text{MeV}$ . These assumptions yield $1.872 \times 10^{20}$ detectable neutrinos per $\text{GW}{\text{th}}$ 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 $300\text{MW}{\text{th}}$ 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\text{ MW}{\text{th}}$ 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 $1063\text{GW}{\text{th}}$ ) and 35 reactor cores (total $93\text{GW}{\text{th}}$ ) 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 $10^{32}$ protons per year (assuming 100% efficiency) due solely to these active reactor cores can be seen in Figure 11.

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.Figure 11: IAEA known reactor background for a $10^{32}\text{p}^+$ detector, saturated at 100 events per year.

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 $^{238}\text{U}$ , $^{232}\text{Th}$ , and $^{40}\text{K}$ 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 $^{238}\text{U}$ , $^{232}\text{Th}$ , and $^{40}\text{K}$ 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.

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)

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 $^{238}\text{U}$ and $^{232}\text{Th}$ .

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 SectionFigure 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.

Figure 13: Crust background for a $10^{32}\text{p}^+$ detector, saturated at 100 events per year.Figure 14: Crust+Mantle background for a $10^{32}\text{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.Figure 16: CRUST 2.0 + Mantle Adjoining Layer density, in $\text{g}/\text{cm}^3$ . Density is seen to generally increase as depth increases. The top ice layer naturally shows the lowest density, slightly under $\text{g}/\text{cm}^3$ , while the bottom Mantle Adjoining Layer density is about $3.4\text{g}/\text{cm}^3$ . Zero-thickness tiles are omitted from this figure.

4.1 these source spectra are combined with the detector cross section, defined in Equation 37 and shown in Figure 29, to create the observed geo-neutrino spectrum for each tile and isotope.

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}\text{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-Atlantic⁸.

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.

⁸The 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.Figure 17: CRUST 2.0 + Mantle Adjoining Layer emissions ( $^{238}\text{U}$ in upper, $^{232}\text{Th}$ 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 $^{232}\text{Th}$ half life of $20.3 \times 10^9$ years and multiplicity of $4\bar{\nu}_e$ per $^{232}\text{Th}$ decay, and $^{238}\text{U}$ half life of $6.45 \times 10^9$ years and multiplicity of $6\bar{\nu}_e$ per $^{238}\text{U}$ 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.

	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 1028	5.30 1028	6.24 1028
C/m3	3.35 1028	3.97 1028	3.79 1028
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
$\backslash mu$-Section (cm2)	1.3 106	0.57 106	4.4 107
$\backslash mu$-Flux (cm-2s-1)	1.6 10-7	0.3 10-7	1.4 10-8
$\backslash mu$-Energy (MeV)	219	276	295
$\backslash mu$-Rate (s-1)	2.13 10-1	1.6 10-2	6.2 10-1
$\backslash mu$-DT (200 $\backslash 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 1032	6.02 1030	1034
$E\_\{\backslash text\{vis\}\}$ threshold [MeV]	0.9	1	0.9	1	2.6
Accidental Rate	80.5$\backslash pm$0.1	0.080$\backslash pm$0.001	277$\backslash pm$3	183$\backslash pm$2	0.3$\backslash pm$0.003
9Li/8He Rate	13.6$\backslash pm$1.0	0.03$\backslash pm$0.02	20$\backslash pm$1	20$\backslash pm$1	18$\backslash pm$1
Fast n Rate	<9 (4.5$\backslash pm$2.6)	<0.05 (0.025$\backslash pm$0.014)	31$\backslash pm$18	31$\backslash pm$18	26$\backslash pm$15
Geo-$\backslash nu$	69.7	2.5$\backslash pm$0.2	2,111	2,088	18
Known Reactor-$\backslash 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].

1. 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.
2. 3. Cosmogenic ⁹Li/⁸He. Cosmogenic isotopes produced by showering muons fall into this category, and are most problematic due to the long lifetimes of ⁸He and ⁹Li 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 ⁹Li/⁸He 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_{\text{vis}} \geq 2.6\text{MeV}$ energy cut due to blurring from the energy measurement resolution.

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 ${}^9\text{Li}/{}^8\text{He}$ count rate. Lasserre et al. scaled the expected cosmogenic ${}^9\text{Li}/{}^8\text{He}$ count rate from KamLAND without accounting for the fact that KamLAND uses a different veto time (2s) than the $\sim 600\text{ms}$ 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 ${}^9\text{Li}/{}^8\text{He}$ 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 37\text{m}$ ) 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 ${}^9\text{Li}/{}^8\text{He}$ 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\text{ m}^3$ fiducial volume of TREND. The discontinuity below $113\text{ m}^3$ occurs because the detector becomes too small to use the $r=3\text{m}$ 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,500\text{m}$ ).

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\text{ }\mu\text{s}$ , $\Delta R=1\text{ m}$ ) than the filter they proposed ( $\tau=30\text{ }\mu\text{s}$ , $\Delta V=1\text{ m}^3$ ). 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\text{ }\mu\text{s}$ , $\Delta R=2\text{ 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}\text{p}^+\text{-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 10) 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_{\text{vis}})=(2.25\times {10}^5)e^{-6.185E_{\text{vis}}}(10)p(E\_\{\backslash text\{vis\}\})\; =\; (2.25\; \backslash times\; 10^5)\; e^\{-6.185\; E\_\{\backslash text\{vis\}\}\}\; \backslash quad\; (10)$$

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 intrinsicFigure 20: Exponential least squares fit to experimental KamLAND data from [79]. Fit was performed using KamLAND data from $E_{\text{vis}} = 0.9\text{MeV}$ to $E_{\text{vis}} = 1.7\text{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_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.

Source	Uncertainty $\backslash pm\; 1\backslash sigma\_\{\backslash text\{syst,\; class\}\}\; \backslash pm\; 1\backslash sigma\_\{\backslash text\{syst,\; source\}\}$	Number of Sources	Comment
Known reactor flux	$\backslash pm\; 2.0\backslash \%\; \backslash pm\; 3.4\backslash \%$	471 reactor cores	does not include sterile $\backslash nu$ uncertainty
Crust flux	$\backslash pm\; 20.0\backslash \%\; \backslash pm\; 8.0\backslash \%$	$\backslash sim\; 129600$ tiles	100% correlated across $^\{238\}\backslash text\{U\}$ and $^\{232\}\backslash text\{Th\}$
Mantle flux	$\backslash pm\; 50.0\backslash \%\; \backslash pm\; 0\backslash \%$	$\backslash sim\; 50396$ tiles	100% correlated across $^\{238\}\backslash text\{U\}$ and $^\{232\}\backslash text\{Th\}$
Cosmogenic $^9\backslash text\{Li\}/^8\backslash text\{He\}$	$\backslash pm\; 3.3\backslash \%\; \backslash pm\; 0\backslash \%$	4584/detector	KamLAND [50] extrapolation
Fast neutron flux	$\backslash pm\; 10.0\backslash \%\; \backslash pm\; 0\backslash \%$	4584/detector	KamLAND [50] (capped at 10.0%)
Accidental	$\backslash pm\; 1.3\backslash \%\; \backslash pm\; 0\backslash \%$	4584/detector	Borexino [72, 56] isopurity extrapolation
Detector proton count	$\backslash pm\; 0.0\backslash \%\; \backslash pm\; 2.4\backslash \%$	1/detector	i.e. fiducial volume uncertainty
$\backslash Delta\; m\_\{12\}^2$	$\backslash begin\{matrix\}\; +2.4\backslash \%\; \backslash \backslash \; -2.8\backslash \%\; \backslash end\{matrix\}\; \backslash pm\; 0\backslash \%$	$\backslash sim\; 180467$	per Fogli et al. [63]
$\backslash Delta\; m\_\{13\}^2$	$\backslash begin\{matrix\}\; +3.4\backslash \%\; \backslash \backslash \; -2.6\backslash \%\; \backslash end\{matrix\}\; \backslash pm\; 0\backslash \%$	$\backslash sim\; 180467$	per Fogli et al. [63]
$\backslash sin^2\; \backslash theta\_\{12\}$	$\backslash begin\{matrix\}\; +5.2\backslash \%\; \backslash \backslash \; -5.7\backslash \%\; \backslash end\{matrix\}\; \backslash pm\; 0\backslash \%$	$\backslash sim\; 180467$	per Fogli et al. [63]
$\backslash sin^2\; \backslash theta\_\{13\}$	$\backslash begin\{matrix\}\; +10.0\backslash \%\; \backslash \backslash \; -9.8\backslash \%\; \backslash end\{matrix\}\; \backslash pm\; 0\backslash \%$	$\backslash 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.

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_{\text{syst, class}}$ ) as well as a systematic uncertainty which is independent for each source in the class ( $\sigma_{\text{syst, source}}$ ). For continuous background sources (non-neutrinos), where $\sigma_{\text{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 TREND⁹.

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 optimal¹⁰ estimator is presented, as well as two tractable suboptimal approximations appropriate for real

⁹“Poor” directional resolution of TREND defined as $\sim 0.045$ SNR per Section 2.4.3.

¹⁰Optimal 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.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 \mathrm{\mid }Z)=\frac{p(Z\mathrm{\mid }\theta )p(\theta )}{p(Z)}=\frac{1}{c}p(Z\mathrm{\mid }\theta )p(\theta )(11)p(\backslash theta|Z)\; =\; \backslash frac\{p(Z|\backslash theta)p(\backslash theta)\}\{p(Z)\}\; =\; \backslash frac\{1\}\{c\}p(Z|\backslash theta)p(\backslash theta)\; \backslash quad\; (11)$$

where $c = p(Z)$ acts as a normalization constant¹¹, 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)(12)Z\; =\; H(\backslash theta,\; d)\; +\; v(d)\; \backslash quad\; (12)$$

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.

$$\theta =\left[\begin{array}{c}{\theta }_u\\ {\theta }_k\\ {\theta }_c\end{array}\right](13)\backslash theta\; =\; \backslash begin\{bmatrix\}\; \backslash theta\_u\; \backslash \backslash \; \backslash theta\_k\; \backslash \backslash \; \backslash theta\_c\; \backslash end\{bmatrix\}\; \backslash quad\; (13)$$

$${\theta }_u=\left[\begin{array}{c}{\theta }_{u_1}\\ \mathrm{⋮}\\ {\theta }_{u_{n_u}}\end{array}\right](14)\backslash theta\_u\; =\; \backslash begin\{bmatrix\}\; \backslash theta\_\{u\_1\}\; \backslash \backslash \; \backslash vdots\; \backslash \backslash \; \backslash theta\_\{u\_\{n\_u\}\}\; \backslash end\{bmatrix\}\; \backslash quad\; (14)$$

$${\theta }_k=\left[\begin{array}{c}{\theta }_{k_1}\\ \mathrm{⋮}\\ {\theta }_{k_{n_k}}\end{array}\right](15)\backslash theta\_k\; =\; \backslash begin\{bmatrix\}\; \backslash theta\_\{k\_1\}\; \backslash \backslash \; \backslash vdots\; \backslash \backslash \; \backslash theta\_\{k\_\{n\_k\}\}\; \backslash end\{bmatrix\}\; \backslash quad\; (15)$$

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$ ¹². For $n_u$ unknown reactors:

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

¹¹In 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.

¹²Note 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.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$ . $^{232}\text{Th}$ and $^{238}\text{U}$ isotope geo-neutrino sources in the mantle and crust are considered here, though potassium $^{40}\text{K}$ 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\text{MeV}$ . The combined flux from all neutrino sources for a $10^{32}\text{p}^+$ detector are shown in Figure 22.

We have eight known categories of background in $\theta_k$ :

$${\theta }_k=\left[\begin{array}{c}{\theta }_{\text{IAEA}}\\ {\theta }_{{\text{crust}}^{238}\text{U}}\\ {\theta }_{{\text{crust}}^{232}\text{Th}}\\ {\theta }_{{\text{mantle}}^{238}\text{U}}\\ {\theta }_{{\text{mantle}}^{232}\text{Th}}\\ {\theta }_{\text{accidental}}\\ {\theta }_{\text{fast }{n}^0}\\ {\theta }_{\text{cosmogenic}}\end{array}\right](17)\backslash theta\_k\; =\; \backslash begin\{bmatrix\}\; \backslash theta\_\{\backslash text\{IAEA\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{crust\}^\{238\}\backslash text\{U\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{crust\}^\{232\}\backslash text\{Th\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{mantle\}^\{238\}\backslash text\{U\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{mantle\}^\{232\}\backslash text\{Th\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{accidental\}\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{fast\; \}\; n^0\}\; \backslash \backslash \; \backslash theta\_\{\backslash text\{cosmogenic\}\}\; \backslash end\{bmatrix\}\; \backslash quad\; (17)$$

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 frames¹³.

Some sources in $\theta_k$ are assumed to be point sources, such as the known reactors in $\theta_{\text{IAEA}}$ . Others are clearly volume sources, such as the mantle and crust $^{238}\text{U}$ and $^{232}\text{Th}$ 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 submerged in the waters of the North Atlantic.

¹³Note that these low level coordinate system transformations are implicitly assumed in all the equations in this paper.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 $(4 \times 0.5^{-30})^\circ \times (4 \times 0.5^{-30})^\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 options¹⁴ 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 }_{-\mathrm{\infty }}^{\mathrm{\infty }}p(x,i)p(i)di(18)p(x)\; =\; \backslash int\_\{-\backslash infty\}^\{\backslash infty\}\; p(x,\; i)p(i)di\; \backslash quad\; (18)$$

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.

¹⁴Approximation is a third option, used in our suboptimal approximations in Section 4.2.