of 14
June 28, 2022
Ionization yield measurement in a germanium CDMSlite detector using
photo-neutron sources
M.F. Albakry,
1, 2
I. Alkhatib,
3
D.W.P. Amaral,
4
T. Aralis,
5
T. Aramaki,
6
I.J. Arnquist,
7
I. Ataee Langroudy,
8
E. Azadbakht,
8
S. Banik,
9, 10
C. Bathurst,
11
D.A. Bauer,
12
L.V.S. Bezerra,
1, 2
R. Bhattacharyya,
8
M.A. Bowles,
13
P.L. Brink,
14
R. Bunker,
7
B. Cabrera,
15
R. Calkins,
16
R.A. Cameron,
14
C. Cartaro,
14
D.G. Cerde ̃no,
4, 17
Y.-Y. Chang,
18
M. Chaudhuri,
9, 10
R. Chen,
19
N. Chott,
13
J. Cooley,
16
H. Coombes,
11
J. Corbett,
20
P. Cushman,
21
F. De Brienne,
22
M.L. di Vacri,
7
M.D. Diamond,
3
E. Fascione,
20, 2
E. Figueroa-Feliciano,
19
C.W. Fink,
18
K. Fouts,
14
M. Fritts,
21
G. Gerbier,
20
R. Germond,
20, 2
M. Ghaith,
20
S.R. Golwala,
5
J. Hall,
23, 24
B.A. Hines,
25
M.I. Hollister,
12
Z. Hong,
3
E.W. Hoppe,
7
L. Hsu,
12
M.E. Huber,
25, 26
V. Iyer,
9, 10
A. Jastram,
8
V.K.S. Kashyap,
9, 10
M.H. Kelsey,
8
A. Kubik,
23
N.A. Kurinsky,
14
R.E. Lawrence,
8
M. Lee,
8
A. Li,
1, 2
J. Liu,
16
Y. Liu,
1, 2
B. Loer,
7
P. Lukens,
12
D. MacDonell,
1, 2
D.B. MacFarlane,
14
R. Mahapatra,
8
V. Mandic,
21
N. Mast,
21
A.J. Mayer,
2
H. Meyer zu Theenhausen,
27, 28
́
E. Michaud,
22
E. Michielin,
1, 2
N. Mirabolfathi,
8
B. Mohanty,
9, 10
J.D. Morales Mendoza,
8
S. Nagorny,
20
J. Nelson,
21
H. Neog,
21
V. Novati,
19
J.L. Orrell,
7
M.D. Osborne,
8
S.M. Oser,
1, 2
W.A. Page,
18
R. Partridge,
14
D.S. Pedreros,
22
R. Podviianiuk,
29
F. Ponce,
7
S. Poudel,
29
A. Pradeep,
1, 2
M. Pyle,
18, 30
W. Rau,
2
E. Reid,
4
R. Ren,
19
T. Reynolds,
3
A. Roberts,
25
A.E. Robinson,
22
T. Saab,
11
B. Sadoulet,
18, 30
I. Saikia,
16
J. Sander,
29
A. Sattari,
3
A. Scarff,
1, 2
B. Schmidt,
19
R.W. Schnee,
13
S. Scorza,
23, 24
B. Serfass,
18
D.J. Sincavage,
21
C. Stanford,
15
J. Street,
13
F.K. Thasrawala,
28
D. Toback,
8
R. Underwood,
20, 2
S. Verma,
8
A.N. Villano,
25
B. von Krosigk,
27, 28
S.L. Watkins,
18
O. Wen,
5
Z. Williams,
21
M.J. Wilson,
27, 3
J. Winchell,
8
K. Wykoff,
13
S. Yellin,
15
B.A. Young,
31
T.C. Yu,
14
B. Zatschler,
3
S. Zatschler,
3
A. Zaytsev,
27, 28
E. Zhang,
3
L. Zheng,
8
and S. Zuber
18
1
Department of Physics & Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada
2
TRIUMF, Vancouver, BC V6T 2A3, Canada
3
Department of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada
4
Department of Physics, Durham University, Durham DH1 3LE, UK
5
Division of Physics, Mathematics, & Astronomy,
California Institute of Technology, Pasadena, CA 91125, USA
6
Department of Physics, Northeastern University,
360 Huntington Avenue, Boston, MA 02115, USA
7
Pacific Northwest National Laboratory, Richland, WA 99352, USA
8
Department of Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, TX 77843, USA
9
School of Physical Sciences, National Institute of Science Education and Research, Jatni - 752050, India
10
Homi Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai 400094, India
11
Department of Physics, University of Florida, Gainesville, FL 32611, USA
12
Fermi National Accelerator Laboratory, Batavia, IL 60510, USA
13
Department of Physics, South Dakota School of Mines and Technology, Rapid City, SD 57701, USA
14
SLAC National Accelerator Laboratory/Kavli Institute for Particle Astrophysics and Cosmology, Menlo Park, CA 94025, USA
15
Department of Physics, Stanford University, Stanford, CA 94305, USA
16
Department of Physics, Southern Methodist University, Dallas, TX 75275, USA
17
Instituto de F ́ısica Te ́orica UAM/CSIC, Universidad Aut ́onoma de Madrid, 28049 Madrid, Spain
18
Department of Physics, University of California, Berkeley, CA 94720, USA
19
Department of Physics & Astronomy, Northwestern University, Evanston, IL 60208-3112, USA
20
Department of Physics, Queen’s University, Kingston, ON K7L 3N6, Canada
21
School of Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
22
D ́epartement de Physique, Universit ́e de Montr ́eal, Montr ́eal, Qu ́ebec H3C 3J7, Canada
23
SNOLAB, Creighton Mine #9, 1039 Regional Road 24, Sudbury, ON P3Y 1N2, Canada
24
Laurentian University, Department of Physics, 935 Ramsey Lake Road, Sudbury, Ontario P3E 2C6, Canada
25
Department of Physics, University of Colorado Denver, Denver, CO 80217, USA
26
Department of Electrical Engineering, University of Colorado Denver, Denver, CO 80217, USA
27
Institute for Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT), 76344, Germany
28
Institut f ̈ur Experimentalphysik, Universit ̈at Hamburg, 22761 Hamburg, Germany
29
Department of Physics, University of South Dakota, Vermillion, SD 57069, USA
30
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
31
Department of Physics, Santa Clara University, Santa Clara, CA 95053, USA
(Dated: June 28, 2022)
Two photo-neutron sources,
88
Y
9
Be and
124
Sb
9
Be, have been used to investigate the ionization
yield of nuclear recoils in the CDMSlite germanium detectors by the SuperCDMS collaboration.
This work evaluates the yield for nuclear recoil energies between 1 keV and 7 keV at a temperature
arXiv:2202.07043v2 [physics.ins-det] 27 Jun 2022
2
of
50 mK. We use a Geant4 simulation to model the neutron spectrum assuming a charge yield
model that is a generalization of the standard Lindhard model and consists of two energy dependent
parameters. We perform a likelihood analysis using the simulated neutron spectrum, modeled
background, and experimental data to obtain the best fit values of the yield model. The ionization
yield between recoil energies of 1 keV and 7 keV is shown to be significantly lower than predicted by
the standard Lindhard model for germanium. There is a general lack of agreement among different
experiments using a variety of techniques studying the low-energy range of the nuclear recoil yield,
which is most critical for interpretation of direct dark matter searches. This suggests complexity
in the physical process that many direct detection experiments use to model their primary signal
detection mechanism and highlights the need for further studies to clarify underlying systematic
effects that have not been well understood up to this point.
I. INTRODUCTION
Many direct detection experiments search for dark
matter particles that elastically scatter off the nucleus
of a target material. The detectors themselves are sen-
sitive to the recoiling nucleus and the dark matter in-
teraction is inferred through detection of the nuclear re-
coil [1]. The measured nuclear recoil energy (
E
r
) from
such an interaction depends on the momentum transfer
between the dark matter particle and the nucleus. A
precise measurement of the recoil spectrum provides in-
formation about the mass of the putative dark matter
particle. In semiconductor targets, such as germanium
(Ge) or silicon (Si), the recoiling nucleus transfers its
energy to ionization (
e
/h
+
pairs) and phonons. The ra-
tio of the ionization energy to the total recoil energy is
commonly referred to as the “ionization yield”. For de-
tectors that rely on deriving the nuclear recoil energy
from the ionization energy, a precise understanding of
the ionization yield is needed [2]. A deeper knowledge
of ionization yield is not only important for direct dark
matter searches, but also for calibrating the nuclear re-
coil signal induced by coherent elastic neutrino nucleus
scattering [3–5].
Figure 1 shows the results from a number of experi-
ments that have measured the ionization yield in Ge as
a function of nuclear recoil energy. The figure also pro-
vides a comparison between the experimental results and
a theoretical model of the ionization yield in Ge which
is derived by Lindhard et al. [22] and which is the ref-
erence model in the field of direct dark matter searches.
Experiments often use a value of the
k
parameter (ex-
plained in Sec. V A), whose physical origin is related to
the electronic stopping power [24], to better match the
Lindhard model to their experimental results. The Lind-
hard model, however, does not account for some factors
that may be important at low energy, such as atomic
binding [7], electric field dependence [25], and tempera-
ture dependence [26]. Thus the Lindhard model is not
expected to be an accurate prediction of the ionization
yield at the low recoil energies (
<
10 keV) expected from
low mass dark matter particles (
<
5 GeV/c
2
). Besides
the shortcoming of the Lindhard model at low energies,
it can also be noted from Fig. 1 that there exists some dis-
agreement in the ionization yield measurements by vari-
ous experiments in this recoil energy regime.
1
10
2
10
3
10
Nuclear Recoil Energy (keV)
0
0.1
0.2
0.3
0.4
0.5
0.6
Ionization Yield
Chasman (77 K)
Messous (77 K)
Sattler (77 K)
Baudis (77 K)
Texono (77 K)
CoGeNT (77 K)
Jones (77 K)
CDMS-II (50 mK)
Simon (35 mK)
Shutt (25 mK)
EDELWEISS ’07 (17 mK)
Lindhard k=0.157
Scholz
EDELWEISS ’17
yield model
}
Collar ’21
FIG. 1. Literature data on ionization yield measurements as a
function of recoil energy in germanium [6–21]. Measurements
by Shutt, Simon, EDELWEISS, and CDMS-II are in the mK
temperature range, while all others are at 77 K. The CDMS II
yield measurements are from a single detector [11] that is simi-
lar to the averaged yield over all detectors. However, variation
across individual detectors was more than the statistical un-
certainty of the measurement shown here [12]. EDELWEISS
has reported two measurements: (i) discrete yield measure-
ments using a
252
Cf neutron source [15], and (ii) continuous
measurement using an
241
Am
9
Be neutron source fitted by a
power law function considered as their yield model [20]. The
work by Collar
et al.
reports three different approaches to
measuring ionization yield in germanium: (i) using photo-
neutrons from an
88
Y
9
Be source indicated by the red band,
(ii) using nuclear recoils from thermal neutron capture in-
dicated by the red hollow square marker, and (iii) using an
iron-filtered low energy neutron beam [21]. The prediction by
Lindhard et. al. [22, 23] is also shown for comparison.
The SuperCDMS SNOLAB experiment will search for
low mass dark matter using Ge and Si high-voltage (HV)
detectors operated at temperatures below 40 mK [2].
These detectors measure phonon signals that are propor-
tional to the bias voltage and the amount of ionization
in the detector due to a recoil event (see Sec. V A for
3
details). The cryogenic HV detectors have a low recoil
energy threshold, but cannot make a direct measurement
of the ionization yield on an event-by-event basis as the
ionization and phonon signals are measured together as
a combined phonon energy response of the detector [2].
A key goal of this study is to provide an ionization yield
measurement in Ge that is applicable to the SuperCDMS
SNOLAB detectors at low recoil energies [2].
In this work, we studied the nuclear recoil ionization
yield in Ge using the SuperCDMS iZIP detectors [27] in
HV mode [28] called CDMSlite [29], which were oper-
ated as part of the SuperCDMS Soudan experiment [30].
To generate nuclear recoils in the detectors, we exposed
them to photo-neutron sources, which produce quasi-
monoenergetic neutrons. Such sources can be fabricated
by pairing a radioactive source that has an intense, high
energy gamma with Be wafers to induce photoproduction
of neutrons from
9
Be [31, 32]. To extract the ionization
yield, we performed a Geant4 simulation to model the
recoil spectrum in the detectors. We then performed a
maximum likelihood fit using the input simulation spec-
trum and a two-parameter yield model to extract the ion-
ization yield from the data. We discuss the experimental
setup and data sets in Section II, the data selection cri-
teria and their selection efficiencies in Section III, the
simulation of the nuclear recoil spectrum in Section IV,
the yield extraction using the likelihood function in Sec-
tion V, the results of this work are discussed in Section
VI, and our concluding remarks are presented in Section
VII.
II. EXPERIMENTAL SETUP
The SuperCDMS Soudan experiment, operated from
2012 to 2015, was located at the Soudan Underground
Laboratory in northern Minnesota. The dark matter
experiment consisted of 15 high-purity germanium de-
tectors that were arranged in a compact array of five
towers with three detectors per tower, as shown in
Fig. 2 (b). Each detector was roughly cylindrical in
shape and approximately 600 g in mass. The top and
bottom faces of each detector were instrumented with
tungsten transition-edge sensors (TES) coupled to alu-
minum energy-collecting “fins” used to read-out ather-
mal phonon signals and provide electrodes for measur-
ing ionization signals. A dilution refrigerator cooled the
payload to a baseline operating temperature that varied
between 40 and 50 mK. Further details about the exper-
imental setup and the detector design can be found in
Ref. [28, 33].
To perform the photo-neutron calibration measure-
ment, the experimental configuration was modified
slightly after the main run that was used for dark matter
searches. Figure 2 (a) shows a schematic of the exper-
iment as it was configured for the photo-neutron mea-
surement. The detectors were operated in two modes:
(i) iZIP [27], and (ii) CDMSlite [34]. In the CDMSlite
mode, the detectors were biased at higher voltages (up
to 70 V) than in iZIP mode to take advantage of the
Neganov-Trofimov-Luke (NTL) effect. The NTL effect
is the production of secondary phonons in proportion to
an ionization energy deposit and an applied bias volt-
age [35, 36]. The higher voltage applied to CDMSlite
detectors thus served to amplify the phonon signal from
the primary ionization, effectively lowering the energy
threshold of the detector.
Two of the detectors, labelled as T5Z2 and T2Z1 (see
Fig. 2 (b)), were operated in CDMSlite mode. They were
biased at 70 V and 25 V, respectively. These detectors
were chosen because they showed adequate signal to noise
performance when biased at higher voltages. The photo-
neutron source was placed above the respective towers
(see Fig. 2 (a)).
Two gamma sources,
88
Y (half-life 106.6 days) and
124
Sb (half-life 60.2 days), were deployed sequentially for
the photo-neutron calibration. Each source was placed
on top of a solid
9
Be disk.
9
Be is the only stable isotope of
Be that occurs naturally in non-trace amounts. The ac-
tivity of each of the gamma sources at the start of the cal-
ibration period was 1 mCi. The most prominent gamma
above the photo-disassociation threshold (1.66 MeV) that
is emitted by
88
Y is at 1.8 MeV. In the case of
124
Sb, the
most prominent gamma above threshold is 1.69 MeV. In
the photoproduction process, the
9
Be nucleus absorbs a
gamma and emits a neutron to become a
8
Be resonant
state. The resulting neutrons emitted from
9
Be, over
multiple occurrences of this process, are nearly mono-
energetic, with an average energy of 152 keV [19] for
88
Y
and 24 keV [37] for
124
Sb. The diameter of the
9
Be disc
was 5 cm with a thickness of 0.2 cm. The purity level of
the Be disc was 98.5%
9
Be. A source holder was made to
constrain the position of the gamma source with respect
to the
9
Be disc.
There is a small spread in the energy of the emit-
ted neutrons, which depends on the angle at which they
were emitted with respect to the direction of the inci-
dent gamma [31]. The spread in neutron energies for
the
88
Y
9
Be source is
8 keV, and for the
124
Sb
9
Be
source is
1
.
3 keV. Only neutrons predominantly mov-
ing in the forward direction (towards the detector) are
likely to reach the detectors, due to the small solid an-
gle subtended by the detector. The angle-energy cor-
relation was not included in the Geant4 simulation; in-
stead we used a mono-energetic energy consistent with
the forward-direction. A cross-check was performed us-
ing the backward-direction energy, and the difference was
assigned as a systematic uncertainty.
The ratio of all gammas to one neutron produced from
the photo-neutron source is approximately 10
5
. Such a
high rate of events had the potential to induce a pro-
hibitive amount of dead-time for the data acquisition
(DAQ). To suppress the high rate of these gammas, 13
to 15 cm of lead shielding was placed between the photo-
neutron source and the detectors at different stages of
the run (see Fig. 2 (a)). The thickness of the lead shield-
4
FIG. 2. (a) The experimental setup with the shielding config-
uration. Lead was used to shield against gammas. Polyethy-
lene shielding was required for the dark matter configuration
and only the top portion was removed from the experiment in
order to perform the photo-neutron calibration. (b) Detectors
were arranged into 5 towers with 3 detectors in each tower.
Each face of the detectors had four phonon channels. The
phonon channel layout is indicated by the different colours on
the upper face of the detector. The photo-neutron source was
placed above and between towers T2 and T5. The T2Z1 and
T5Z2 detectors used in this analysis have been highlighted
in pink. The detectors T2Z2 and T5Z1, highlighted in blue,
were used for a systematic study of the event rate due to the
position of the source over the two towers.
ing between the source and the detectors was optimized
to balance degradation of the neutron energy spectrum
against the desire to adequately reduce the DAQ dead-
time. Additionally, the thickness of the lead shielding
between the photo-neutron source and the detectors was
decreased as the decay rate of the source diminished. A
TABLE I. Summary of the photo-neutron data-taking runs.
Source
Duration
Detector
Voltage
124
Sb and
124
Sb
9
Be
62 days
T5Z2
70 V
88
Y and
88
Y
9
Be
42 days
T5Z2
70 V
88
Y and
88
Y
9
Be
38 days
T2Z1
25 V
special “window” trigger was implemented at the soft-
ware level. This served to further reduce the rate of
stored events and manage the DAQ dead-time by ve-
toing high-energy electron recoils that produced energy
deposits well above the nuclear recoil energy spectrum.
The photo-neutron data-taking took place for
144 days between June 5th, 2015 and October 26th,
2015. Data were recorded with the
9
Be disk (neutron-
ON) and without the
9
Be disk (neutron-OFF). The latter
served to measure the gamma background. A summary
of the data-taking condition is shown in Table I. The data
was first taken with the T5Z2 detector, and subsequently
with the T2Z1 detector. As the activity of the
124
Sb
source (
τ
1
/
2
= 60.2 days) had become much weaker by
the time of T2Z1 operation in CDMSlite mode, no data
were taken with this source in that operation period.
III. DATA RECONSTRUCTION, SELECTION
AND EFFICIENCIES
The phonon pulse shape observed in the various chan-
nels features a steeply rising edge (few microseconds) and
a slowly falling, exponential tail (with a time constant of
a few milliseconds) [38]. The rising edge is governed by
the timescale for the earliest phonons to reach the sen-
sors and provides information about the position of the
interaction. The phonon sensor coverage is relatively low
(
5%), and thus the majority of phonons are reflected
multiple times from uninstrumented surfaces before they
are detected, contributing to the slowly falling tail. This
slow component of the signal provides information about
the total event energy and the tail shape is largely inde-
pendent of position throughout the detector.
The data readout scheme and reconstruction are de-
scribed in detail in Ref. [38]. As described there, the
energy estimation was done using a variation on an op-
timal filter (OF), referred to as the “non-stationary” op-
timal filter (NSOF) which used a template to fit the
recorded pulses. The template was made by averaging
over a selection of recorded pulses that were deemed to
be of high quality and representative of real, particle-
induced events. The NSOF treats the differences be-
tween the template and the fitted pulse as non-stationary
noise. Since variation in the rising edge of the pulse arises
due to position dependence in the phonon response, the
NSOF had the effect of deweighting the fit to position-
dependent systematics and yielded a more precise energy
estimation. In addition to fitting the recorded pulses to
a standard pulse template, the pulses were also fit with
templates that matched “glitches” and typical low fre-
5
quency noise (LFN). The glitches are electronic in origin
and are characterised by pulses with an unusually fast
rise and fall time. The LFN (
<
1 kHz) originated from
vibrations created by a cryocooler, which provided sup-
plemental cooling power to the experiment. A correlation
between the cryocooler activity and the LFN in CDM-
Slite detectors was observed and reported in detail in
Ref. [38]. The
χ
2
goodness of fit from the NSOF was
stored for each fit and used for the selection of events as
described in more detail below.
Prior to determining the ionization yield for nuclear re-
coils, it was necessary to calibrate the response of the de-
tectors to electron recoils. To do this, the germanium de-
tectors were exposed to neutrons from
252
Cf. This led to
a neutron capture process, creating unstable
71
Ge.
71
Ge
decays via electron capture to
71
Ga, releasing predom-
inantly x-rays with a total energy corresponding to the
binding energy of K-, L- or M-shell electrons. The energy
scale of each detector was calibrated using the 10.3 keV
and 1.3 keV peaks from the K- and L-shell electron cap-
ture lines in
71
Ge [38–40]. The calibration values from a
prior dark matter search were adopted for this work, fol-
lowing the above described technique, and are described
in further detail in Ref. [38]. For this dataset, the sta-
bility of the adopted calibration values was then checked
using the electron capture lines as observed in dedicated
low background datasets (i.e. no source placed during
the run) taken at periodic intervals before, after, and be-
tween the photo-neutron runs. These checks showed that
the energy calibration was stable to within 2% through-
out the several-months period when photo-neutron data
was collected.
To optimize data quality, several event selection cri-
teria (cuts) were applied in the analysis. A “livetime”
cut was used to remove specific time periods of data
from the analysis. “Quality cuts” were applied on an
event-by-event basis using the recorded pulse shapes. A
“threshold” cut was used to set the low-energy analysis
threshold.
Livetime cuts:
Cuts were applied to ensure data
were taken with the correct running conditions. These
included checks on the base temperature, detector bias
voltage, and time between LED flashes (near-IR LEDs
were used to release trapped charges from the detector)
[41]. Additionally
1% of livetime was removed during
periods of unusually high event rates, based on the aver-
age trigger rate calculated every 30 s of runtime. A few
runs (totaling less than 12 hours) that showed an excep-
tionally low noise environment were also removed from
the dataset. The “pre-trigger” baseline for each phonon
pulse was recorded for approximately 1 ms before each
event. Excessive fluctuation or large RMS in the pre-
trigger baseline is indicative of noise or event pileup. A
cut was applied to remove events for which the standard
deviation in baseline of the phonon pre-trigger region de-
viated by more than 4
σ
from its mean value, as calculated
over the course of the run. The livetime cuts removed
10% of the events from the 70 V datasets, and
5% of
the events from the 25 V dataset.
Quality cuts:
As described above, the energy of each
event was reconstructed from the phonon pulses using a
variant on an optimal filter routine called NSOF. Using
the NSOF template fit to phonon pulses, a
χ
2
goodness
of fit was calculated for each event. The
χ
2
as a function
of fit energy was constant over the neutron scattering
energy range, indicating that the fit was well-behaved
over the energy range of interest. Events with very large
χ
2
or with a rise time outside a 230
μ
s window around
the trigger time were removed from the analysis. Such
events typically consisted of pileup events or other poorly
reconstructed events.
In the CDMSlite detectors, the signals from the ioniza-
tion channels cannot be used in the reconstruction of the
photo-neutron recoil spectrum, as they do not recover the
entire charge energy from an ER or NR event. However,
positive correlations were found between the charge and
phonon goodness of fit for extremely high charge OF
χ
2
values (
χ
2
>
6250 for data from T2Z1, and
χ
2
>
50000 for
data from T5Z2). Thus, events with a high charge
χ
2
,
which typically displayed a sharp glitch in the waveform,
were rejected from the analysis to remove bad events that
survived phonon-based selection criteria.
In addition to the standard pulse template, each
phonon pulse was also fit with both an LFN and a glitch
template. Based on the goodness of the respective fits,
a variable discriminating between signal-like and noise-
like events called ∆
OFχ
2
noise
=
OFχ
2
pulse
OFχ
2
noise
was
constructed, where the subscript “noise” stands for either
LFN or glitch. A cut at an appropriate level of signal
acceptance (3
σ
) was applied on the ∆
OFχ
2
noise
distri-
bution as a function of total phonon energy to remove
both classes of noise events. To best optimize the cuts,
the photo-neutron dataset was divided into eleven time
blocks such that the behavior of the LFN associated with
the cryocooler within each time block was similar. The
cuts to reject LFN and glitches were defined for each of
the time blocks separately.
In the data acquisition system, charge triggers were
recorded when the charge pulse exceeded a discriminator
threshold. Only phonon triggers were used to make a de-
cision on recording the data, however the charge triggers
were recorded and used for data quality studies. A good
pulse would ideally issue a trigger in both the charge and
phonon channels, but glitches on a particular channel
could create a difference in the number of triggers issued
between the channels. Thus, glitch events were further
rejected by requiring that selected events had coincident
triggers from both the phonon and charge channels.
Threshold cut:
For a given detector, the phonon
trigger efficiency calculation relied on events triggered
by the two other detectors in the same tower. Only
data from the triggering detector were read out in the
photo-neutron data-taking mode, hence a separate set of
data had to be utilized to measure the trigger efficiency.
These data required that an entire tower of detectors be
read out whenever a single detector passed the phonon
6
5
10 15 20 25 30 35 40 45 50
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Efficiency
Be (70 V)
9
Sb+
124
Be (70 V)
9
Y+
88
Be (25 V)
9
Y+
88
(a)
5
10 15 20 25 30 35 40 45 50
Total Phonon Energy (keV)
0.86
0.88
0.9
0.92
0.94
0.96
0.98
1
Efficiency
Sb at (70 V)
124
Y (70 V)
88
Y (25 V)
88
(b)
FIG. 3. Combined and smoothed cut efficiency as a function
of energy above the analysis threshold for each data set with
a 1
σ
uncertainty. (a) With the Be wafer i.e. neutron-ON, and
(b) without the Be wafer i.e. neutron-OFF.
trigger condition. This general technique has been used
by SuperCDMS in the past and is described in refer-
ences [38] [42]. The trigger efficiency for this analysis is
based on a combination of data taken with a
252
Cf source
and data taken in the low background mode. These
datasets were selected because they had trigger thresh-
olds matching those in the photo-neutron runs. They
also provided data with a continuous neutron energy dis-
tribution up to a few MeV and real physical events that
induced multiple triggering events in the full tower of
detectors. To minimize the potential impact of the sys-
tematic uncertainty in the trigger efficiency calculation
on our results, the analysis energy thresholds were set
to just above where the trigger efficiency was uniformly
100%. Lastly, in order to interpret the threshold values
in terms of recoil energy we applied a correction for the
NTL phonons. For the case of electron recoils (yield
1),
the threshold for this analysis was 74 eV for T5Z2 and
236 eV for T2Z1.
We calculated the efficiencies of the quality-cuts de-
scribed previously by relying on a data driven approach
following the procedures described in Ref. [33]. To com-
pute the efficiencies, we generated artificial events by
adding a scaled pulse template to random noise traces
gathered during the data acquisition. To take position
dependencies into account the pulse templates were gen-
erated as a linear combination of two sub-templates mod-
eling the fast and slow components of the absorption sig-
nal described above and in more detail in Ref. [38]. The
relative weights of the sub-templates were drawn from the
fast-to-slow amplitude ratio distribution observed in the
photo-neutron data. These artificial data, generated for
each data set in a time block, reproduced similar distri-
butions as observed for the measured total phonon energy
and the phonon OF
χ
2
in the real data. The efficiency
for each time block was calculated as the ratio of the
events that pass the cuts over the total number of gen-
erated events, as a function of energy. For the final cal-
culation, the efficiency functions of the constituent time
periods were combined together, weighting each period
by both its livetime and the decreasing activity of the
source (see Sec. V). The efficiency functions were then
smoothed with a Gaussian filter to reduce the effects of
statistical fluctuations. Figure 3 shows the resulting com-
bined and smoothed cut efficiency functions for each data
set. As seen in the figure, the quality cuts were generally
above 90% efficient and showed relatively little energy
dependence.
Between data sets, where the source configuration was
changed, we observed discrete jumps down in rate by
as much as 5%, which is not compatible with statistical
fluctuations in the rate or known changes in operating
conditions. The same jumps were present in the data
when no selection cuts were applied. We also measured a
smooth downward trend in rate that was consistent with
the exponential decay of the sources. A systematic study
on the effect of a slight deviation in the source position
on the overall event rate revealed that this effect was
negligible and not consistent with the rate jumps.
We used a
χ
2
test to determine whether the jumps in
rate affected the shape of the neutron-ON energy dis-
tribution. To test the consistency of the energy spectra
over time, a
χ
2
was calculated between pairs of data sets
drawn from the eleven time periods and the same photo-
neutron source. The
χ
2
tests returned a p-value which
informed whether the data were drawn from the same
parent distribution. Out of the 62 tests we performed,
only one returned a p-value lower than 0.01. Based on
this result, we concluded that the shape of the energy
spectrum is consistent over time and the fit should not
be affected by the rate fluctuations seen between data
sets.
7
IV. SIMULATIONS
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Nuclear recoil energy (keV)
1
10
2
10
3
10
4
10
5
10
Counts/0.05 keV
Neutron induced ER
Multiple scatter NR
Single scatter NR
All deposits
(a)
Be
9
Sb
124
0
2
4
6
8
10
12
14
Nuclear recoil energy (keV)
1
10
2
10
3
10
4
10
5
10
Counts/0.1 keV
Neutron induced ER
Multiple scatter NR
Single scatter NR
All deposits
(b)
Be
9
Y
88
FIG. 4.
The simulated recoil energy distributions from (a)
124
Sb
9
Be, and (b)
88
Y
9
Be in T5Z2. The energy depositions
from electron recoils (ER) that arise from secondary gammas
is shown in blue, the energy depositions from multiple-scatter
nuclear recoils (NR) by neutrons in the detector is shown
in green. The energy distribution of single-scatter NR by
neutrons is shown in red.
A primary input to the likelihood calculation is the
neutron energy probability distribution function (PDF).
Geant4 10.6 [43] was used to simulate the neutron energy
spectrum from the different source configurations. The
Geant4 NeutronHP physics model [44] along with the
G4NDL4.6 cross-section packages were implemented for
the simulation of 1.2
×
10
9
neutrons propagating through
the experimental geometry. G4NDL4.6 is the JEFF3.3
library. It produces results closer to MCNP, has more
isotopes than ENDF/B, and has an overall lower error in
the energies of secondary neutrons [45, 46]. The Compton
spectrum from the source gammas was modeled with an
analytical spectrum and cross-checked with the neutron-
OFF data, as described in detail in Sec. V B. Therefore,
the only simulated data used in this analysis was the neu-
tron sample, which included a very small contribution of
secondary gammas originating from neutron interactions
in the surrounding materials.
Figure 4 shows the recoil energy deposition in T5Z2
of neutrons from the
124
Sb
9
Be and
88
Y
9
Be sources in
panels (a) and (b) respectively. Single and multiple scat-
ter nuclear recoil (NR) energy deposits from neutrons
interacting directly in the detector and electron recoils
(ER) from neutron-induced gammas are shown sepa-
rately. The sharp drop-off in the single-scatter neutron
spectrum (“endpoint”) corresponds to the maximum en-
ergy a neutron of the given input energy can transfer to
a Ge nucleus in a single interaction (
1.3 keV for the 24
keV neutrons from
124
Sb
9
Be and
8.2 keV for the 152
keV neutrons from
88
Y
9
Be; see Fig. 4).
There are a large number of multiple neutron scatters,
which obscure the endpoint of the single scatter neutrons
in the total energy distribution. This makes it difficult to
determine the exact position of the maximum recoil for
single-scatters in the actual data distribution, and thus to
directly determine the ionization yield with an integral-
based method as described in Ref. [37]. Additionally, the
neutrons pass through layers of lead, polyethylene, and
copper before reaching the germanium detectors. This
can produce neutron induced gammas from these mate-
rials giving rise to an ER background. As seen in Fig. 4,
the rate of the ER background is low compared to the
neutron scattering rate and is almost flat over the en-
tire energy range. Furthermore it is subdominant to the
Compton scattering spectrum that results from the gam-
mas directly produced by the source (as described in the
next section). The various peaks seen in Fig. 4 (b) arise
from the electron capture processes in germanium.
We studied two systematic uncertainties in the neutron
PDF arising from the simulation. The first originates
from the choice of physics model used in Geant4. To
study this uncertainty, the simulation data were regener-
ated with the same sample size as used for the primary
analysis with the LEND [47] package in place of the Neu-
tronHP [44] and G4NDL4.6 packages [48]. The likelihood
analysis was re-done using this alternative simulation
dataset (see Sec. V C 1). The second uncertainty is in-
troduced by the neutron-nucleus scattering cross-section
information available for germanium through the Geant4
package. Figure 5 shows the cross-section in germanium
as a function of neutron energy for three different ex-
ample realizations of germanium cross-section files from
8
the TENDL-2017 nuclear data library [49]. The TENDL
data library is created via the Total Monte Carlo method
[50], wherein input parameters to the TALYS nuclear
reaction simulation are randomly varied around known
central values, which are in turn determined from a com-
bination of experimental data and computational mod-
els. The resulting cross-sections show some differences,
particularly in the neutron energy region below 50 keV.
Using NeutronHP, we re-ran the simulation with 100 dif-
ferent realizations of the neutron cross-sections. Each
simulation was for 6
×
10
8
neutrons. These simulation
datasets were used for determining the systematic uncer-
tainty introduced by this effect as described further in
Sec. V C 1. The two systematics are comparable in size,
with the one from using the LEND library instead of the
NeutronHP and G4NDL4.6 being slightly larger.
1
10
2
10
3
10
Neutron Energy (keV)
1
10
1
10
2
10
3
10
4
10
Cross section (barn)
realization 1
realization 2
realization 3
10
20
1
10
2
10
Zoomed
FIG. 5.
Three examples of realizations of neutron-nucleus
scattering cross-section over a large neutron energy range
from the TENDL-2017 nuclear data library [49]. Inset is
zoomed into the neutron energy range of interest where the
differences between the cross sections are considerable be-
tween the three realizations. See text body for description
of how TENDL cross-section realizations are generated.
V. YIELD EXTRACTION
In order to extract the nuclear recoil ionization yield,
we parameterized the conversion from total phonon en-
ergy to nuclear recoil energy following a Lindhard-like
model where the parameters were determined by fitting
the model to the neutron-ON data. The fit minimized
the summed negative log likelihood of all events given a
PDF, which is the sum of a PDF of the gamma back-
ground model and a PDF for the neutron component.
The neutron component is constructed from the values
of the yield modeling parameters, the simulated neutron
spectrum, and a model of the detector response. The
signal modeling is described in Sec. V A, and the back-
ground model in Sec. V B. Section V C discusses the de-
tails of the likelihood fit, the propagation of statistical
and systematic uncertainties, and the final results.
A. Signal modeling
The recoil energy (
E
r
) distributions from the two
photo-neutron sources were derived from the simulated
neutron data described in Sec. IV. In order to compare
the simulation to real data, the simulated energy depo-
sition was converted into the total phonon energy (
E
p
),
the observable measured by the detectors.
The neutron signal is modeled for the three source con-
figurations: (i)
124
Sb
9
Be with the detector biased at 70 V,
(ii)
88
Y
9
Be with the detector biased at 70 V, and (iii)
88
Y
9
Be with the detector biased at 25 V. When a bias
voltage is applied to the detectors, the phonon energy
generated through the NTL effect must be taken into ac-
count. Work done by the electric field (
E
NTL
) in drifting
the
e
/h
+
pairs to the charge electrodes is converted into
NTL phonons which are detected by the TES sensors
along with the phonons generated by the primary recoil
interaction. The total phonon energy is the sum of the
recoil and NTL energies. The NTL energy is calculated
as
E
NTL
=
N
e/h
e
V
=
Y
(
E
r
)
E
r

