of 14
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ño,
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
É. 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 and Astronomy, University of British Columbia,
Vancouver, British Columbia V6T 1Z1, Canada
2
TRIUMF, Vancouver, British Columbia V6T 2A3, Canada
3
Department of Physics, University of Toronto, Toronto, Ontario M5S 1A7, Canada
4
Department of Physics, Durham University, Durham DH1 3LE, United Kingdom
5
Division of Physics, Mathematics, and Astronomy, California Institute of Technology,
Pasadena, California 91125, USA
6
Department of Physics, Northeastern University, 360 Huntington Avenue,
Boston, Massachusetts 02115, USA
7
Pacific Northwest National Laboratory, Richland, Washington, D.C. 99352, USA
8
Department of Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, Texas 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, Florida 32611, USA
12
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
13
Department of Physics, South Dakota School of Mines and Technology,
Rapid City, South Dakota 57701, USA
14
SLAC National Accelerator Laboratory/Kavli Institute for Particle Astrophysics and Cosmology,
Menlo Park, California 94025, USA
15
Department of Physics, Stanford University, Stanford, California 94305, USA
16
Department of Physics, Southern Methodist University, Dallas, Texas 75275, USA
17
Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049 Madrid, Spain
18
Department of Physics, University of California, Berkeley, California 94720, USA
19
Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208-3112, USA
20
Department of Physics, Queen
s University, Kingston, Ontario K7L 3N6, Canada
21
School of Physics and Astronomy, University of Minnesota, Minneapolis, Minnesota 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, Ontario 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, Colorado 80217, USA
26
Department of Electrical Engineering, University of Colorado Denver, Denver, Colorado 80217, USA
27
Institute for Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT), 76344 Eggenstein-
Leopoldshafen, Germany
28
Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany
PHYSICAL REVIEW D
105,
122002 (2022)
2470-0010
=
2022
=
105(12)
=
122002(14)
122002-1
© 2022 American Physical Society
29
Department of Physics, University of South Dakota, Vermillion, South Dakota 57069, USA
30
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
31
Department of Physics, Santa Clara University, Santa Clara, California 95053, USA
(Received 15 February 2022; accepted 20 May 2022; published 17 June 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 and 7 keV at a temperature of
50
mK. We use a
GEANT
4 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 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.
DOI:
10.1103/PhysRevD.105.122002
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 sensitive to the
recoiling nucleus, and the dark matter interaction is inferred
through detection of the nuclear recoil
[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 information 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 ratio of the ionization energy to the
total recoil energy is commonly referred to as the
ioniza-
tion yield.
For detectors that rely on deriving the nuclear
recoil energy from the ionization energy, a precise under-
standing 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 recoil signal induced by coherent elastic neutrino
nucleus scattering
[3
5]
.
Figure
1
shows the results from a number of experiments
that have measured the ionization yield in Ge as a function
of nuclear recoil energy. The figure also provides a com-
parison between the experimental results and a theoretical
model of the ionization yield in Ge, which is derived by
Lindhard
et al.
[6]
and which is the reference model in
the field of direct dark matter searches. Experiments often
use a value of the
k
parameter (explained in Sec.
VA
),
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
[8,11
25]
. 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
[15]
that is similar
to the averaged yield over all detectors. However, variation across
individual detectors was more than the statistical uncertainty of
the measurement shown here
[16]
. EDELWEISS has reported
two measurements: (i) discrete yield measurements using a
252
Cf
neutron source
[19]
, and (ii) continuous measurement using an
241
Am
9
Be neutron source fitted by a power law function con-
sidered as their yield model
[24]
. The work by Collar
et al.
reports
three different approaches to measuring ionization yield in
germanium: (i) using photoneutrons from an
88
Y
9
Be source
indicated by the red band, (ii) using nuclear recoils from thermal
neutron capture indicated by the red hollow square marker, and
(iii) using an iron-filtered low energy neutron beam
[25]
. The
prediction by Lindhard
et al.
[6,26]
is also shown for comparison.
*
Corresponding author.
vijayiyer@niser.ac.in
M. F. ALBAKRY
et al.
PHYS. REV. D
105,
122002 (2022)
122002-2
whose physical origin is related to the electronic stopping
power
[7]
, to better match the Lindhard model to their
experimental results. The Lindhard model, however, does
not account for some factors that may be important at
low energy, such as atomic binding
[8]
, electric field
dependence
[9]
, and temperature dependence
[10]
. 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 disagreement in the ionization yield
measurements by various experiments in this recoil energy
regime.
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 proportional to
the bias voltage and the amount of ionization in the detector
due to a recoil event (see Sec.
VA
for 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 operated
as part of the SuperCDMS Soudan experiment
[30]
.To
generate nuclear recoils in the detectors, we exposed them
to photo-neutron sources, which produce quasimonoener-
getic 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
GEANT
4
simulation to model the recoil spectrum in the
detectors. We then performed a maximum likelihood fit
using the input simulation spectrum and a two-parameter
yield model to extract the ionization yield from the data. We
discuss the experimental setup and data sets in Sec.
II
,
the data selection criteria and their selection efficiencies in
Sec.
III
, the simulation of the nuclear recoil spectrum in
Sec.
IV
, the yield extraction using the likelihood function
in Sec.
V
, the results of this work are discussed in Sec.
VI
,
and our concluding remarks are presented in Sec.
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 experi-
ment consisted of 15 high-purity germanium detectors 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 aluminum energy-collecting
fins
used to read
out athermal phonon signals and provide electrodes for
measuring ionization signals. A dilution refrigerator cooled
the payload to a baseline operating temperature that varied
FIG. 2. (a) The experimental setup with the shielding configu-
ration. Lead was used to shield against gammas. Polyethylene
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 ar-
ranged into five towers with three 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.
IONIZATION YIELD MEASUREMENT IN A GERMANIUM
...
PHYS. REV. D
105,
122002 (2022)
122002-3
between 40 and 50 mK. Further details about the exper-
imental setup and the detector design can be found in
Refs.
[28,33]
.
To perform the photo-neutron calibration measurement,
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 experiment as it
was configured for the photo-neutron measurement. 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 secon-
dary phonons in proportion to an ionization energy deposit
and an applied bias voltage
[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, labeled as T5Z2 and T2Z1 [see
Fig.
2(b)
], were operated in CDMSlite mode. They were
biased at 70 Vand 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 nontrace amounts. The activity of each
of the gamma sources at the start of the calibration 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 photoproduc-
tion 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 monoenergetic, with an average
energy of 152 keV
[23]
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 emitted
neutrons, which depends on the angle at which they were
emitted with respect to the direction of the incident 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 moving in the forward direction
(towards the detector) are likely to reach the detectors, due
to the small solid angle subtended by the detector. The
angle-energy correlation was not included in the
GEANT
4
simulation; instead, we used a monoenergetic energy
consistent with the forward direction. A cross-check was
performed using 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
prohibitive 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 shielding
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 special
window
trigger was implemented at the software level.
This served to further reduce the rate of stored events and
manage the DAQ dead time by vetoing 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 con-
dition is shown in Table
I
. The data were 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 sensors 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
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
M. F. ALBAKRY
et al.
PHYS. REV. D
105,
122002 (2022)
122002-4
energy and the tail shape is largely independent of position
throughout the detector.
The data readout scheme and reconstruction are
described in detail in Ref.
[38]
. As described there, the
energy estimation was done using a variation on an optimal
filter (OF), referred to as the
nonstationary
optimal filter
(NSOF) that 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 between the template and the
fitted pulse as nonstationary noise. Since variation in
the rising edge of the pulse arises due to position depend-
ence 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 frequency noise (LFN). The glitches are
electronic in origin and are characterized by pulses with
an unusually fast rise and fall time. The LFN (
<
1
kHz)
originated from vibrations created by a cryocooler, which
provided supplemental cooling power to the experiment. A
correlation between the cryocooler activity and the LFN in
CDMSlite 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
recoils, it was necessary to calibrate the response of the
detectors to electron recoils. To do this, the germanium
detectors 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 predomi-
nantly 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 keVand
1.3 keV peaks from the K- and L-shell electron capture
lines in
71
Ge
[38
40]
. The calibration values from a prior
dark matter search were adopted for this work, following
the above described technique and are described in further
detail in Ref.
[38]
. For this dataset, the stability of the
adopted calibration values was then checked using the
electron capture lines as observed in dedicated low back-
ground datasets (i.e., no source placed during the run)
taken at periodic intervals before, after, and between the
photo-neutron runs. These checks showed that the energy
calibration was stable to within 2% throughout the several-
months period when photo-neutron data were collected.
To optimize data quality, several event selection criteria
(cuts) were applied in the analysis. A
live time
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.
Live time 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 live time was removed during
periods of unusually high event rates, based on the average
trigger rate calculated every 30 s of runtime. A few runs
(totaling less than 12 hours) that showed an exceptionally
low noise environment were also removed from the
dataset. The
pretrigger
baseline for each phonon pulse
was recorded for approximately 1 ms before each event.
Excessive fluctuation or large RMS in the pretrigger 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 pretrigger region deviated by more than
4
σ
from
its mean value, as calculated over the course of the run. The
live time 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 signallike and noiselike 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
distribution 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
IONIZATION YIELD MEASUREMENT IN A GERMANIUM
...
PHYS. REV. D
105,
122002 (2022)
122002-5
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
decision 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 trigger condition. This general
technique has been used by SuperCDMS in the past and is
described in Refs.
[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
thresholds matching those in the photo-neutron runs.
They also provided data with a continuous neutron energy
distribution 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 system-
atic 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
described previously by relying on a data driven approach
following the procedures described in Ref.
[33]
.To
compute 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
generated as a linear combination of two subtemplates
modeling the fast and slow components of the absorption
signal described above and in more detail in Ref.
[38]
. The
relative weights of the subtemplates 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 distributions 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 generated events, as a function
of energy. For the final calculation, the efficiency functions
of the constituent time periods were combined together,
weighting each period by both its live time 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 combined 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
5 101520253035404550
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 101520253035404550
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)
Total Phonon Energy (keV)
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) with-
out the Be wafer, i.e., neutron-OFF.
M. F. ALBAKRY
et al.
PHYS. REV. D
105,
122002 (2022)
122002-6
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 distribution.
To test the consistency of the energy spectra over time, a
χ
2
was calculated between pairs of datasets 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 datasets.
IV. SIMULATIONS
A primary input to the likelihood calculation is the
neutron energy probability distribution function (PDF).
GEANT
4 10.6
[43]
was used to simulate the neutron energy
spectrum from the different source configurations. The
GEANT
4
N
eutron
HP
physics model
[44]
along with the
G
4
NDL
4.6
cross section packages were implemented for
the simulation of
1
.
2
×
10
9
neutrons propagating through
the experimental geometry.
G
4
NDL
4.6
is the
JEFF
3.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.
VB
. Therefore, the only
simulated data used in this analysis was the neutron 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 scatter nuclear
recoil (NR) energy deposits from neutrons interacting
directly in the detector and electron recoils (ER) from
neutron-induced gammas are shown separately. The sharp
dropoff in the single-scatter neutron spectrum (
end point
)
corresponds to the maximum energy 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 materials
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 0.5
1 1.5
2 2.5
3 3.5
4 4.5
5
02468101214
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 are
shown in blue, the energy depositions from multiple-scatter
nuclear recoils (NR) by neutrons in the detector are shown in
green. The energy distribution of single-scatter NR by neutrons is
shown in red.
IONIZATION YIELD MEASUREMENT IN A GERMANIUM
...
PHYS. REV. D
105,
122002 (2022)
122002-7
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 entire energy
range. Furthermore, it is subdominant to the Compton
scattering spectrum that results from the gammas 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
GEANT
4
. To study this
uncertainty, the simulation data were regenerated with the
same sample size as used for the primary analysis with the
LEND
[47]
package in place of the
N
eutron
HP
[44]
and
G
4
NDL
4.6
packages
[48]
. The likelihood analysis was
redone using this alternative simulation dataset (see
Sec.
VC1
). The second uncertainty is introduced by the
neutron-nucleus scattering cross section information avail-
able for germanium through the
GEANT
4
package. Figure
5
shows the cross section in germanium as a function of
neutron energy for three different example realizations of
germanium cross-section files from 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 combination of experimental
data and computational models. The resulting cross sec-
tions show some differences, particularly in the neutron
energy region below 50 keV. Using
N
eutron
HP
, we reran the
simulation with 100 different realizations of the neutron
cross sections. Each simulation was for
6
×
10
8
neutrons.
These simulation datasets were used for determining the
systematic uncertainty introduced by this effect as
described further in Sec.
VC1
. The two systematics are
comparable in size, with the one from using the
LEND
library instead of the
N
eutron
HP
and
G
4
NDL
4.6
being slightly
larger.
V. YIELD EXTRACTION
In order to extract the nuclear recoil ionization yield,
we parametrized the conversion from total phonon energy
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 background 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.
VA
and the background model in Sec.
VB
.
Section
VC
discusses the details 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 deposition
was converted into the total phonon energy (
E
p
), the
observable measured by the detectors.
The neutron signal is modeled for the three source
configurations: (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
account. 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 inter-
action. 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 ratio between
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 between the three
realizations. See text body for description of how TENDL cross
section realizations are generated.
M. F. ALBAKRY
et al.
PHYS. REV. D
105,
122002 (2022)
122002-8
the ionization signal and the primary recoil energy, and
ε
γ
¼
3
eV
[51
53]
is the average energy required 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
[6,54
56]
, supported by earlier
measurements
[26]
, provides a semiempirical prediction for
Y
ð
E
r
Þ
of NR interactions for a material of mass number
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.
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
energy 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 detector
resolution to the simulated neutron spectrum. The reso-
lution 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 Table
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
Þ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
σ
2
B
þ
F
ε
γ
E
þ
AE
2
q
:
ð
7
Þ
To determine the values of
F
and
A
in Eq.
(7)
, the
resolutions 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 Sec.
VC1
.
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 inhomo-
geneous, leading to a reduction in the observed NTL
phonon contribution. To account for this a 3D electric
field map was calculated and taken into account when
determining the total energy of the simulated events.
Details 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 fluctuation
effects, we smoothed the energy distributions using a
Savitzky-Golay filter
[57]
before interpolating them using
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 back-
ground 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]
, corre-
sponding to each of the Ge electron shells. The Compton
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
TABLE II. Values of the parameters derived from the energy
resolution model. The energy scale is electron-equivalent (eV
ee
),
where all events are assumed to be electron recoils.
Bias voltage
σ
B
(eV
ee
)
FA
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
IONIZATION YIELD MEASUREMENT IN A GERMANIUM
...
PHYS. REV. D
105,
122002 (2022)
122002-9