# Ionization yield measurement in a germanium CDMSlite detector using photo-
neutron sources

M.F. Albakry  Department of Physics & Astronomy, University of British
Columbia, Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T 2A3, Canada
I. Alkhatib  Department of Physics, University of Toronto, Toronto, ON M5S
1A7, Canada  D.W.P. Amaral  Department of Physics, Durham University, Durham
DH1 3LE, UK  T. Aralis  Division of Physics, Mathematics, & Astronomy,
California Institute of Technology, Pasadena, CA 91125, USA  T. Aramaki
Department of Physics, Northeastern University, 360 Huntington Avenue, Boston,
MA 02115, USA  I.J. Arnquist  Pacific Northwest National Laboratory, Richland,
WA 99352, USA  I. Ataee Langroudy  Department of Physics and Astronomy, and
the Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M
University, College Station, TX 77843, USA  E. Azadbakht  Department of
Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and
Astronomy, Texas A&M University, College Station, TX 77843, USA  S. Banik
School of Physical Sciences, National Institute of Science Education and
Research, Jatni - 752050, India  Homi Bhabha National Institute, Training
School Complex, Anushaktinagar, Mumbai 400094, India  C. Bathurst  Department
of Physics, University of Florida, Gainesville, FL 32611, USA  D.A. Bauer
Fermi National Accelerator Laboratory, Batavia, IL 60510, USA  L.V.S. Bezerra
Department of Physics & Astronomy, University of British Columbia, Vancouver,
BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T 2A3, Canada  R. Bhattacharyya
Department of Physics and Astronomy, and the Mitchell Institute for
Fundamental Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  M.A. Bowles  Department of Physics, South Dakota School of Mines
and Technology, Rapid City, SD 57701, USA  P.L. Brink  SLAC National
Accelerator Laboratory/Kavli Institute for Particle Astrophysics and
Cosmology, Menlo Park, CA 94025, USA  R. Bunker  Pacific Northwest National
Laboratory, Richland, WA 99352, USA  B. Cabrera  Department of Physics,
Stanford University, Stanford, CA 94305, USA  R. Calkins  Department of
Physics, Southern Methodist University, Dallas, TX 75275, USA  R.A. Cameron
SLAC National Accelerator Laboratory/Kavli Institute for Particle Astrophysics
and Cosmology, Menlo Park, CA 94025, USA  C. Cartaro  SLAC National
Accelerator Laboratory/Kavli Institute for Particle Astrophysics and
Cosmology, Menlo Park, CA 94025, USA  D.G. Cerdeño  Department of Physics,
Durham University, Durham DH1 3LE, UK  Instituto de Física Teórica UAM/CSIC,
Universidad Autónoma de Madrid, 28049 Madrid, Spain  Y.-Y. Chang  Department
of Physics, University of California, Berkeley, CA 94720, USA  M. Chaudhuri
School of Physical Sciences, National Institute of Science Education and
Research, Jatni - 752050, India  Homi Bhabha National Institute, Training
School Complex, Anushaktinagar, Mumbai 400094, India  R. Chen  Department of
Physics & Astronomy, Northwestern University, Evanston, IL 60208-3112, USA  N.
Chott  Department of Physics, South Dakota School of Mines and Technology,
Rapid City, SD 57701, USA  J. Cooley  Department of Physics, Southern
Methodist University, Dallas, TX 75275, USA  H. Coombes  Department of
Physics, University of Florida, Gainesville, FL 32611, USA  J. Corbett
Department of Physics, Queen’s University, Kingston, ON K7L 3N6, Canada  P.
Cushman  School of Physics & Astronomy, University of Minnesota, Minneapolis,
MN 55455, USA  F. De Brienne  Département de Physique, Université de Montréal,
Montréal, Québec H3C 3J7, Canada  M.L. di Vacri  Pacific Northwest National
Laboratory, Richland, WA 99352, USA  M.D. Diamond  Department of Physics,
University of Toronto, Toronto, ON M5S 1A7, Canada  E. Fascione  Department of
Physics, Queen’s University, Kingston, ON K7L 3N6, Canada  TRIUMF, Vancouver,
BC V6T 2A3, Canada  E. Figueroa-Feliciano  Department of Physics & Astronomy,
Northwestern University, Evanston, IL 60208-3112, USA  C.W. Fink  Department
of Physics, University of California, Berkeley, CA 94720, USA  K. Fouts  SLAC
National Accelerator Laboratory/Kavli Institute for Particle Astrophysics and
Cosmology, Menlo Park, CA 94025, USA  M. Fritts  School of Physics &
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  G. Gerbier
Department of Physics, Queen’s University, Kingston, ON K7L 3N6, Canada  R.
Germond  Department of Physics, Queen’s University, Kingston, ON K7L 3N6,
Canada  TRIUMF, Vancouver, BC V6T 2A3, Canada  M. Ghaith  Department of
Physics, Queen’s University, Kingston, ON K7L 3N6, Canada  S.R. Golwala
Division of Physics, Mathematics, & Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA  J. Hall  SNOLAB, Creighton Mine #9, 1039
Regional Road 24, Sudbury, ON P3Y 1N2, Canada  Laurentian University,
Department of Physics, 935 Ramsey Lake Road, Sudbury, Ontario P3E 2C6, Canada
B.A. Hines  Department of Physics, University of Colorado Denver, Denver, CO
80217, USA  M.I. Hollister  Fermi National Accelerator Laboratory, Batavia, IL
60510, USA  Z. Hong  Department of Physics, University of Toronto, Toronto, ON
M5S 1A7, Canada  E.W. Hoppe  Pacific Northwest National Laboratory, Richland,
WA 99352, USA  L. Hsu  Fermi National Accelerator Laboratory, Batavia, IL
60510, USA  M.E. Huber  Department of Physics, University of Colorado Denver,
Denver, CO 80217, USA  Department of Electrical Engineering, University of
Colorado Denver, Denver, CO 80217, USA  V. Iyer  School of Physical Sciences,
National Institute of Science Education and Research, Jatni - 752050, India
Homi Bhabha National Institute, Training School Complex, Anushaktinagar,
Mumbai 400094, India  A. Jastram  Department of Physics and Astronomy, and the
Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M
University, College Station, TX 77843, USA  V.K.S. Kashyap  School of Physical
Sciences, National Institute of Science Education and Research, Jatni -
752050, India  Homi Bhabha National Institute, Training School Complex,
Anushaktinagar, Mumbai 400094, India  M.H. Kelsey  Department of Physics and
Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, TX 77843, USA  A. Kubik  SNOLAB,
Creighton Mine #9, 1039 Regional Road 24, Sudbury, ON P3Y 1N2, Canada  N.A.
Kurinsky  SLAC National Accelerator Laboratory/Kavli Institute for Particle
Astrophysics and Cosmology, Menlo Park, CA 94025, USA  R.E. Lawrence
Department of Physics and Astronomy, and the Mitchell Institute for
Fundamental Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  M. Lee  Department of Physics and Astronomy, and the Mitchell
Institute for Fundamental Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  A. Li  Department of Physics & Astronomy, University
of British Columbia, Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T
2A3, Canada  J. Liu  Department of Physics, Southern Methodist University,
Dallas, TX 75275, USA  Y. Liu  Department of Physics & Astronomy, University
of British Columbia, Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T
2A3, Canada  B. Loer  Pacific Northwest National Laboratory, Richland, WA
99352, USA  P. Lukens  Fermi National Accelerator Laboratory, Batavia, IL
60510, USA  D. MacDonell  Department of Physics & Astronomy, University of
British Columbia, Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T
2A3, Canada  D.B. MacFarlane  SLAC National Accelerator Laboratory/Kavli
Institute for Particle Astrophysics and Cosmology, Menlo Park, CA 94025, USA
R. Mahapatra  Department of Physics and Astronomy, and the Mitchell Institute
for Fundamental Physics and Astronomy, Texas A&M University, College Station,
TX 77843, USA  V. Mandic  School of Physics & Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  N. Mast  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  A.J. Mayer  TRIUMF,
Vancouver, BC V6T 2A3, Canada  H. Meyer zu Theenhausen  Institute for
Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT), 76344,
Germany  Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg,
Germany  É. Michaud  Département de Physique, Université de Montréal,
Montréal, Québec H3C 3J7, Canada  E. Michielin  Department of Physics &
Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada
TRIUMF, Vancouver, BC V6T 2A3, Canada  N. Mirabolfathi  Department of Physics
and Astronomy, and the Mitchell Institute for Fundamental Physics and
Astronomy, Texas A&M University, College Station, TX 77843, USA  B. Mohanty
School of Physical Sciences, National Institute of Science Education and
Research, Jatni - 752050, India  Homi Bhabha National Institute, Training
School Complex, Anushaktinagar, Mumbai 400094, India  J.D. Morales Mendoza
Department of Physics and Astronomy, and the Mitchell Institute for
Fundamental Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  S. Nagorny  Department of Physics, Queen’s University, Kingston,
ON K7L 3N6, Canada  J. Nelson  School of Physics & Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  H. Neog  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  V. Novati  Department of
Physics & Astronomy, Northwestern University, Evanston, IL 60208-3112, USA
J.L. Orrell  Pacific Northwest National Laboratory, Richland, WA 99352, USA
M.D. Osborne  Department of Physics and Astronomy, and the Mitchell Institute
for Fundamental Physics and Astronomy, Texas A&M University, College Station,
TX 77843, USA  S.M. Oser  Department of Physics & Astronomy, University of
British Columbia, Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T
2A3, Canada  W.A. Page  Department of Physics, University of California,
Berkeley, CA 94720, USA  R. Partridge  SLAC National Accelerator
Laboratory/Kavli Institute for Particle Astrophysics and Cosmology, Menlo
Park, CA 94025, USA  D.S. Pedreros  Département de Physique, Université de
Montréal, Montréal, Québec H3C 3J7, Canada  R. Podviianiuk  Department of
Physics, University of South Dakota, Vermillion, SD 57069, USA  F. Ponce
Pacific Northwest National Laboratory, Richland, WA 99352, USA  S. Poudel
Department of Physics, University of South Dakota, Vermillion, SD 57069, USA
A. Pradeep  Department of Physics & Astronomy, University of British Columbia,
Vancouver, BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T 2A3, Canada  M. Pyle
Department of Physics, University of California, Berkeley, CA 94720, USA
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA  W. Rau
TRIUMF, Vancouver, BC V6T 2A3, Canada  E. Reid  Department of Physics, Durham
University, Durham DH1 3LE, UK  R. Ren  Department of Physics & Astronomy,
Northwestern University, Evanston, IL 60208-3112, USA  T. Reynolds  Department
of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada  A. Roberts
Department of Physics, University of Colorado Denver, Denver, CO 80217, USA
A.E. Robinson  Département de Physique, Université de Montréal, Montréal,
Québec H3C 3J7, Canada  T. Saab  Department of Physics, University of Florida,
Gainesville, FL 32611, USA  B. Sadoulet  Department of Physics, University of
California, Berkeley, CA 94720, USA  Lawrence Berkeley National Laboratory,
Berkeley, CA 94720, USA  I. Saikia  Department of Physics, Southern Methodist
University, Dallas, TX 75275, USA  J. Sander  Department of Physics,
University of South Dakota, Vermillion, SD 57069, USA  A. Sattari  Department
of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada  A. Scarff
Department of Physics & Astronomy, University of British Columbia, Vancouver,
BC V6T 1Z1, Canada  TRIUMF, Vancouver, BC V6T 2A3, Canada  B. Schmidt
Department of Physics & Astronomy, Northwestern University, Evanston, IL
60208-3112, USA  R.W. Schnee  Department of Physics, South Dakota School of
Mines and Technology, Rapid City, SD 57701, USA  S. Scorza  SNOLAB, Creighton
Mine #9, 1039 Regional Road 24, Sudbury, ON P3Y 1N2, Canada  Laurentian
University, Department of Physics, 935 Ramsey Lake Road, Sudbury, Ontario P3E
2C6, Canada  B. Serfass  Department of Physics, University of California,
Berkeley, CA 94720, USA  D.J. Sincavage  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  C. Stanford  Department
of Physics, Stanford University, Stanford, CA 94305, USA  J. Street
Department of Physics, South Dakota School of Mines and Technology, Rapid
City, SD 57701, USA  F.K. Thasrawala  Institut für Experimentalphysik,
Universität Hamburg, 22761 Hamburg, Germany  D. Toback  Department of Physics
and Astronomy, and the Mitchell Institute for Fundamental Physics and
Astronomy, Texas A&M University, College Station, TX 77843, USA  R. Underwood
Department of Physics, Queen’s University, Kingston, ON K7L 3N6, Canada
TRIUMF, Vancouver, BC V6T 2A3, Canada  S. Verma  Department of Physics and
Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, TX 77843, USA  A.N. Villano  Department
of Physics, University of Colorado Denver, Denver, CO 80217, USA  B. von
Krosigk  Institute for Astroparticle Physics (IAP), Karlsruhe Institute of
Technology (KIT), 76344, Germany  Institut für Experimentalphysik, Universität
Hamburg, 22761 Hamburg, Germany  S.L. Watkins  Department of Physics,
University of California, Berkeley, CA 94720, USA  O. Wen  Division of
Physics, Mathematics, & Astronomy, California Institute of Technology,
Pasadena, CA 91125, USA  Z. Williams  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  M.J. Wilson  Institute
for Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT),
76344, Germany  Department of Physics, University of Toronto, Toronto, ON M5S
1A7, Canada  J. Winchell  Department of Physics and Astronomy, and the
Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M
University, College Station, TX 77843, USA  K. Wykoff  Department of Physics,
South Dakota School of Mines and Technology, Rapid City, SD 57701, USA  S.
Yellin  Department of Physics, Stanford University, Stanford, CA 94305, USA
B.A. Young  Department of Physics, Santa Clara University, Santa Clara, CA
95053, USA  T.C. Yu  SLAC National Accelerator Laboratory/Kavli Institute for
Particle Astrophysics and Cosmology, Menlo Park, CA 94025, USA  B. Zatschler
Department of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada  S.
Zatschler  Department of Physics, University of Toronto, Toronto, ON M5S 1A7,
Canada  A. Zaytsev  Institute for Astroparticle Physics (IAP), Karlsruhe
Institute of Technology (KIT), 76344, Germany  Institut für
Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany  E. Zhang
Department of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada  L.
Zheng  Department of Physics and Astronomy, and the Mitchell Institute for
Fundamental Physics and Astronomy, Texas A&M University, College Station, TX
77843, USA  S. Zuber  Department of Physics, University of California,
Berkeley, CA 94720, USA