γ
e
V,
(1)
where
N
e/h
is the number of
e
/h
+
pairs created by the
recoil, ∆
V
is the bias voltage, and
e
is the electronic
charge.
Y
(
E
r
) is the ionization yield, defined as the ra-
tio between the ionization signal and the primary recoil
energy, and

γ
= 3 eV [51–53] is the average energy re-
quired per generated
e
/h
+
pair in an ER interaction in
germanium.
The total phonon energy can now be written as,
E
p
=
E
r
(1 +
Y
(
E
r
)
e
V/
γ
)
,
(2)
where the dependence of
E
p
on the parameter of interest,
Y
(
E
r
), appears clearly.
The Lindhard theory [22, 54–56], supported by earlier
measurements [23], provides a semi-empirical prediction
for
Y
(
E
r
) of NR interactions for a material of mass num-
ber
A
and atomic number
Z
:
Y
(
E
r
) =
k
g
(

)
1 +
kg
(

)
,
(3)
where:
g
(

) = 3

0
.
15
+ 0
.
7

0
.
6
+

;

= 11
.
5
E
r
Z
7
/
3
.
(4)
The parameter
k
describes the electronic energy loss and
nominally ranges from 0.156 to 0.160 for stable isotopes
of Ge.
9
TABLE II. Values of the parameters derived from the en-
ergy resolution model. The energy scale is electron-equivalent
(eV
ee
), where all events are assumed to be electron-recoils.
Bias voltage
σ
B
(eV
ee
)
F
A
70 V
8.31
±
0.08
0.19
±
0.07
0.0095
±
0.0015
25 V
18.71
±
0.16
0.27
±
0.02
0.0107
±
0.0005
We introduce a modification to the Lindhard model in
which we allow the
k
value to vary linearly with the recoil
energy:
k
(
E
r
) =
k
low
+
k
high
k
low
E
high
E
low
(
E
r
E
low
)
.
(5)
E
low
is equal to 0.39 keV and
E
high
is equal to 7.0 keV.
These are the minimum and maximum nuclear recoil en-
ergy that the fit was sensitive to. In the following, we
refer to
k
low
and
k
high
as the two components of a vector
~
k
.
In order to compare the simulated neutron signal to
the collected data, we included effects due to the detec-
tor resolution to the simulated neutron spectrum. The
resolution model (
σ
T
) used was:
σ
2
T
(
E
) =
σ
2
B
+
σ
2
F
(
E
) +
σ
2
D
(
E
)
.
(6)
where
σ
B
is the baseline noise due to the electronics,
σ
F
is the variance in the number of
e
/h
+
pairs produced in
a recoil event, and
σ
D
is an empirical term that models
detector effects and is given by
σ
2
D
=
AE
2
.
The values of
σ
B
are obtained from the energy distri-
bution of events randomly triggered in the detector and
are listed in Tab. II. The second term can be written as
σ
2
F
=
FN
e/h

