Also at ]Tel Aviv Univ. formerly of ]The Univ. of Tennessee

# Rock Neutron Backgrounds from FNAL Neutrino Beamlines
in the $\nu$BDX-DRIFT Detector

D. Aristizabal Sierra<EMAIL_ADDRESS>Universidad Técnica Federico
Santa María - Departamento de Física
Casilla 110-V, Avda. España 1680, Valparaíso, Chile  J. L. Barrow
<EMAIL_ADDRESS>The Massachusetts Institute of Technology, Department of
Physics, 77 Massachusetts Avenue, Building 4, Room 304, Cambridge, MA 02139,
USA  [  [  B. Dutta<EMAIL_ADDRESS>Mitchell Institute for Fundamental
Physics and Astronomy, Department of Physics and Astronomy, Texas A&M
University, College Station, TX  D. Kim<EMAIL_ADDRESS>Mitchell
Institute for Fundamental Physics and Astronomy, Department of Physics and
Astronomy, Texas A&M University, College Station, TX  D. Snowden-Ifft
<EMAIL_ADDRESS>Physics Department, Occidental College, 1600 Campus Rd., Los
Angeles, CA 90041  L. Strigari<EMAIL_ADDRESS>Mitchell Institute for
Fundamental Physics and Astronomy, Department of Physics and Astronomy, Texas
A&M University, College Station, TX  M. H. Wood<EMAIL_ADDRESS>Department
of Quantitative Sciences, Canisius College, 2001 Main St., Buffalo, NY

###### Abstract

The $\nu$BDX-DRIFT collaboration seeks to detect low-energy nuclear recoils
from CE$\nu$NS or BSM interactions at FNAL. Backgrounds due to rock neutrons
are an important concern. We present a GENIE and GEANT4 based model to
estimate backgrounds from rock neutrons produced in neutrino-nucleus
interactions within the rock walls surrounding the underground halls. This
model was bench-marked against the $2009$ COUPP experiment performed in the
MINOS hall in the NuMI neutrino beam, and agreement is found between
experimental results and the modeled result to within $30\%$. Working from
this validated model, a similar two-stage simulation was performed to estimate
recoil backgrounds in the $\nu$BDX-DRIFT detector across several beamlines. In
the first stage utilizing GEANT4, neutrons were tallied exiting the walls of a
rectangular underground hall utilizing four different neutrino beam
configurations. These results are presented for use by other underground
experiments requiring estimations of their rock neutron backgrounds. For
$\nu$BDX-DRIFT, the second stage propagated neutrons from the walls and
recorded energy deposited within a scintillator veto surrounding the detector
and nuclear recoils within the detector’s fiducial volume. The directional
signal from the $\nu$BDX-DRIFT detector allows additional background
subtraction. A sample calculation of a $10\,$m${}^{3}\cdot\,$yr exposure to
the NuMI Low Energy (LE) beam configuration shows a CE$\nu$NS signal-to-noise
ratio of $\sim$2.5.

††preprint: APS/123-QED

##  I Introduction

The $\nu$BDX-DRIFT detector is a directional time projection chamber (TPC)
suitable for measurements of nuclear recoils produced by coherent elastic
neutrino-nucleus scattering (CE$\nu$NS) [1, 2] and by new physics interactions
within the neutrino and dark-sectors, including those such as light (MeV) dark
matter (DM) [3]. Its directional capabilities offer a unique environment for
the identification of beyond Standard Model (BSM) signals [4]. The detector
can operate with a variety of target nuclei, e.g. H, C, S and possibly Pb [3].
Studies of the performance of the detector using decay-in-flight neutrinos
produced in the Long Baseline Neutrino Facility (LBNF) beamline at Fermi
National Accelerator Laboratory (FNAL) [5] have been presented in Ref. [3].
These results have demonstrated that, with reasonable exposures
($10\,\text{m}^{3}$ for $7$ years of data taking), the detector will be able
to measure $\sim 300$-$400$ CE$\nu$NS events across various target materials.
The resulting large statistics will in turn enable measurements of Standard
Model (SM) electroweak and nuclear parameters, as well as searches for
neutrino non-standard interactions (NSI), among others.

After the first measurements of CE$\nu$NS using CsI and liquid argon (LAr)
detectors by the COHERENT collaboration [6, 7] at Oak Ridge National
Laboratory’s Spallation Neutron Source (ORNL SNS), an effort to undertake
further measurements across other target nuclei and different energy spectra
utilizing various neutrino sources continues globally [8]. Low energy
experiments using reactor neutrinos are underway [9, 10, 11, 12, 13, 14, 15],
as well as further experiments at the ORNL SNS [16]; this includes planning
stages for the SNS Second Target Station, along with the European Spallation
Source [17]. As a part of this global effort, the $\nu$BDX-DRIFT detector can
provide a new and complementary avenue if it was to be based at FNAL: it would
utilize decay-in-flight neutrinos and thereby observe higher energy regimes
than the other global suite of experiments. Further, its technology offers
measurements of the angular spectrum, in addition to the recoil energy
spectrum; thus, in principle, cross section measurements in kinematic
variables pertaining to the nuclear recoil are possible.

However, the viability of all the above depends critically on background
levels. Neutrino-induced neutrons produced in the rock, so-called “rock
neutrons”, produce recoil-like backgrounds which are problematic and occupy
the majority of discussions within this paper. The rock neutrons can be
produced directly from neutrino-nucleus collisions, or when other neutrino-
nucleus end-state particles interact in the surrounding material, generating
still more neutrons. As will be shown below, rock neutrons produced in these
ways have energies up to $\sim 100\,$MeV and can produce nuclear recoils $\sim
100\,$keV which themselves are expected from CE$\nu$NS and BSM interactions
[3]. Recoils produced inside the shielding material around $\nu$BDX-DRIFT were
considered in Ref. [3], where it was shown that an expected signal-to-
background ratio of better than $23$ could be achieved. Rock neutrons produced
in the much larger volume of rock surrounding the underground facilities at
FNAL are harder to estimate as the calculation must convolve the neutrino
energy spectrum and interaction cross section on a variety of nuclei, the
propagation of all end-state particles through the rock to the experimental
hall, the possible interactions with shielding surrounding the detector, and,
finally, the generation of nuclear recoils inside the fiducial volume of the
detector.

The procedure presented here relies first upon a Monte Carlo neutrino event
generator package, GENIE [18], accounting for interactions of the neutrino
beam with the rock material in the surrounding walls of the FNAL underground
MINOS experimental hall [19]. This first step is followed by a GEANT4 [20]
simulation, which accounts for the propagation of the end-state particles
generated in the GENIE calculation and which potentially can enter the
detector fiducial volume. The procedure is bench-marked with the aid of the
COUPP beam-tagged data, which provides information on neutron-induced nuclear
recoils. Four independent simulations will be presented based on four
different neutrino flux configurations (NuMI LE and HE modes [21] as well as
DUNE on-axis and $39\,$m off-axis [22]), and so collectively provide
information not only valuable for a potential $\nu$BDX-DRIFT physics program
but also for future neutrino detectors at FNAL. The results to be presented
here can thus be understood as being aligned with and complementary to current
efforts at the Accelerator Neutrino Neutron Interaction Experiment (ANNIE) at
FNAL [23]. Finally, results will be presented for rock neutron backgrounds in
the fiducial volume of the $\nu$BDX-DRIFT with strong background protections
afforded from the surrounding scintillator and the directionality of the
interaction.

The remainder of this paper is organized as follows. In Sec. II we provide a
detailed discussion of the physics capabilities of the $\nu$BDX-DRIFT
detector. In Sec. III, details of the beam-tagged COUPP data are presented. In
Sec. III.1, the inputs used in the GENIE-GEANT4 Monte Carlo simulations are
given. Results of the GENIE output for final state particles are presented,
along with the nuclear recoil spectrum in the COUPP detector’s fiducial
volume. In Sec. IV, the neutron energy, zenith and azimuth spectra are
provided for all four simulations, while in Sec. V these results will be used
as input for the determination of the neutron background in the $\nu$BDX-DRIFT
detector fiducial volume. Finally, in Sec. VII, a summary and conclusions will
be presented.

##  II Physics capabilities of the $\nu$BDX-DRIFT detector

Measurements of CE$\nu$NS within the $\nu$BDX-DRIFT detector will provide data
enabling: (i) the determination of SM parameters, and (ii) searches for new
interactions in the neutrino sector. These measurements can also enable
searches for MeV-scale DM candidates produced in collisions of a proton beam
on a fixed target. Detection proceeds by observation of the nuclear recoils
produced by either of these progenitors within the fiducial volume of the
detector.

Focusing on (i), the measurements which can be carried out include a precision
determination of the weak mixing angle at $\sqrt{Q^{2}}\simeq 100\,$MeV, and
the determination of the neutron root-mean-square (rms) radius of nuclides for
which no data is yet available. As for (ii), searches include NSIs,
interactions mediated by light vectors and scalars, along with sterile
neutrinos. Analysis of these types of interactions have been completed using
COHERENT and other reactor CE$\nu$NS data (see e.g. [24, 25, 26, 27, 28]).
Results from $\nu$BDX-DRIFT will thus prove complementary, while testing these
hypotheses in a different energy domain and with different detector
technologies.

As a function of detector operation pressure, CE$\nu$NS event rates in CS2
peak at about $400\,$Torr. For a $10\,\text{m}^{3}$ detector operating over
$7\,$years, the expected rate is on the order of $400\,$events. For CF4 and
utilizing the same operation pressure, the event yield increases by about a
factor of two. With C8H20Pb, although with a lead target, the event yield is
smaller because of the rapid loss of coherence. However, the statistics
combined with the detector features are still large enough for the analysis of
a few physics cases. Demanding isolation of lead-induced events, to study lead
nuclear properties, fixes the operation pressure in that case to $\sim 5$ Torr
[3].

Using CF4 (C8H20Pb) as material target, a $10\,\text{m}^{3}$ detector operated
at the pressures mentioned above will be able to measure the carbon and
fluorine (lead) neutron rms with a $\sim 3$% ($\sim 5$%) precision. Ref. [3]
has reported the following $1\sigma$ measurements

$\displaystyle r_{\text{rms}}^{n}|_{\text{C}}$
$\displaystyle=2.84^{+0.13}_{-0.15}\,\text{fm}\ ,$ $\displaystyle
r_{\text{rms}}^{n}|_{\text{Pb}}$
$\displaystyle=5.50^{+0.30}_{-0.29}\,\text{fm}\ .$ (1)

Measurements for carbon and fluorine through electroweak neutral current
processes do not exist, so these results provide valuable information for a
better understanding of nuclear properties of light nuclide. For lead the
result is not as competitive as that derived from PREX measurements [29, 30],
but can be understood as complementary to it.

Figure 1: The 2009 COUPP bubble formation data tagged to the beam pulse.
Published here with the permission of the COUPP collaboration.

Studies of the weak mixing angle in CS2 and CF4 result in the following
$1\sigma$ measurements

$\displaystyle\sin^{2}\theta_{W}|_{\text{CS}_{2}}$
$\displaystyle=0.238^{+0.020}_{-0.016}\ ,$
$\displaystyle\sin^{2}\theta_{W}|_{\text{CF}_{4}}$
$\displaystyle=0.238^{+0.021}_{-0.017}\ ,$ (2)

both for $\sqrt{Q^{2}}\subset[78,397]\,$MeV, a renormalization scale for which
at present no data is available. Interestingly enough, these results exceed
what so far COHERENT measurements have achieved (see e.g. [24, 31]) and are
competitive with those expected from DUNE using the electron channel [32].

Searches for NSI in CS2 can explore muon flavor related effective couplings.
Sensitivities can improve by about a factor 2-3 upon current limits. To a
certain extent they are not very sensitive to backgrounds (assuming reasonable
amounts) nor to quark flavor. The $1\sigma$ measurements that can be achieved
are given by [3],

$\displaystyle\epsilon_{\mu\mu}$
$\displaystyle=[-0.013,0.011]\oplus[0.30,0.32]\ ,$
$\displaystyle\epsilon_{e\mu}$ $\displaystyle=[-0.064,0.064]\ .$ (3)

As has been emphasized, in order to achieve these goals a detailed
understanding of rock neutron backgrounds becomes mandatory. The following
sections focus on that.

Composition in rock at FNAL
---
Isotope |  ${}^{1}_{1}$H |  ${}^{12}_{6}$C |  ${}^{16}_{8}$O |  ${}^{23}_{11}$Na |  ${}^{27}_{13}$Al |  ${}^{28}_{14}$Si |  ${}^{39}_{19}$K |  ${}^{40}_{20}$Ca |  ${}^{56}_{26}$Fe
Composition [$\%$] | 1.5 | 1.1 | 56.4 | 0.3 | 9.5 | 24.2 | 0.9 | 4.3 | 1.8

Input parameters used in the simulations
---
Beamline & Mode | (POT/Pulse)$\times 10^{13}$ | (Inter/Pulse/$\text{m}^{3})\times 10^{-4}$ | Period [s]
NuMI LE (c. 2009) | $2.88$ | $204.42$ | $2.43$
NuMI LE | $4.00$ | $283.92$ | $1.3$
NuMI HE | $4.00$ | $1277.69$ | $1.3$
DUNE On-Axis at $1.2\,$MW | $7.5$ | $1142.23$ | $1.2$
DUNE $39\,$m Off-Axis at $1.2\,$MW | $7.5$ | $9.89$ | $1.2$

Table 1: Upper: The percentages of various nuclear isotopes in the rock, taken
from discussions with FNAL experts. Lower: Summary of the input parameters for
the models considered in this paper. The numbers of POT per pulse for NuMI and
DUNE have been taken from Refs. [5, 21].

##  III COUPP

In order to present reliable results for nuclear recoil background predictions
within the $\nu$BDX-DRIFT detector, any simulation used to predict such
backgrounds requires bench-marking against data. Fortunately, such data
exists. In 2009, the COUPP DM collaboration performed an experiment in the
MINOS hall on-axis to an active NuMI beam [33] at FNAL. COUPP was a bubble
chamber experiment with a $15$-$20\,$keV threshold for detecting nuclear
recoils filled with 3.5 kg of CF3I [33]. As discussed in [33], COUPP was a
threshold detector providing no information on recoil energy or particle
(nucleus) identification. Additionally, COUPP had no sensitivity to $\beta$,
$\gamma$, or minimum ionizing particles. Using acoustic information $\alpha$
particle discrimination was possible [33].

In 2009, events were tagged as occurring when the beam was on or not. For the
DM data analysis, only events uncorrelated with the beam were analyzed and
published. However, unpublished, beam-tagged data from the COUPP collaboration
was obtained [34]; a summary of these findings can be seen in Fig. 1. The pink
data points are single, fiducial events not tagged as $\alpha$ particles and
are interpreted here as nuclear recoil events. The average of these data–taken
from September 27, 2009 to November 8, 2009–is $4.65\pm
0.19\,$events$/$kg$\cdot$day. During this running period, the cosmic veto was
not operational; thus, some fraction of these events were caused by non-beam-
related particles. To estimate this background, non-beam-related, background
data taken during this time were averaged. Using a $100\,$ms timing window;
the background rate due to random coincidences was estimated to be $0.0863\pm
0.0074\,$events$/$kg$\cdot$day. Subtracting this from the observed rate gives,
a true, beam-related nuclear recoil rate of $4.56\pm
0.19\,$events$/$kg$\cdot$day to be compared to predictions.

###  III.1 The Model

The parameters and model for backgrounds in the COUPP 2009 exposure to the
neutrino beam are presented here. The composition of the rock can be seen in
Table 1 (upper Table), and was assumed to be at a density of $2.33\,$g/cm3.
From the FNAL Data Logger [35], the average number of protons on target (POT)
per pulse was $2.88\times 10^{13}$ with an average period of $2.43\,$s. These
parameters as well as other assumed parameters are summarized in Table 1
(lower Table). The neutrino flux at the COUPP location was taken from [36] and
increased by a factor of $(1040/939)^{2}$ due to the upstream location of the
COUPP experiment relative to the originally assumed location [36]. Fig. 2
shows the resultant flux, alongside several others to be discussed below.

According to MINOS logs [37], the NuMI beam was in reverse horn current mode
during the COUPP 2009 run, implying predominately antineutrino production
during the run period. Given the on-axis nature of the COUPP detector, it is
expected that few differences exist between the $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ fluxes (horn current settings) across the various NuMI
beam energy settings [38]. Despite $\nu_{\mu}$ contamination of the
$\overline{\nu}_{\mu}$ beam at high energies, we consider this single
neutrino-type approximation robust, especially given the comparative lack of
neutrons (which yield the most background events) entering the final state via
charged current $\nu_{\mu}$ interactions.

Figure 2: The $\nu_{\mu}$ energy spectra for various locations at FNAL: Fluxes
at $1040\,$m downstream at NuMI in the LE and HE mode and at DUNE at $574\,$m
downstream as well for the on-axis and the off-axis $39\,$m configurations.
Results for NuMI are adapted from Ref. [36], while for DUNE from the DUNE
Technical Design Report (Fig. 4.9) [22]. For the purposes of this study, small
deviations in shape and rate between the $\nu_{\mu}$ and
$\overline{\nu}_{\mu}$ horn modes spectra are ignored and are utilized
identically.

###  III.2 GENIE Event Generation

Given the previously discussed inputs, simulation of primary particle
production via NuMI $\overline{\nu}_{\mu}$ interactions within the rock
surrounding the COUPP detector could be undertaken. Neutral and charged
current processes across the whole range of energies of the NuMI flux
resulting from $\overline{\nu}_{\mu}$ scattering were considered, providing
predictions for final state neutrons, protons, charged and neutral pions, and
antimuons. Fig. 3 shows energy distributions of the six different final state
particles considered in this model for the NuMI LE neutrino flux employed in
the COUPP simulation.

Figure 3: Energy spectra for $n$, $p$, $\pi^{-}$, $\pi^{+}$, $\mu^{+}$ and
$\pi^{0}$ end-states of $\overline{\nu}_{\mu}$-nucleus interactions obtained
by a GENIE Monte Carlo simulation, for the NuMI LE neutrino flux. These
spectra are used as input for the GEANT4 simulation of the COUPP result.

These primary particle production simulations were completed using the GENIE
Monte Carlo event generator [18], a staple within the FNAL neutrino community.
The G18_10a GENIE tune [39] was used as a baseline, and cross section splines
for all constituent elements were produced across the whole NuMI LE energy
range. The chosen tune utilizes the hA2018 final state interaction
(intranuclear cascade) model [40, 41], which uses a table-based method to
predict full final states. A similar simulation was undertaken using the
hN2018 final state interaction model, which employs a fully stochastic
intranuclear cascade and generally provides final state predictions with
higher final state nucleon multiplicities. The mixture of elements making up
the rock served as a direct input to GENIE for event production, creating
single samples; generally, the samples used throughout the studies discussed
here were $\sim 10^{6}$ events in size. Histograms with $\sim 50\,$MeV/c
binning were constructed for the 6 most abundant final state particle types,
$n$, $p$, $\pi^{-}$, $\pi^{+}$, $\mu^{+}$ and $\pi^{0}$. As an example Figure
3 shows the energy distributions for these 6 end-state particles for the NuMI
LE configuration and the hA GENIE model. These distributions were use to as
inputs for GEANT4 111Correlated, event-by-event simulation of primary
interaction products is indeed possible, and future work will utilize such
techniques..

Figure 4: Upper graph: The labeled geometry of the underground experimental
hall. Lower graph: A GEANT4 simulation showing the location of the detector
relative to the walls. The dimensions of the underground hall are 480/1070/427
cm in $x/y/z$. The aqua color shows the fiducial volume of the $\nu$BDX-DRIFT
detector. The white frames show the location of the scintillator. Purple lines
show neutrons trajectories. Yellow shows electron trajectories.

###  III.3 GEANT4 Propagation

GEANT4 [42] was used to propagate the end-state particles from GENIE through
the rock and into the experimental hall and detector shown in Fig. 4. The
dimensions of this hall (chosen to roughly approximate the size of the hallway
where the COUPP experiment occurred) were considered small enough that uniform
generation of end-state particles was assumed. The source considered in these
simulations was taken as the rock walls, whose thickness was increased up to
$2\,$m, at which point the observed rates in the detector stabilized. The
COUPP detector was modeled as a cylindrical fiducial volume 15 cm in diameter
and 12 cm high filled with CF3I. This was surrounded on almost all sides with
propylene glycol (C3H8O$2$) the exception being a water filled region above
the CF3I. The outer dimensions of the these elements were 30 cm in diameter
and 44 cm high. Again we thank members of the COUPP collaboration for
providing this information [34]. All massive nuclear recoils in the CF3I were
analyzed. Fig. 5 shows the resulting nuclear recoil spectrum in nuclear mass.

Figure 5: The spectrum of recoiling nuclei with kinetic energies ($E_{r}$)
greater than 16.8 keV. The small number of isotopes at masses other than C, F
or I natural abundances are due to inelastic collisions between, mostly
neutrons, and the target nuclei. Three regions in recoil masses are
identified. Recoil masses in the region labelled “$\alpha$ discrimination”
were not counted because of $\alpha$ discrimination. Recoils in the region
labelled “C and F recoil region” were treated similarly, see text. Recoils in
the third region, “Iodide recoil region”, were treated similarly as well. The
shaded C, F and I regions are largely arbitrary, but there exist effectively
no events within them beyond those at and slightly below the expected masses
of these species. See text for further details.

The nucleation efficiency for bubble formation following nuclear recoil within
the COUPP detector is given [43] as,

$\epsilon(E)=1-e^{-\alpha[(E-E_{T})/E_{T})]}\quad(E>E_{T})\ ,$ (4)

where $E_{T}$ is a universal threshold while $\alpha$ depends on the recoil
type; $\alpha_{\text{CF}}$ (for Carbon and Fluorine recoils) was determined to
be $0.15$ from AmBe neutron exposures, while
$\alpha_{\text{I}}=2.8^{+1.6}_{-0.8}$ (for Iodide recoils) and
$E_{T}=16.8^{+0.8}_{-1.1}\,\text{keV}$ were determined using a $12\,$GeV
$\pi^{-}$ beam [43]. For this work, the mean values of these quantities were
employed; note that no uncertainty was given for $\alpha_{CF}$.

GEANT4 events in which multiple bubbles were removed as the COUPP data reports
only single events in the fiducial volume. GENIE’s input simulation to GEANT4
utilizing the hA2018 model yields a predicted rate of $2.930\pm
0.039\,$events$/$kg$\cdot$day. As a check on the effect of the geometry of the
experimental hall on this result the length of the experimental hall was
increased by a factor of 3. The result was $2.890\pm
0.046\,$events$/$kg$\cdot$day in agreement with the previous result. For
clarity these results, and the ones discussed below, are summarized in Table
2. The GENIE hN model yields a rate of $3.081\pm
0.025\,$events$/$kg$\cdot$day. These were averaged together to produce a
predicted rate of $3.006\pm 0.023\,$events$/$kg$\cdot$day. These events were
created by, largely, rock neutrons entering the COUPP detector from the walls,
thus creating recoils which nucleated a bubble.

Recoils can of course also be created directly inside the COUPP fiducial
volume by direct neutrino scatters, the dominant component of these being non-
CE$\nu$NS events such as neutrino-nucleon quasi-elastic scattering, a
subdominant contribution from neutrino-nucleon scattering, resonant single
pion production and by products of deep-inelastic scattering. To better
understand this, GENIE was run with CF3I, instead of rock, as the target, and
an overall rate for such scatters was $0.35\,$events$/$kg$\cdot$day. However,
this total event count ignores the fact that not all such events will nucleate
a bubble. For some events, no large remnant nuclei survive; for those that do
survive, there is a less than 100% chance of nucleating a bubble given their
momentum222Note that GENIE is currently unable to record all the properties of
remnant nuclei; similarly, for all but one nucleus (Oxygen), no photonic de-
excitation occurs. There is motion within the community to include more of
this necessary microphysics [44, 45], and we look forward to more updates to
such tools.. We therefore bracket our modeled results as $(3.006,3.356)\pm
0.023\,$events$/$kg$\cdot$day. These event rates are to be compared to the
experimental rate of $4.56\pm 0.19\,$events$/$kg$\cdot$day.

The predicted rate of this study sits roughly $30\%$ lower than the observed
experimental rate. There are, however, a large number of systematics which
could explain this difference. The bubble formation model has systematics
associated with the assumptions discussed above, though these appear to be
relatively small. For instance, varying the bubble formation parameters such
as $\alpha_{I}$ and $E_{T}$ gives a $0.09\,$events/kg$\cdot$day systematic
variation to the rate. GENIE and GEANT have systematics associated with the
particular models chosen, and are largely unknown to this study without the
use of a universe style approach. Slight changes in the geometric
configuration of the detector can also contribute to the uncertainty.
Similarly, the neutrino flux model is known to have large normalization
uncertainties which have not been considered for this study.

Rate Comparison Summary
---
Source | Rate [events/kg$\cdot$day]
GENIE hA | $2.930\pm 0.039$
GENIE hA w/3$\times$longer exp. hall | $2.890\pm 0.046$
GENIE hN | $3.081\pm 0.025$
GENIE hA, hN average | $3.006\pm 0.023$
Unshielded in-situ |  $0$ to $0.35$
Prediction | $(3.006,3.356)\pm 0.023$
Experiment | $4.56\pm 0.19$
Table 2: This table summarizes the rates from various sources and, at the end,
the final prediction range in comparison with the COUPP data. Number of
simulated particles
---
Beamline & Mode | Stage I [$\times$106] | Walls [$\times$106] | Stage II [$\times$109]
NuMI LE | 207 | 17.2 | 2.36
NuMI HE | 130 | 2.66 | 1.70
DUNE On-Axis at $1.2\,$MW | 434 | 8.26 | 2.36
DUNE $39\,$m Off-Axis at $1.2\,$MW | 1660 | 5.51 | 2.10
Table 3: Output table number of particles simulated at various stages. Column
2 shows the number of end-state particles simulated in Stage I (see Sec. IV).
Column 3 shows the number of neutrons entering the experimental hall from the
walls. These neutrons were used to generate the distributions for the Stage II
simulations (see Sec. V). Column 4 shows the number of neutrons simulated in
Stage II restarted on the walls of the experimental hall. (see Sec. V).

##  IV Stage I: Rock Neutron Results

With the bench-marked model in hand we now turn to predicting backgrounds in
future, planned experiments. As the COUPP results show, backgrounds due to
rock neutrons in an unshielded detector are high, too high to accomplish the
goals of the $\nu$BDX-DRIFT collaboration. We therefore include a
scintillating veto around the simulated $\nu$BDX-DRIFT detector. The COUPP
collaboration installed a scintillating veto around most of their detector
with a resulting drop in un-vetoed rate after the period of unshielded running
described above and shown in Figure 1. That the rate did not drop further was
the result of lack of shielding around the bottom of the detector; the
shielding was designed to veto cosmic-ray generated events not beam events.
For purposes of simulation we will assume the $\nu$BDX-DRIFT detector is
surrounded by 75 cm of BC-521 organic scintillator on all sides, similar to
the veto COUPP utilized. As will be shown below use of this veto drastically
reduces the rate of events in the $\nu$BDX-DRIFT detector.

But as a result the simple, single-stage simulation used for the COUPP
background calculation is impractical. A two-stage strategy was therefore
adopted in which neutrons were recorded exiting the walls of the experimental
hall. The hall was assumed to have an upstream and downstream wall
perpendicular to the neutrino beamline and 4 walls parallel to the beamline as
shown in Fig. 4. For each wall the energy and angular distributions of
neutrons exiting the walls for the first time were recorded and smoothed. In a
second stage, neutrons were restarted at the walls with the same energy and
angular distributions with a resulting increase in simulation speed of roughly
two orders of magnitude. The computed energy and angular distributions for all
simulations are shown below for use in other applications.

To bracket the range of possibilities at FNAL four simulations were done.
Table 1 (lower) summarizes the main input parameters for these simulations.
The neutrino energy spectra for all simulations are shown in Fig. 2. All
simulations assumed that the horn currents were set to predominantly produce
$\overline{\nu}_{\mu}$s. $\overline{\nu}_{\mu}$s produce more neutrons than
$\nu_{\mu}$s due to the nature of the charge current interaction, and, in
terms of background, therefore represent a worst case scenario. The location
of the COUPP detector was on the far upstream end of the MINOS hall, 939 m
from the target. All NuMI simulations were done at this location. As before
the fluxes for NuMI, from [36] assuming 1040 m from target, were increased by
(1040/939)2 to correct for this assumption. For the DUNE simulations the
experimental hall, shown in Fig. 4, was located 574 m from the DUNE target at
the location of the DUNE near detector hall. Two positions were chosen, on-
axis and 39 m off-axis, to bracket the possibilities there. As shown in Fig. 2
these positions have very different fluxes and energy spectra. Note that for
the DUNE simulations it was assumed that the experimental hall shown in Fig. 4
was completely surrounded by rock which is not what is planned for the near
detector hall. The DUNE simulations, therefore, are more indicative of
backgrounds generated on either side of the DUNE near detector hall. The total
number of neutrino interactions per m3 of rock per pulse is shown in column 3
of Table 1 (lower) for each beam and mode.

Figure 6: The energy distribution of rock-neutrons generated in the 4
simulations of Table 1. The blue lines show the spectra coming from the
upstream wall. The orange lines show the spectra coming from the side walls.
And the green lines show the spectra coming from the downstream wall. In the
left graph, the solid (dashed) curves correspond to results obtained with the
NuMI HE (LE) neutrino mode. In the right graph—instead—to results derived with
the DUNE on-axis (off-axis)

configuration.

Simulations output for neutron flux from the walls
---
Beamline & Mode | Upstream [$n^{0}$/s/m2] | Sides [$n^{0}$/s/m2] | Downstream [$n^{0}$/s/m2] | Background [events/m3/year]
NuMI LE | 0.0355 | 0.0204 | 0.0110 | 8.61 $\pm$ 0.62
NuMI HE | 0.209 | 0.131 | 0.0727 | 54.9 $\pm$ 3.8
DUNE On-Axis at $1.2\,$MW | 0.101 | 0.0276 | 0.0524 | 23.3 $\pm$ 1.3
DUNE $39\,$m Off-Axis at $1.2\,$MW | 0.000381 | 0.0000831 | 0.000162 | 0.0396 $\pm$ 0.0031
Table 4: Output table shows neutron flux from different walls and background
in the signal region. For details see Sec. IV.

Figure 7: Upper row graphs: The zenith angle distribution of rock-neutrons
generated in the NuMI HE (left graph) and the DUNE on-axis (right graph)
simulations of Table 1. The blue lines show the spectra coming from the
upstream wall. The yellow lines show the spectra coming from the side walls.
And the green lines show the spectra coming from the downstream wall. With
rather small variations, results for the NuMI LE (DUNE off-axis) resemble
those of the NuMI HE (DUNE on-axis) as so are not displayed. Lower row graphs:
Same as for those on top, but for azimuth angle distribution. Results are
presented for the same simulations as we have found that differences as well
with the other two are negligible.

End-state particles from these interactions were propagated, using GEANT4, to
the walls of the experimental hall where, as discussed above, neutron
characteristics were recorded and saved. Charged particles exiting the walls
of the experimental hall were not saved as they would either range out in the
scintillator or be vetoed there. Table 3 shows the number of particles
simulated at each stage of the simulation. The smoothed, rock-neutron energy
distributions for the four simulations are shown in Fig. 6. As expected the
flux of neutrons exiting the walls is higher on the upstream wall than the
downstream wall with a harder spectrum. The sides fall somewhere in between.
Also for the same POT/pulse, see Table 1, higher energy configurations produce
higher fluxes of rock-neutrons. Table 4 shows a summary of the output from the
simulations. Columns 2, 3 and 4 show the rates for various surfaces relative
to the beam. These numbers are nothing more than the integral of the
differential flux, see Fig. 6, with energy but they provide a simple way of
comparing the various beamlines and modes.

Fig. 7 (upper row) shows the spectra of zenith angles, measured from the
$z$-axis, for each of the walls. As expected the upstream wall shows a more
pronounced peak than does the downstream wall. Results are shown only for the
NuMI HE mode and the DUNE on-axis configuration. Results for the NuMI LE (DUNE
off-axis 39 m) resemble rather closely those of the NuMI HE mode (DUNE on-
axis) and so are not displayed. Fig. 7 (lower row) shows as well the spectra
of azimuth angles, measured from the $x$-axis, for each of the walls. The
zenith and azimuth angle specifies a vector which, adopting the GEANT
convention, points in a direction from which the particle came. The upstream
wall therefore emits particles with azimuth angles from 0 to $\pi$, vectors
which point into the rock, while the downstream wall emits particles from
$\pi$ to $2\pi$, vectors which point into the experimental hall. Once again
the upstream wall exhibits a more concentrated distribution as the emission of
neutrons from the downstream wall would entail multiple bounces before
emission from the wall. Finally the sides, right hand wall shown here, shows
an asymmetric distribution skewed towards smaller azimuth angles indicating a
preference for emission from the beam direction. In summary the angular
distributions show a preference for neutron emission from the direction to the
target which decreases from the upstream wall to the sides to the downstream
wall.

Figure 8: Plots showing the distribution of recoil energies vs energy
deposited in the scintillator with the neutrino-induced end-state particle
responsible for the recoil shown in different colors. The left graph shows the
results for C recoils while the right graph shows the S recoils. Both are
heavily dominated by neutron end-states (about 63% for both target nuclei).
The vertical dashed black lines indicate the recoil thresholds, 75 keV for C
and 200 keV for S. The horizontal dashed black lines show the threshold for
the scintillator veto, 1 MeV; events with larger energies are vetoed. The
lower right region therefore shows the signal region where either CE$\nu$NS or
BSM recoils events would occur. The background rate, in events per m3 per
year, are shown there. The sum is shown in the fifth column of Table 4.

##  V Stage II: $\nu$BDX-DRIFT Results

As discussed above the main motivation for this work is the reliable
prediction of backgrounds for the $\nu$BDX-DRIFT experiment. To that end a
Stage II simulation was set up and run to predict backgrounds. Neutrons were
fired from the walls of the experimental hall with energy and angular spectra
such as shown in Figures 6 and 7. From the outside in, the detector consisted
of a $75\,$cm thick BC-521 scintillator veto surrounding the entire detector
with outer dimensions of $3\,$m, a $0.5\,$inch thick stainless-steel, cubic
vacuum vessel with outer dimensions of $1.5\,$m and a cubic fiducial volume
for recoils composed of CS2 at a density 2.44 times higher than $400\,$Torr.
This increased pressure increases the efficiency for recording recoils while
minimizing double recoils [46]; final results are corrected at the end.

GEANT recorded any energy deposited in the scintillator veto and in the
fiducial volume. Fig. 8 shows the results for the NuMI LE beamline and mode.
On the horizontal axis is the recoil kinetic energy for C and S. On the
vertical axis is the amount of energy deposited in the scintillator. The
different colors represent the end-state particles from
$\overline{\nu}$-nucleus interactions which produced neutrons which entered
the experimental hall and created C or S recoils in the fiducial volume of the
$\nu$BDX-DRIFT detector. Neutron end-state particles from
$\overline{\nu}$-nucleus interactions dominate the recoil rate. The vertical
dashed line shows the kinetic energy threshold for recoil detection after [3].
As can be seen in these graphs a huge number of recoils are predicted above
threshold. However the vast majority of nuclear recoils above threshold also
come with an enormous deposition of energy in the scintillator, on order
$100\,$MeV 333It should be noted that the benchmarked COUPP 2009 experiment
was mostly sensitive to 1-10 MeV neutrons while $\nu$BDX-DRIFT is mostly
sensitive to 10-100 MeV neutrons due to the necessity of penetrating the
scintillator.. These large energy depositions occur due to showers produced as
the neutrons traverse the detector and resulting charged particle interactions
in the scintillator veto. The horizontal dashed line indicates a $1\,$MeV
threshold on the veto; events with energy greater than this are vetoed. Signal
events, CE$\nu$NS events or BSM interactions, would appear in the lower right
corner of these graphs. Backgrounds, in this context, means events due to beam
neutrons appearing in this lower right corner. The rate of recoils, and
errors, for C and S appear in this lower right corner in Fig. 8 in units of
events per m3 per year. The fifth column in Table 4 shows the background rates
for each of the beamlines and modes studied in this paper. As can be seen the
highest backgrounds occur in the NuMI beamline in the HE mode.

As a check of the Stage I of the simulation for this high background
configuration, a run  _with the scintillating veto in place_ was completed.
After firing $2.3\times 10^{9}$ end-state neutrons from the walls the result
was in statistical agreement with the Stage II neutron results to within
$15\%$ validating the use of the multi-stage procedure.

There remains a question as to how these beam-related backgrounds compare to
their non-beam-related cousins. While this question has not been studied in
detail, an indication can be found when again considering the 2009 COUPP
results [33]. The COUPP collaboration found a neutron background of $3$ events
across a $28.1\,$kg$\cdot$day exposure for a rate of about
$0.1\,$events/kg$\cdot$day; this rate was measured with lower thresholds and
while maintaining a scintillating shield similar to that described in this
work. However, this rate was not in coincidence with the beam. We can estimate
to an order of magnitude that $10\,\mu$s timing resolution is possible, giving
an approximate additional 10-5 reduction in background from non-beam-related
sources occurring during a beam-spill for a total rate of about $\sim
10^{-6}\,$events/kg$\cdot$day, or $\sim 6\times
10^{-4}\,$events/m${}^{3}\cdot$yr. This rate is much smaller than any of those
predicted in Table 4.

Figure 9: Left graph: Neutron and CE$\nu$NS zenith angle distribution as a
function of zenith angle in degrees. The result has been derived assuming the
NuMI LE neutrino flux, with parameters as specified in Table 1. As expected,
the CE$\nu$NS signal peaks at $90^{\circ}$ while the neutron-induce recoils
have a much wider spread (see text in Sec. VI for details). The histograms for
different maximum recoil energies show that events pile up with increasing
energy. Right graph: Neutron and CE$\nu$NS recoil energy spectra as a function
of nuclear recoil energy. The result has been derived with the same
assumptions that those used for the left graph. The different energy lines are
correlated with the zenith angle histogram in the left graph and graphically
indicate the number of events that for that energy have been piled up in the
zenith angle distribution peak.

##  VI Signal and rock neutrons backgrounds

Our results demonstrate that un-vetoed, rock neutron backgrounds can be
substantial, in particular for the NuMI HE mode and the DUNE on-axis
configuration. Further discrimination of the CE$\nu$NS signal against this
background would be helpful. To do so the directional capabilities of the
detector can be employed. Information from the neutron and CE$\nu$NS zenith
angle distribution spectra combined with their recoil energy spectra provide
information that allows—in principle—efficient background discrimination. The
CE$\nu$NS angular distribution is expected to peak in the direction
perpendicular to the neutrino flux. This can be readily understood from the
fact that the recoil (zenith) angle $\theta_{r}$ and recoil energy $E_{r}$ are
related through [4]

$\cos\theta_{r}=\sqrt{\frac{m_{N}\,E_{r}}{2}}\left(\frac{1}{E_{\nu}}+\frac{1}{m_{N}}\right)\,,$
(5)

where $\theta_{r}$ is the recoil angle relative to the direction of the
neutrino, $m_{N}$ is the mass of the nucleus and $E_{\nu}$ is the energy of
the neutrino. For the typical recoil energies ($<1\,$MeV), induced by a “high-
energy” neutrino beam ($\sim$GeV) as those we have consider in these
simulations, lead to small $\cos\theta_{r}$. For CE$\nu$NS this translates
into most events clustering at $90^{\circ}$, independent of the neutrino beam
we choose. To exploit this fact the neutron zenith angle distribution has to
be as well categorized. Its exact morphology, in contrast to the CE$\nu$NS
signal, does depend on the neutrino flux and so for concreteness we have
performed calculations for the NuMI LE mode.

The left graph in Fig. 9 shows the results for both spectra for a 10 m3 year
exposure. The neutron recoil angular distribution has a mild tendency to
cluster at about $90^{\circ}$ due to a tendency of rock neutrons to preserve
the forward direction of the beam. However, their spectrum has a much wider
spread in comparison to that of neutrinos recoils. This result thus shows that
with a reasonable angular resolution further discrimination ($\sim$104:1
altogether, scintillating veto plus angular cuts) of background events is
possible. At $90^{\circ}$ the signal-to-background ratio is estimated by
comparing the number of events at peak, is $\sim$2.5.

The recoil energy spectra provide, as well, useful discrimination power. To
determine the degree to which by itself, or through its interplay with zenith
angle spectral information this can be done, we have calculated the CE$\nu$NS
recoil energy signal as well as neutron recoil energy spectra for the same
neutrino flux configuration. Results are shown in the right graph in Fig. 9.
The CE$\nu$NS signal spreads over a wider energy range (compared to its
clustering at $90^{\circ}$) but does peak towards lower recoil energies. The
rock neutron background peaks as well at low recoil energies, but in contrast
to the CE$\nu$NS signal does populate the full energy range suggesting a
different spectrum which could be exploited.

In addition some amount of C and S recoil discrimination is present. The
difference in these spectra could be used to further discriminate the signals.
More work is needed to fully exploit the background rejection capability of
these signatures.

Other backgrounds could be considered and studied. The decay-in-flight
neutrino beam energies extend up to and even beyond $\sim 10\,$GeV; thus, in
addition to CE$\nu$NS, other higher-energy processes such as quasielastic,
resonance, and deep-inelastic scattering will occur, see e.g. [47]. The cross
sections for these higher-energy interactions (wherein the constituent
nucleons become the system’s dominant degrees of freedom) are sizable at
higher $Q^{2}$. As discussed above for COUPP, these type of events occur at a
rate of $0.35\,$events$/$kg$\cdot$day. The $\nu$BDX-DRIFT detector with a mass
of 1.6 kg will see these events on the order of 1 per day. In terms of
backgrounds to $\nu$BDX-DRIFT in searches for CE$\nu$NS and BSM nuclear
recoils though, the large neutrino energies generally imply high particle
multiplicities and are comparatively unique in their topologies. For instance,
charged particles produced in conjunction with nuclear recoils can be rejected
as signal events. As shown above, the scintillating veto is extremely
effective at rejecting neutrals at these large energies. Additionally, as
events like this will be present in the data, their characteristics can be
measured and studied themselves, an interesting topic it’s own right.

##  VII Conclusion

In this paper we have studied rock neutron backgrounds in the $\nu$BDX-DRIFT
detector. Rock neutrons are produced by the interaction of neutrinos with the
rock surrounding the underground hall where the detector is deployed. End-
state particles produced in these interactions come from a GENIE Monte Carlo
calculation which uses four possible neutrino fluxes (NuMI LE and HE modes and
DUNE on-axis and off-axis 39 m configurations) interacting with the rock
composed mainly of Oxygen, Silicon, Aluminum and Iron. The energy spectra of
the final state particles produced in these interactions serve then as an
input for a GEANT4 Monte Carlo simulation, which propagates these states
throughout the rock and so allows the characterization of the neutrons
emerging from the walls of the hall. These neutrons are then used to study the
possible backgrounds to which the $\nu$BDX-DRIFT detector will be subject to
while being operated at the FNAL.

The simulation is bench-marked against the 2009 beam-tagged COUPP data,
obtained by the COUPP collaboration during operation in the MINOS hall while
the NuMI beamline was operated in the LE mode. Agreement between the simulated
and actual data is found within 30%. After this validation, results for
energy, zenith and azimuth spectra for the neutrons emitted by the walls are
reported. These results, crucial for the determination of rock neutron
backgrounds in the $\nu$BDX-DRIFT detector, are as well useful for future
neutrino experiments at the FNAL. They add to undergoing efforts by the ANNIE
collaboration, which aims to characterize neutron backgrounds at the FNAL.

With the “morphology” of the emitted neutrons at hand, rock neutron
backgrounds within the $\nu$BDX-DRIFT fiducial volume have been determined. By
assuming the detector to be fully surrounded by a 75 cm thick BC-521
scintillator veto, for the four different neutrino flux configurations we have
found that the DUNE off-axis 39 m provides the most background-suppressed
experimental scenario. Rock neutron backgrounds gradually increase from the
NuMI LE to the DUNE on-axis to the NuMI HE, with the latter being the
configuration leading to the largest background. Detailed results have been
reported in Table 4.

Finally we have discussed discrimination of rock neutron backgrounds against
CE$\nu$NS signals. Using NuMI LE as a representative case, we have compared
neutron and CE$\nu$NS zenith and recoil energy spectra. The results
demonstrate that discrimination against rock neutron backgrounds is possible.
Firstly, the CE$\nu$NS signal peaks at $90\degree$, in contrast to the neutron
background that spreads more uniformly. At peak, the signal-to-background
ratio has been roughly estimated to be $\sim 2.5$. Information from the recoil
energy spectra shows that background-free energy windows exist, thus offering
an experimental avenue for CE$\nu$NS measurements as well as for BSM searches.

## Acknowledgements

We thank Andrew Sonnenschein and Jeter Hall for useful conversations and
providing us the the COUPP 2009 results upon which this work is based. We
thank Eric Vázquez-Jáuregui for carefully reading the manuscript and useful
comments, particularly concerning the actual geometry of the COUPP 2009
detector. D.A.S. is supported by ANID grant “Fondecyt Regular”
N${}^{\text{o}}$ 1221445\. BD and LES acknowledge support from DOE Grant de-
sc0010813.

## References

  * Freedman [1974] D. Z. Freedman, Phys. Rev. D 9, 1389 (1974).
  * Freedman  _et al._ [1977] D. Z. Freedman, D. N. Schramm, and D. L. Tubbs, Ann. Rev. Nucl. Part. Sci. 27, 167 (1977).
  * Aristizabal Sierra  _et al._ [2021] D. Aristizabal Sierra, B. Dutta, D. Kim, D. Snowden-Ifft, and L. E. Strigari, Phys. Rev. D 104, 033004 (2021).
  * Abdullah  _et al._ [2020] M. Abdullah, D. Aristizabal Sierra, B. Dutta, and L. E. Strigari, Phys. Rev. D 102, 015009 (2020), arXiv:2003.11510 [hep-ph] .
  * Strait  _et al._ [2016] J. Strait _et al._ (DUNE), (2016), arXiv:1601.05823 [physics.ins-det] .
  * Akimov  _et al._ [2017] D. Akimov _et al._ (COHERENT), Science 357, 1123 (2017), arXiv:1708.01294 [nucl-ex] .
  * Akimov  _et al._ [2021] D. Akimov _et al._ (COHERENT), Phys. Rev. Lett. 126, 012002 (2021), arXiv:2003.10630 [nucl-ex] .
  * Abdullah  _et al._ [2022] M. Abdullah _et al._ , (2022), arXiv:2203.07361 [hep-ph] .
  * Aguilar-Arevalo  _et al._ [2019] A. Aguilar-Arevalo _et al._ (CONNIE), Phys. Rev. D 100, 092005 (2019), arXiv:1906.02200 [physics.ins-det] .
  * Agnolet  _et al._ [2017] G. Agnolet _et al._ (MINER), Nucl. Instrum. Meth. A 853, 53 (2017), arXiv:1609.02066 [physics.ins-det] .
  * Strauss  _et al._ [2017] R. Strauss _et al._ , Eur. Phys. J. C 77, 506 (2017), arXiv:1704.04320 [physics.ins-det] .
  * Akimov  _et al._ [2020] D. Y. Akimov _et al._ (RED-100), JINST 15, P02020 (2020), arXiv:1910.06190 [physics.ins-det] .
  * Colaresi  _et al._ [2022] J. Colaresi, J. I. Collar, T. W. Hossbach, C. M. Lewis, and K. M. Yocum, (2022), arXiv:2202.09672 [hep-ex] .
  * Billard  _et al._ [2017] J. Billard _et al._ , J. Phys. G 44, 105101 (2017), arXiv:1612.09035 [physics.ins-det] .
  * Bonet  _et al._ [2021] H. Bonet _et al._ (CONUS), Phys. Rev. Lett. 126, 041804 (2021), arXiv:2011.00210 [hep-ex] .
  * Akimov  _et al._ [2022] D. Akimov _et al._ , in  _2022 Snowmass Summer Study_ (2022) arXiv:2204.04575 [hep-ex] .
  * Garoby  _et al._ [2018] R. Garoby _et al._ , Phys. Scripta 93, 014001 (2018).
  * Andreopoulos  _et al._ [2010] C. Andreopoulos _et al._ , Nucl. Instrum. Meth. A 614, 87 (2010), arXiv:0905.2517 [hep-ph] .
  * Ambats  _et al._ [1998] I. Ambats _et al._ (MINOS), NUMI-L-337, FERMILAB-DESIGN-1998-02 (1998), 10.2172/1861363.
  * Agostinelli  _et al._ [2003a] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003a).
  * Adamson  _et al._ [2016] P. Adamson _et al._ , Nucl. Instrum. Meth. A 806, 279 (2016), arXiv:1507.06690 [physics.acc-ph] .
  * Abi  _et al._ [2020] B. Abi _et al._ (DUNE), (2020), arXiv:2002.03005 [hep-ex] .
  * Anghel  _et al._ [2015] I. Anghel _et al._ (ANNIE), (2015), arXiv:1504.01480 [physics.ins-det] .
  * Papoulias and Kosmas [2018] D. K. Papoulias and T. S. Kosmas, Phys. Rev. D 97, 033003 (2018), arXiv:1711.09773 [hep-ph] .
  * Aristizabal Sierra  _et al._ [2019a] D. Aristizabal Sierra, V. De Romeri, and N. Rojas, JHEP 09, 069 (2019a), arXiv:1906.01156 [hep-ph] .
  * Aristizabal Sierra  _et al._ [2019b] D. Aristizabal Sierra, B. Dutta, S. Liao, and L. E. Strigari, JHEP 12, 124 (2019b), arXiv:1910.12437 [hep-ph] .
  * Aristizabal Sierra  _et al._ [2022] D. Aristizabal Sierra, V. De Romeri, and D. K. Papoulias, (2022), arXiv:2203.02414 [hep-ph] .
  * Coloma  _et al._ [2020] P. Coloma, I. Esteban, M. C. Gonzalez-Garcia, and M. Maltoni, JHEP 02, 023 (2020), [Addendum: JHEP 12, 071 (2020)], arXiv:1911.09109 [hep-ph] .
  * Abrahamyan  _et al._ [2012] S. Abrahamyan _et al._ , Phys. Rev. Lett. 108, 112502 (2012), arXiv:1201.2568 [nucl-ex] .
  * Horowitz  _et al._ [2014] C. J. Horowitz, K. S. Kumar, and R. Michaels, Eur. Phys. J. A 50, 48 (2014), arXiv:1307.3572 [nucl-ex] .
  * Miranda  _et al._ [2020] O. G. Miranda, D. K. Papoulias, G. Sanchez Garcia, O. Sanders, M. Tórtola, and J. W. F. Valle, JHEP 05, 130 (2020), [Erratum: JHEP 01, 067 (2021)], arXiv:2003.12050 [hep-ph] .
  * de Gouvea  _et al._ [2020] A. de Gouvea, P. A. N. Machado, Y. F. Perez-Gonzalez, and Z. Tabrizi, Phys. Rev. Lett. 125, 051803 (2020), arXiv:1912.06658 [hep-ph] .
  * Behnke  _et al._ [2011] E. Behnke, J. Behnke, S. J. Brice, D. Broemmelsiek, J. I. Collar, P. S. Cooper, M. Crisler, C. E. Dahl, D. Fustin, J. Hall, J. H. Hinnefeld, M. Hu, I. Levine, E. Ramberg, T. Shepherd, A. Sonnenschein, M. Szydagis, and C. Collaboration, Phys. Rev. Lett. 106 (2011), 10.1103/PhysRevLett.106.021303.
  * COUPP Collaboration [2020] COUPP Collaboration, Private Communication (2020).
  * [35] Z. Yuan, “D44 - data logger plotter,” .
  * Kopp [2007] S. E. Kopp, (2007), arXiv:0709.2737 [hep-ex] .
  * [37] Collaboration (MINOS), “Minos run logs,” .
  * Aliaga Soplin [2016] L. Aliaga Soplin,  _Neutrino Flux Prediction for the NuMI Beamline_ , Ph.D. thesis, William-Mary Coll. (2016).
  * Andreopoulos  _et al._ [2015] C. Andreopoulos, C. Barry, S. Dytman, H. Gallagher, T. Golan, R. Hatcher, G. Perdue, and J. Yarba, (2015), arXiv:1510.05494 [hep-ph] .
  * Niewczas and Sobczyk [2019] K. Niewczas and J. T. Sobczyk, Phys. Rev. C 100, 015505 (2019), arXiv:1902.05618 [hep-ex] .
  * Golan  _et al._ [2012] T. Golan, C. Juszczak, and J. T. Sobczyk, Phys. Rev. C 86, 015505 (2012), arXiv:1202.4197 [nucl-th] .
  * Agostinelli  _et al._ [2003b] S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003b).
  * Behnke  _et al._ [2013] E. Behnke, T. Benjamin, S. J. Brice, D. Broemmelsiek, J. I. Collar, P. S. Cooper, M. Crisler, C. E. Dahl, D. Fustin, J. Hall, C. Harnish, I. Levine, W. H. Lippincott, T. Moan, T. Nania, R. Neilson, E. Ramberg, A. E. Robinson, M. Ruschman, A. Sonnenschein, E. Vazquez-Jauregui, R. A. Rivera, L. Uplegger, and C. Collaboration, Phys. Rev. D 88 (2013), 10.1103/PhysRevD.88.021101.
  * Gardiner [2021a] S. Gardiner, Phys. Rev. C 103, 044604 (2021a), arXiv:2010.02393 [nucl-th] .
  * Gardiner [2021b] S. Gardiner, Comput. Phys. Commun. 269, 108123 (2021b), arXiv:2101.11867 [nucl-th] .
  * Burgos  _et al._ [2007] S. Burgos, J. Forbes, C. Ghag, M. Gold, V. A. Kudryavtsev, T. B. Lawson, D. Loomba, P. Majewski, D. Muna, A. S. Murphy, G. G. Nicklin, S. M. Paling, A. Petkov, S. J. S. Plank, M. Robinson, N. Sanghi, N. J. T. Smith, D. P. Snowden-Ifft, N. J. C. Spooner, T. J. Sumner, J. Turk, and E. Tziaferi, ASTROPARTICLE PHYSICS 28, 409 (2007).
  * Formaggio and Zeller [2012] J. A. Formaggio and G. P. Zeller, Rev. Mod. Phys. 84, 1307 (2012), arXiv:1305.7513 [hep-ex] .