###### Abstract

Two photo-neutron sources, 88Y9Be and 124Sb9Be, 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 of $\sim$ 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.

††preprint:

##  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 “ionization
yield”. For detectors 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 recoil signal
induced by coherent elastic neutrino nucleus scattering [3, 4, 5].

Figure 1: Literature data on ionization yield measurements as a function of
recoil energy in germanium [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
19, 20, 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 similar to the averaged
yield over all detectors. However, variation across individual detectors was
more than the statistical uncertainty of the measurement shown here [12].
EDELWEISS has reported two measurements: (i) discrete yield measurements using
a 252Cf neutron source [15], and (ii) continuous measurement using an 241Am9Be
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
88Y9Be 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 [21]. The prediction by
Lindhard et. al. [22, 23] is also shown for comparison.

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 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 reference model in the field of direct dark matter searches.
Experiments often use a value of the $k$ parameter (explained in Sec. V.1),
whose physical origin is related to the electronic stopping power [24], 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 [7], electric field dependence [25], and
temperature 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/c2). 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. V.1 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 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 9Be [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 spectrum and a two-parameter yield model to extract the
ionization yield from the data. We discuss the experimental setup and data
sets in Section II, the data selection criteria 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 Section 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 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 between 40 and 50
mK. Further details about the experimental setup and the detector design can
be found in Ref. [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 secondary 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, 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)).