2
γ
=
F
γ
E
, where
F
is the Fano factor. The
resolution model can thus be rewritten as:
σ
T
(
E
) =
σ
2
B
+
F
γ
E
+
AE
2
,
(7)
To determine the values of
F
and
A
in Eq. 7, the res-
olutions of each of the K-, L- and M-shell peaks from
electron-captures (EC) decays of
71
Ge were determined
by first fitting them with a Gaussian and obtaining their
1
σ
values. Next, the model described in Eq. 7 was fit to
the resolutions of these EC peaks. The values of
F
and
A
were determined from this fit. Table II shows the values
of the different parameters determined for the resolution
model at the operating voltages of 70 V and 25 V. Since
the Fano factor for nuclear recoils vs electron recoils can
be significantly different, a systematic uncertainty on F
is discussed in V C 1.
Spatial variations in the uniformity of the electric field
caused a variation in the amount of NTL gain achieved.
Near the edge of the detector, the electric field is inho-
mogeneous, leading to a reduction in the observed NTL
phonon contribution. To account for this a 3-D electric
field map was calculated and taken into account when
determining the total energy of the simulated events. De-
tails of this technique are described in Ref. [38].
Finally, to determine the neutron PDFs used in fitting
the data, the resolution-smeared energies were converted
to the total phonon energy using the parametrized model
of Eq. 5. In order to mitigate the statistical fluctua-
tion effects, we smoothed the energy distributions using
a Savitzky-Golay filter [57] before interpolating them us-
ing a cubic spline to determine the PDF.
B. Modeling electron recoil backgrounds
The dominant background that affected this measure-
ment was photons from the source Compton-scattering
off the electrons in the Ge crystal. Another source of
background was the K-, L- and M-shell electron capture
x-rays.
The Compton continuum is characterized by steps in
the scattering rate (i.e., x-ray absorption edges) [58], cor-
responding to each of the Ge electron shells. The Comp-
ton model in the region of interest is:
f
C
(
E
) =
a
0
[1 +
S
M
Θ(
E,μ
M
)]
+
S
L
a
0
[1 +
S
M
Θ(
E,μ
M
)]Θ(
E,μ
L
)
+
S
K
a
0
[1 +
S
M
Θ(
E,μ
M
)]Θ(
E,μ
L
)Θ(
E,μ
K
);
(8)
here
S
i
is the fractional step amplitude where
i
can be
the K-, L-, or M-shell, and
a
0
is the flat background level.
The step sizes and background level were free parameters
that were determined by fitting the neutron-OFF data.
The parameter
μ
is the central energy of the step, which
was determined via a
Geant4
simulation. Θ is the error
function:
Θ(
E,μ
) =
1
2
[
1 + erf
(
E
μ
2
σ
(
μ
)
)]
,
(9)
with the detector resolution
σ
(
μ
).
The electron capture peaks were modeled using Gaus-
sian functions. Following Ref. [59], we assumed the frac-
tional electron capture probabilities for each shell to be
f
K
= 87
.
6%,
f
L
= 10
.
5%, and
f
M
= 1
.
8% for the K-, L-,
and M-shells, respectively. The only parameter deter-
mined with a fit to the neutron-OFF data for this back-
ground was the amplitude of the K-shell peak. We then
constrained the other amplitudes using the fractional ra-
tio. While data collected at 70 V are sensitive to the L-,
K-, and M-shells, for the 25 V data the M-shell step oc-
curs below the analysis threshold, and was thus excluded
from the model of the two backgrounds.
The combined background model from the EC peaks
and Compton steps was fit to the neutron-OFF energy
distribution. In order to include any outstanding effect
which is not modeled, we added the fit residual, smoothed
by applying a Gaussian filter, to the model. Fig. 6 shows
the neutron-OFF energy spectrum with the background
model overlaid. The Compton steps are clearly visible
and described by the model for the Yttrium at 70 and
25 V. A systematic uncertainty for the background model
is presented in Sec. V C 1.
10
0 2 4 6 8 10 12 14 16 18 20 22
Total Phonon Energy (keV)
0
50
100
150
200
250
300
350
400
Counts/0.3 keV
Antimony (70 V)
Compton steps
10)
×
EC peaks (
Compton steps+EC peaks+residual
data
Analysis threshold
(a)
0
50
100
150
200
250
Total Phonon Energy (keV)
0
100
200
300
400
500
600
700
Counts/2.75 keV
Yttrium (70 V)
Compton steps
10)
×
EC peaks (
Compton steps+EC peaks+residual
data
Analysis threshold
(b)
0
20
40
60
80
100 120 140
Total Phonon Energy (keV)
0
0.2
0.4
0.6
0.8
1
1.2
3
10
×
Counts/1.5 keV
Yttrium (25 V)
Compton steps
10)
×
EC peaks (
Compton steps+EC peaks+residual
data
Analysis threshold
(c)
FIG. 6. Neutron-OFF energy spectrum for the three data sets: (a)
124
Sb
9
Be with T5Z2 biased at 70 V, (b)
88
Y
9
Be with T5Z2
biased at 70 V and (c)
88
Y
9
Be with T2Z1 biased at 25 V. The contribution from each component in the background model is
shown. The contribution of the electron capture peaks have been scaled up by a factor of ten for better visualization.
Finally, the energy-dependent cut efficiencies for the
quality cuts described in Sec. III were applied to the
signal and background PDF in order to represent the
expected shape of the experimental data with all cuts
applied. For the background PDF, since it was obtained
by exploiting the neutron-OFF data, we first applied the
inverse of the neutron-OFF cut efficiency function to re-
move any effects on its shape due to the neutron-OFF
cut efficiency.
C. Likelihood analysis
The background and signal PDFs allowed us to define
the negative log likelihood function as:
ln
L
=
3
D
=1
N
D
i
=1
ln(
f
D
ν
D
(
E
i
,
~
k
) + (1
f
D
)
b
D
(
E
i
))
,
(10)
where
N
D
is the number of events in the data set
D
,
f
D
is the fractional contribution of the neutron signal,
ν
D
(
E,
~
k
) are the parameter-dependent signal PDFs, and
b
D
(
E
) are the background PDFs. The free parameters
of the negative log likelihood function, which were min-
imized using the MINUIT [60] package, were the three
neutron contribution fractions
f
D
and the Lindhard pa-
rameters
~
k
. The relative goodness of fit, which takes into
account only the statistical uncertainty, was evaluated by
binning the experimental data and calculating the fit
χ
2
per degree of freedom for each data set. Figure 7 shows
the energy spectra with the best fit result, for the 2-
parameter model, overlaid on the three data sets. The
χ
2
/
(degrees of freedom)
was found to be 137
/
95, 141
/
95
and 129
/
95 for
124
Sb at 70 V,
88
Y at 70 V and
88
Y at
25 V respectively. The best fit values of
~
k
are shown in
Tab. III.
In order to study whether the 1-parameter (standard
Lindhard) or 2-parameter (see Eq. 5)
k
model better de-
scribed the data, we performed a likelihood ratio test. We
generated 3000 Monte Carlo (MC) data sets sampling a
number of events equal to the experimental data, from
the background model plus the simulated neutron spec-
trum, which was converted to total phonon energy using
the 1-parameter Lindhard, see Eq. 3. The ratio between
the number of background and neutron events and the
value of the
k
parameter in the MC data generation were
set to the best-fit values. We fit the MC data both using
the 1-parameter and the 2-parameter modified Lindhard
model. The likelihoods obtained in the two fits were used
to construct the likelihood ratio, which was used as a test
statistic to determine which model was a better fit to the
data. Since the two models were nested, we used the
likelihood ratio test to compare the goodness of fits in
the two cases. Given that the 1-parameter model was
used to generate the MC data, the distribution of the
likelihood ratio obtained with the MC data samples pro-
vides a quantitative measurement of how much better
the 2-parameter model fit performed compared to the 1-
parameter model, based on only statistical fluctuations.
Finally, we calculated the likelihood ratio using the
measured data and we compared it to the distribution
obtained using the MC data samples. We found no single
likelihood ratio value calculated using the simulated data
that had a value greater than the collected data. We thus
inferred that the 2-parameter hypothesis, which allowed
for a linear energy dependence of the
k
parameter of the
Lindhard model, was preferred with a significance greater
than 3
σ
. Based on this result, we report only the result
for the 2-parameter modified Lindhard model.
1. Calculation of uncertainties
We took into account the statistical uncertainties due
to the finite size of the simulated neutron spectrum and of
11
2 4
6 8 10 12 14 16 18 20 22
Total Phonon Energy (keV)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
3
10
×
Counts/0.2 keV
Be (70 V) data
9
Sb+
124
Best fit function
Background model
Analysis threshold
(a)
0
50
100
150
200
250
Total Phonon Energy (keV)
1
1.5
2
2.5
3
3.5
4
3
10
×
Counts/2.4 keV
Be (70 V) data
9
Y+
88
Best fit function
Background model
Analysis threshold
(b)
0
20
40
60
80 100 120 140
Total Phonon Energy (keV)
8
10
12
14
16
3
10
×
Counts/1.4 keV
Be (25 V) data
9
Y+
88
Best fit function
Background model
Analysis threshold
(c)
FIG. 7.
Energy spectrum for the three data sets, (a)
124
Sb
9
Be with T5Z2 biased at 70 V, (b)
88
Y
9
Be with T5Z2
biased at 70 V and (c)
88
Y
9
Be with T2Z1 biased at 25 V. The
best fit result for each dataset is overlaid in red. The blue
curve shows the fit gamma background contribution.
TABLE III. Summary of the
k
low
and
k
high
fit results along
with their statistical and systematic uncertainties
Best fit value
Stat. uncertainty
Sys. uncertainty
k
low
0.040
0.005
0.008
k
high
0.142
0.011
0.026
the experimental neutron-ON and neutron-OFF spectra
by performing fits to simulated experiments. The dis-
tributions of the fit results were then fit with Gaussian
distributions, and the standard deviation of each Gaus-
sian was propagated as a statistical uncertainty for the
respective contribution. We propagated the uncertainty
on the cut efficiency functions by repeating the analysis
twice, using the efficiencies shifted by one standard devi-
ation up and down with respect to their central values.
We took the difference in the results as estimate for the
uncertainty in the yield.
The value of the Fano factor used in the energy reso-
lution model was determined from a fit to data (predom-
inantly electron recoils) and it was affected by statistical
and systematic uncertainties. In particular, a previous
measurement [61] showed that the Fano factor for nu-
clear recoils could be significantly higher than in elec-
tron recoils. To account for these uncertainties, the fit
was repeated twice, once forcing a downward shift for
F
compatible with its statistical uncertainty, and once forc-
ing an upward shift to a value of 10 for nuclear recoils,
which is large enough to include any measurements in
literature [62, 63]. The resulting shift in the fit results
were taken as the estimate for systematic uncertainty as-
sociated with the Fano factor. In order to evaluate the
uncertainty on the background model shape, described
in Sec. V B, we recalculated the yield by using as back-
ground PDF the analytical model only, without adding
the residuals from the fit to the neutron-OFF data. The
difference with respect to the result obtained using the
standard background model was used as an estimate for
this systematic uncertainty.
We studied the systematic uncertainty arising from
the limited knowledge of the neutron elastic scattering
cross section used to generate the signal model. The
neutron simulation was repeated using several different
realizations of cross section files described in Sec. IV. The
negative log likelihood fit was repeated for each of the
simulations, and the standard deviations of the resulting
distribution of the fit results were used as a systematic
uncertainty. Another source of systematic uncertainty
arose from the choice of the specific neutron cross sec-
tion libraries, G4NDL4.6 [48]. In order to quantify this
uncertainty, another simulation was generated using the
LEND [47] neutron cross section library. The difference
between the fit results obtained using these two cross sec-
tion libraries was used as an estimate of this systematic
uncertainty. This approach does not take into account
possible correlated biases between the two packages.