Figure 2: (a) The experimental setup with the shielding configuration. 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 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.

Two gamma sources, 88Y (half-life 106.6 days) and 124Sb (half-life 60.2 days),
were deployed sequentially for the photo-neutron calibration. Each source was
placed on top of a solid 9Be disk. 9Be is the only stable isotope of Be that
occurs naturally in non-trace 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
88Y is at 1.8 MeV. In the case of 124Sb, the most prominent gamma above
threshold is 1.69 MeV. In the photoproduction process, the 9Be nucleus absorbs
a gamma and emits a neutron to become a 8Be resonant state. The resulting
neutrons emitted from 9Be, over multiple occurrences of this process, are
nearly mono-energetic, with an average energy of 152 keV [19] for 88Y and 24
keV [37] for 124Sb. The diameter of the 9Be disc was 5 cm with a thickness of
0.2 cm. The purity level of the Be disc was 98.5% 9Be. A source holder was
made to constrain the position of the gamma source with respect to the 9Be
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 88Y9Be source is
$\sim 8$ keV, and for the 124Sb9Be source is $\sim 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 Geant4
simulation; instead we used a mono-energetic 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 $\sim$144 days between June 5th,
2015 and October 26th, 2015. Data were recorded with the 9Be disk (neutron-ON)
and without the 9Be disk (neutron-OFF). The latter served to measure the gamma
background. A summary of the data-taking condition is shown in Table 1. The
data was first taken with the T5Z2 detector, and subsequently with the T2Z1
detector. As the activity of the 124Sb source ($\tau_{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.

Table 1: Summary of the photo-neutron data-taking runs. Source | Duration | Detector | Voltage
---|---|---|---
124Sb and 124Sb9Be | 62 days | T5Z2 | 70 V
88Y and 88Y9Be | 42 days | T5Z2 | 70 V
88Y and 88Y9Be | 38 days | T2Z1 | 25 V

##  III Data reconstruction, selection and efficiencies

The phonon pulse shape observed in the various channels 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 ($\sim 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 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 “non-stationary” optimal 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 between 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 frequency 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 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 $\chi^{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 252Cf. This led
to a neutron capture process, creating unstable 71Ge. 71Ge decays via electron
capture to 71Ga, releasing predominantly 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 capture lines in 71Ge [38, 39, 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 background 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 was collected.

To optimize data quality, several event selection criteria (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 $\sim 1$% of livetime
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 “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 deviated by more than $4\sigma$ from
its mean value, as calculated over the course of the run. The livetime cuts
removed $\sim$10% of the events from the 70 V datasets, and $\sim$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 $\chi^{2}$ goodness of
fit was calculated for each event. The $\chi^{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
$\chi^{2}$ or with a rise time outside a 230 $\mu$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 ionization 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 $\chi^{2}$ values ($\chi^{2}>$6250 for data from
T2Z1, and $\chi^{2}>$50000 for data from T5Z2). Thus, events with a high
charge $\chi^{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 $\Delta OF\chi^{2}_{\rm noise}=OF\chi^{2}_{\rm
pulse}-OF\chi^{2}_{\rm noise}$ was constructed, where the subscript “noise”
stands for either LFN or glitch. A cut at an appropriate level of signal
acceptance (3$\sigma$) was applied on the $\Delta OF\chi^{2}_{\rm 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 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.

Figure 3:  Combined and smoothed cut efficiency as a function of energy above
the analysis threshold for each data set with a $1\sigma$ uncertainty. (a)
With the Be wafer i.e. neutron-ON, and (b) without the Be wafer i.e. neutron-
OFF.

_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 references [38] [42]. The trigger efficiency for this analysis
is based on a combination of data taken with a 252Cf 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 systematic 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$\equiv$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 sub-templates modeling the fast
and slow components of the absorption signal 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 distributions as observed for the measured total phonon
energy and the phonon OF $\chi^{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 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 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 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 $\chi^{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 $\chi^{2}$ was calculated between pairs of data
sets drawn from the eleven time periods and the same photo-neutron source. The
$\chi^{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.

##  IV Simulations

Figure 4:  The simulated recoil energy distributions from (a) 124Sb9Be, and
(b) 88Y9Be 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$\times 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.2. 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
124Sb9Be and 88Y9Be 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 drop-off in the single-scatter neutron
spectrum (“endpoint”) corresponds to the maximum energy a neutron of the given
input energy can transfer to a Ge nucleus in a single interaction ($\sim$1.3
keV for the 24 keV neutrons from 124Sb9Be and $\sim$8.2 keV for the 152 keV
neutrons from 88Y9Be; 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 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
Geant4. 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 NeutronHP [44] and G4NDL4.6 packages [48]. The
likelihood analysis was re-done using this alternative simulation dataset (see
Sec. V.3.1). The second uncertainty is introduced 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 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-sections show
some differences, particularly in the neutron energy region below 50 keV.
Using NeutronHP, we re-ran the simulation with 100 different realizations of
the neutron cross-sections. Each simulation was for 6$\times 10^{8}$ neutrons.
These simulation datasets were used for determining the systematic uncertainty
introduced by this effect as described further in Sec. V.3.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.

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

##  V Yield extraction

In order to extract the nuclear recoil ionization yield, we parameterized 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. V.1, and the background model in Sec. V.2. Section V.3
discusses the details of the likelihood fit, the propagation of statistical
and systematic uncertainties, and the final results.

###  V.1 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)
124Sb9Be with the detector biased at 70 V, (ii) 88Y9Be with the detector
biased at 70 V, and (iii) 88Y9Be 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_{\text{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_{\text{NTL}}=N_{e/h}e\Delta
V=Y(E_{r})\frac{E_{r}}{\epsilon_{\gamma}}e\Delta V,$ (1)

where $N_{e/h}$ is the number of $e^{-}/h^{+}$pairs created by the recoil,
$\Delta V$ is the bias voltage, and $e$ is the electronic charge. $Y(E_{r})$
is the ionization yield, defined as the ratio between the ionization signal
and the primary recoil energy, and $\epsilon_{\gamma}=3$ eV [51, 52, 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\Delta V/\epsilon_{\gamma}),$ (2)

where the dependence of $E_{p}$ on the parameter of interest, $Y(E_{r})$,
appears clearly.

The Lindhard theory [54, 22, 55, 56], supported by earlier measurements [23],
provides a semi-empirical prediction for $Y(E_{r})$ of NR interactions for a
material of mass number $A$ and atomic number $Z$:

$Y(E_{r})=k\frac{g(\epsilon)}{1+kg(\epsilon)},$ (3)

where:

$g(\epsilon)=3\epsilon^{0.15}+0.7\epsilon^{0.6}+\epsilon;~{}~{}\epsilon=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_{\text{{low}}}+\frac{k_{\text{high}}-k_{\text{low}}}{E_{\text{high}}-E_{\text{low}}}(E_{r}-E_{\text{low}}).$
(5)

$E_{\text{low}}$ is equal to 0.39 keV and $E_{\text{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_{\text{low}}$ and
$k_{\text{high}}$ as the two components of a vector $\vec{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 resolution model ($\sigma_{T}$) used was:

$\sigma_{T}^{2}(E)=\sigma_{B}^{2}+\sigma_{F}^{2}(E)+\sigma_{D}^{2}(E).$ (6)

where $\sigma_{B}$ is the baseline noise due to the electronics, $\sigma_{F}$
is the variance in the number of $e^{-}/h^{+}$pairs produced in a recoil
event, and $\sigma_{D}$ is an empirical term that models detector effects and
is given by $\sigma_{D}^{2}=AE^{2}$.

The values of $\sigma_{B}$ are obtained from the energy distribution of events
randomly triggered in the detector and are listed in Tab. 2. The second term
can be written as
$\sigma_{F}^{2}=FN_{e/h}\epsilon_{\gamma}^{2}=F\epsilon_{\gamma}E$, where $F$
is the Fano factor. The resolution model can thus be rewritten as:

$\sigma_{T}(E)=\sqrt{\sigma_{B}^{2}+F\epsilon_{\gamma}E+AE^{2}},$ (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 71Ge were
determined by first fitting them with a Gaussian and obtaining their 1$\sigma$
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 2
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.3.1.

Table 2: Values of the parameters derived from the energy resolution model. The energy scale is electron-equivalent (eVee), where all events are assumed to be electron-recoils. Bias voltage |  $\sigma_{B}$ (eVee) | $F$ | $A$
---|---|---|---
70 V | 8.31$\pm$0.08 | 0.19$\pm$0.07 | 0.0095$\pm$0.0015
25 V | 18.71$\pm$0.16 | 0.27$\pm$0.02 | 0.0107$\pm$0.0005

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 inhomogeneous, 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. 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.

Figure 6: Neutron-OFF energy spectrum for the three data sets: (a) 124Sb9Be
with T5Z2 biased at 70 V, (b) 88Y9Be with T5Z2 biased at 70 V and (c) 88Y9Be
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.

###  V.2 Modeling electron recoil backgrounds

The dominant background that affected this measurement 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], corresponding to each of the Ge electron shells.
The Compton model in the region of interest is:

$\displaystyle f_{C}(E)=a_{0}[1$ $\displaystyle+S_{M}\Theta(E,\mu_{M})]$ (8)
$\displaystyle+S_{L}a_{0}[1$
$\displaystyle+S_{M}\Theta(E,\mu_{M})]\Theta(E,\mu_{L})$
$\displaystyle+S_{K}a_{0}[1$
$\displaystyle+S_{M}\Theta(E,\mu_{M})]\Theta(E,\mu_{L})\Theta(E,\mu_{K});$

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 $\mu$ is the central energy of the step, which
was determined via a Geant4 simulation. $\Theta$ is the error function:

$\Theta(E,\mu)=\frac{1}{2}\left[1+\text{erf}\left(\frac{E-\mu}{\sqrt{2}\sigma(\mu)}\right)\right],$
(9)

with the detector resolution $\sigma(\mu)$.

The electron capture peaks were modeled using Gaussian functions. Following
Ref. [59], we assumed the fractional 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 determined with a fit to the
neutron-OFF data for this background was the amplitude of the K-shell peak. We
then constrained the other amplitudes using the fractional ratio. While data
collected at 70 V are sensitive to the L-, K-, and M-shells, for the 25 V data
the M-shell step occurs 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.3.1.

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 remove any effects on its shape due to the neutron-OFF cut efficiency.

###  V.3 Likelihood analysis

The background and signal PDFs allowed us to define the negative log
likelihood function as:

$-\text{ln}\,\mathcal{L}=-\sum_{D=1}^{3}\sum_{i=1}^{N_{D}}\text{ln}(f_{D}\nu_{D}(E_{i},\vec{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, $\nu_{D}(E,\vec{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 minimized
using the MINUIT [60] package, were the three neutron contribution fractions
${f}_{D}$ and the Lindhard parameters $\vec{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 $\chi^{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
$\chi^{2}$/(degrees of freedom) was found to be $137/95$, $141/95$ and
$129/95$ for 124Sb at 70 V, 88Y at 70 V and 88Y at 25 V respectively. The best
fit values of $\vec{k}$ are shown in Tab. 3.

Figure 7:  Energy spectrum for the three data sets, (a) 124Sb9Be with T5Z2
biased at 70 V, (b) 88Y9Be with T5Z2 biased at 70 V and (c) 88Y9Be 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.

In order to study whether the 1-parameter (standard Lindhard) or 2-parameter
(see Eq. 5) $k$ model better described 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 spectrum, 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 provides 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$\sigma$. Based on this result, we report only the result for the 2-parameter
modified Lindhard model.

####  V.3.1 Calculation of uncertainties

We took into account the statistical uncertainties due to the finite size of
the simulated neutron spectrum and of the experimental neutron-ON and neutron-
OFF spectra by performing fits to simulated experiments. The distributions of
the fit results were then fit with Gaussian distributions, and the standard
deviation of each Gaussian 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 deviation 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 resolution model was
determined from a fit to data (predominantly electron recoils) and it was
affected by statistical and systematic uncertainties. In particular, a
previous measurement [61] showed that the Fano factor for nuclear recoils
could be significantly higher than in electron 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 forcing 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 associated with
the Fano factor. In order to evaluate the uncertainty on the background model
shape, described in Sec. V.2, we recalculated the yield by using as background
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 section 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 section libraries was used as an
estimate of this systematic uncertainty. This approach does not take into
account possible correlated biases between the two packages.

Table 3: Summary of the $k_{\text{{low}}}$ and $k_{\text{{high}}}$ fit results along with their statistical and systematic uncertainties | Best fit value | Stat. uncertainty | Sys. uncertainty
---|---|---|---
klow | 0.040 | 0.005 | 0.008
khigh | 0.142 | 0.011 | 0.026

##  VI Results

The fitted Lindhard model parameters and their respective uncertainties were
propagated to the ionization yield $Y(E_{r})$ vs. nuclear recoil energy plane
in Fig. 8. While the resulting ionization yield was compatible with the
Lindhard model with the standard $k=0.157$ value at the higher end of our
energy range, our data prefer a significantly lower yield at the lower
energies. As indicated by the result of the likelihood ratio test described in
Sec. V.3, the data were better described by a model with a linear energy
variation for the $\vec{k}$ parameters, which departs from the standard
Lindhard model. Although the final result was obtained with a simultaneous fit
to the 70V and 25V datasets, we performed a cross-check where we fit each set
separately. The results were consistent within the statistical uncertainty.
The dominant sources of uncertainty for this measurement were the statistical
uncertainty due to the size of the experimental data set, the systematic
uncertainty due to the modeling of the background as well as due to the
elastic neutron-nucleus scattering cross section for germanium. This
highlights the importance of developing an accurate and precise simulation for
this type of analysis.

Figure 8:  (a) The ionization yield with 1 $\sigma$ uncertainty from the best
fit values of the 2-component $k$-parameter Lindhard model as a function of
the nuclear recoil energy in germanium. The blue line shows the standard
Lindhard model with $k=0.157$. (b) The contribution of various sources of
statistical and systematic uncertainties to the ionization yield in germanium.
The lower limit of the x-axis represents the analysis threshold converted to
recoil energy using the best-fit value.

##  VII Conclusions

This paper presents a nuclear recoil induced ionization yield measurement for
first generation, high-voltage, SuperCDMS germanium detectors operated at
voltages of 25 V and 70 V and a temperature of approximately 50 mK. Our
analysis has robustly shown that under these conditions, the ionization yield
in the 1–7 keV recoil energy region has a behavior that falls off
significantly faster than the standard Lindhard model. These data also differ
significantly from those of other existing measurements in this energy regime.
However, previous measurements exhibit a lack of consistency among themselves,
suggesting that there are systematic effects that remain to be understood
within the field as a whole.

Our model for the yield function is a first order approximation to a generic
function that tends to the Lindhard functional form at high energies. This is
achieved by allowing the $k$ parameter to have a linear energy dependence in
the region below 7 keV. To improve on this technique, future measurements
could employ thinner detectors to reduce the population of multiple neutron
scatters, gather a significantly larger data sample (to reduce the statistical
uncertainty), and achieve a better understanding of low-energy neutron
scattering cross sections implemented in the simulation (thus reducing the
largest source of systematic uncertainty).

In the case of experiments searching for the elusive scattering of DM off
nuclei, the nuclear recoil ionization yield affects the signal spectrum
substantially, thus injecting uncertainty into DM searches. In the low energy
region, measurements in literature for germanium and silicon have not always
yielded consistency with Lindhard or with each other. Furthermore, this low-
energy region is particularly important for recent dark matter searches with
solid state detectors. These inconsistencies may be due to temperature,
electric field, or other effects that are generally unaccounted for in
measurements to date (see Fig. 1). Thus, the situation suggests that multiple
measurements, performed with a variety of experimental approaches and devices,
are required to shed light on possible, heretofore, unquantified systematic
effects in these studies. Until those effects are clearly understood, the
existing measurements will continue to be critical for determining the
ionization yield under the specific operating conditions of individual dark
matter experiments. The current report provides the community a new set of
data via an in-situ photon-neutron measurement method, presented in comparison
to the Lindhard model.

###### Acknowledgements.

The SuperCDMS collaboration gratefully acknowledges technical assistance from
the staff of the Soudan Underground Laboratory and the Minnesota Department of
Natural Resources. The CDMSlite and iZIP detectors were fabricated in the
Stanford Nanofabrication Facility, which is a member of the National
Nanofabrication Infrastructure Network, sponsored and supported by the NSF.
Funding and support were received from the National Science Foundation, the
U.S. Department of Energy (DOE), Fermilab URA Visiting Scholar Grant No.
15-S-33, NSERC Canada, the Canada First Excellence Research Fund, the Arthur
B. McDonald Institute (Canada), the Department of Atomic Energy Government of
India (DAE), the Department of Science and Technology (DST, India) and the DFG
(Germany) - Project No. 420484612 and under Germany’s Excellence Strategy -
EXC 2121 “Quantum Universe” – 390833306. Femilab is operated by Fermi Research
Alliance, LLC, SLAC is operated by Stanford University, and the PNNL is
operated by the Battelle Memorial Institute, each for the U.S. Department of
Energy under contracts DE-AC02-37407CH11359, DE-AC02-76SF00515, and DE-
AC05-76RL01830, respectively.

## References

  * Lewin and Smith [1996] J. Lewin and P. Smith, Astroparticle Physics 6, 87 (1996).
  * Agnese  _et al._ [2017] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 95, 082002 (2017).
  * Freedman [1974] D. Z. Freedman, Phys. Rev. D 9, 1389 (1974).
  * Drukier and Stodolsky [1984] A. Drukier and L. Stodolsky, Phys. Rev. D 30, 2295 (1984).
  * Akimov  _et al._ [2017] D. Akimov _et al._ , Science 357, 1123 (2017).
  * Jones and Kraner [1971] K. W. Jones and H. W. Kraner, Phys. Rev. C 4, 125 (1971).
  * Jones and Kraner [1975] K. W. Jones and H. W. Kraner, Phys. Rev. A 11, 1347 (1975).
  * Aalseth  _et al._ [2013] C. E. Aalseth _et al._ (CoGeNT), Phys. Rev. D88, 012002 (2013), arXiv:1208.5737 [astro-ph.CO] .
  * Xi-chao  _et al._ [2011] R. Xi-chao _et al._ (CDEX-TEXONO), Nuclear recoil quenching factor measurement for HPGe detector, presented at Application of Germanium Detector in fundamental research, Tsinghua University , 24–26 (2011).
  * Messous  _et al._ [1995] Y. Messous _et al._ , Astroparticle Physics 3, 361 (1995).
  * Ahmed  _et al._ [2011] Z. Ahmed _et al._ (CDMS Collaboration), Phys. Rev. Lett. 106, 131302 (2011), arXiv:1011.2482 .
  * Moore [2012] D. C. Moore,  _A Search for Low-Mass Dark Matter with the Cryogenic Dark Matter Search and the Development of Highly Multiplexed Phonon-Mediated Particle Detectors_ , Ph.D. thesis, California Institute of Technology, California (2012).
  * Chasman  _et al._ [1968] C. Chasman _et al._ , Phys. Rev. Lett. 21, 1430 (1968).
  * Shutt  _et al._ [1992] T. Shutt _et al._ , Phys. Rev. Lett. 69, 3425 (1992).
  * Benoit  _et al._ [2007] A. Benoit _et al._ (EDELWEISS), Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 577, 558 (2007).
  * Sattler  _et al._ [1966] Sattler _et al._ , Phys. Rev. 143, 588 (1966).
  * Baudis  _et al._ [1998] L. Baudis _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 418, 348 (1998).
  * Simon  _et al._ [2003] E. Simon _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 507, 643 (2003).
  * Scholz  _et al._ [2016] B. J. Scholz _et al._ , Phys. Rev. D94, 122003 (2016), arXiv:1608.03588 [physics.ins-det] .
  * Armengaud  _et al._ [2017] E. Armengaud _et al._ (EDELWEISS), Journal of Instrumentation 12 (08), P08010.
  * Collar  _et al._ [2021] J. I. Collar _et al._ , Phys. Rev. D 103, 122003 (2021).
  * Lindhard  _et al._ [1963a] J. Lindhard _et al._ , Kgl. Danske Videnskab. Selskab. Mat. Fys. Medd. 33, (1963a).
  * Barker and Mei [2012] D. Barker and D.-M. Mei, Astroparticle Physics 38, 1 (2012).
  * Horsfield  _et al._ [2016] A. P. Horsfield, A. Lim, W. M. C. Foulkes, and A. A. Correa, Phys. Rev. B 93, 245106 (2016).
  * Penn  _et al._ [1996] M. J. Penn, B. L. Dougherty, B. Cabrera, R. M. Clarke, and B. A. Young, Journal of Applied Physics 79, 8179 (1996), https://doi.org/10.1063/1.362552 .
  * Wei and Mei [2017] W.-Z. Wei and D.-M. Mei, Journal of Instrumentation 12 (04), P04022.
  * Agnese  _et al._ [2013] R. Agnese _et al._ (SuperCDMS), Appl. Phys. Lett. 103, 164105 (2013), arXiv:1305.2405 [physics.ins-det] .
  * Agnese  _et al._ [2014a] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 112, 041302 (2014a).
  * Agnese  _et al._ [2016] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 116, 071301 (2016).
  * Agnese  _et al._ [2019a] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D99, 062001 (2019a), arXiv:1808.09098 [astro-ph.CO] .
  * Wattenberg [1947] A. Wattenberg, Phys. Rev. 71, 497 (1947).
  * Robinson [2016] A. E. Robinson, Phys. Rev. C 94, 024613 (2016).
  * Agnese  _et al._ [2019b] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 99, 062001 (2019b).
  * Agnese  _et al._ [2014b] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 112, 041302 (2014b).
  * Neganov  _et al._ [1993] B. Neganov _et al._ , J Low Temp Phys 93, 417–422 (1993).
  * Luke  _et al._ [1990] P. Luke _et al._ , Nucl. Instrum. Meth. A 289, 406 (1990).
  * Chavarria  _et al._ [2016] A. E. Chavarria _et al._ , Phys. Rev. D 94, 082007 (2016).
  * Agnese  _et al._ [2018] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 97, 022002 (2018).
  * Hampel and Remsberg [1985] W. Hampel and L. P. Remsberg, Phys. Rev. C 31, 666 (1985).
  * Bearden and Burr [1967] J. A. Bearden and A. F. Burr, Rev. Mod. Phys. 39, 125 (1967).
  * Chagani  _et al._ [2012] H. Chagani _et al._ , J. Low Temp. Phys. 167, 1125 (2012).
  * Akerib  _et al._ [2010] D. S. Akerib _et al._ (CDMS), Phys. Rev. D 82, 122004 (2010), arXiv:1010.4290 [astro-ph.CO] .
  * Allison  _et al._ [2016] J. Allison _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 835, 186 (2016).
  * Koi [2013] T. Koi, in  _18th GEANT4 collaboration meeting_ (2013).
  * Mendoza and Cano-Ott [2020] E. Mendoza and D. Cano-Ott, Neutron library g4ndl4.6: validation and computing performance, geant4 hadronic meeting (2020).
  * Mendoza and Cano-Ott [2018] E. Mendoza and D. Cano-Ott, Update of the evaluated neutron cross section libraries for the geant4 code (2018).
  * GEANT4 [2020] GEANT4, Low energy nuclear data package (2020).
  * Meeting [2020] G. H. P. W. G. Meeting, Neutron library g4ndl4.6 : validation and computing performance (2020).
  * Koning and Rochman [2012] A. Koning and D. Rochman, Nuclear Data Sheets 113, 2841 (2012), special Issue on Nuclear Reaction Data.
  * Rochman and Koning [2011] D. Rochman and A. Koning, Nuclear Science and Engineering 169, 68 (2011).
  * Antman and all [1966] S. Antman and all, Nuclear Instruments and Methods 40, 272 (1966).
  * Pehl  _et al._ [1968] R. Pehl, F. Goulding, D. Landis, and M. Lenzlinger, Nuclear Instruments and Methods 59, 45 (1968).
  * Emery and Rabson [1965] F. E. Emery and T. A. Rabson, Phys. Rev. 140, A2089 (1965).
  * Lindhard and Scharff [1961] J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).
  * Lindhard  _et al._ [1968] J. Lindhard, V. Nielsen, and M. Scharff,  _Approximation Method in Classical Scattering by Screened Coulomb Fields: Notes on Atomic Collisions, 1_, Matematisk-fysiske meddelelser (Munksgaard, 1968).
  * Lindhard  _et al._ [1963b] Lindhard _et al._ , Kgl. Danske Videnskab., Selskab. Mat. Fys. Medd. 33, (1963b).
  * Savitzky and Golay [1964] A. Savitzky and M. J. E. Golay, Analytical Chemistry 36, 1627 (1964), https://doi.org/10.1021/ac60214a047 .
  * Ramanathan  _et al._ [2017] K. Ramanathan _et al._ , Phys. Rev. D 96, 042002 (2017), arXiv:1706.06053 [physics.ins-det] .
  * Schönfeld [1998] E. Schönfeld, Applied Radiation and Isotopes 49, 1353 (1998).
  * James and Roos [1975] F. James and M. Roos, Computer Physics Communications 10, 343 (1975).
  * Dougherty [1992] B. L. Dougherty, Phys. Rev. A 45, 2104 (1992).
  * Papp  _et al._ [2004] T. Papp, M.-C. Lepy, J. Plagnard, G. Kalinka, and E. Papp-Szabó, X‐Ray Spectrometry 34, 106 (2004).
  * Lowe [1997] B. Lowe, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 399, 354 (1997).