of 7
First Measurement of the Nuclear-Recoil Ionization Yield in Silicon at 100 eV
M. F. Albakry,
1,2
I. Alkhatib,
3
D. Alonso,
4,5
D. W. P. Amaral,
6
T. Aralis,
7
T. Aramaki,
8
I. J. Arnquist,
9
I. Ataee Langroudy,
10
E. Azadbakht,
10
S. Banik,
11
C. Bathurst,
12
R. Bhattacharyya,
10
P. L. Brink,
13
R. Bunker,
9
B. Cabrera,
14
R. Calkins,
15
R. A. Cameron,
13
C. Cartaro,
13
D. G. Cerdeño,
4,5
Y.-Y. Chang,
7
M. Chaudhuri,
11
R. Chen,
16
N. Chott,
17
J. Cooley,
18,15
H. Coombes,
12
J. Corbett,
19
P. Cushman,
20
S. Das,
11
F. De Brienne,
21
M. Rios,
4,5
S. Dharani,
22,23
M. L. di Vacri,
9
M. D. Diamond,
3
M. Elwan,
12
E. Fascione,
19,2
E. Figueroa-Feliciano,
16
C. W. Fink,
24
K. Fouts,
13
M. Fritts,
20
G. Gerbier,
19
R. Germond,
19,2
M. Ghaith,
25
S. R. Golwala,
7
J. Hall,
18,26
S. A. S. Harms,
3
N. Hassan,
21
B. A. Hines,
27
Z. Hong,
3
E. W. Hoppe,
9
L. Hsu,
28
M. E. Huber,
27,29
V. Iyer,
3
V. K. S. Kashyap,
11
M. H. Kelsey,
10
A. Kubik,
18
N. A. Kurinsky,
13
M. Lee,
10
M. Litke,
15
J. Liu,
15
Y. Liu,
1,2
B. Loer,
9
E. Lopez Asamar,
4,5
P. Lukens,
28
D. B. MacFarlane,
13
R. Mahapatra,
10
N. Mast,
20
A. J. Mayer,
2
H. Meyer zu Theenhausen,
22
É. Michaud,
21
E. Michielin,
1,2
N. Mirabolfathi,
10
B. Mohanty,
11
B. Nebolsky,
16
J. Nelson,
20
H. Neog,
20
V. Novati,
16
J. L. Orrell,
9
M. D. Osborne,
10
S. M. Oser,
1,2
W. A. Page,
24
L. Pandey,
30
S. Pandey,
20
R. Partridge,
13
D. S. Pedreros,
21
L. Perna,
3
R. Podviianiuk,
30
F. Ponce,
9
S. Poudel,
30
A. Pradeep,
1,2
M. Pyle,
24,31
W. Rau,
2
E. Reid,
6
R. Ren ,
16
,*
T. Reynolds,
3
E. Tanner,
20
A. Roberts,
27
A. E. Robinson,
21
T. Saab,
12
D. Sadek,
12
B. Sadoulet,
24,31
S. P. Sahoo,
10
I. Saikia,
15
J. Sander,
30
A. Sattari,
3
B. Schmidt,
16
R. W. Schnee,
17
S. Scorza,
18,26
B. Serfass,
24
S. S. Poudel,
9
D. J. Sincavage,
20
P. Sinervo,
3
Z. Speaks,
12
J. Street,
17
H. Sun,
12
G. D. Terry,
30
F. K. Thasrawala,
23
D. Toback,
10
R. Underwood,
19,2
S. Verma,
10
A. N. Villano,
27
B. von Krosigk,
22
S. L. Watkins,
24
O. Wen,
7
Z. Williams,
20
M. J. Wilson,
22
J. Winchell,
10
K. Wykoff,
17
S. Yellin,
14
B. A. Young,
32
T. C. Yu,
13
B. Zatschler,
3
S. Zatschler,
3
A. Zaytsev,
22
A. Zeolla,
12
E. Zhang,
3
L. Zheng,
10
Y. Zheng,
16
A. Zuniga,
3
(SuperCDMS Collaboration)
P. An,
33
P. S. Barbeau,
33
S. C. Hedges,
34
L. Li,
33
and J. Runge
33
1
Department of Physics & 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
Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049 Madrid, Spain
5
Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain
6
Department of Physics, Durham University, Durham DH1 3LE, United Kingdom
7
Division of Physics, Mathematics, & Astronomy, California Institute of Technology, Pasadena, California 91125, USA
8
Department of Physics, Northeastern University, 360 Huntington Avenue, Boston, Massachusetts 02115, USA
9
Pacific Northwest National Laboratory, Richland, Washington 99352, USA
10
Department of Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,
Texas A&M University, College Station, Texas 77843, USA
11
School of Physical Sciences, National Institute of Science Education and Research, HBNI, Jatni - 752050, India
12
Department of Physics, University of Florida, Gainesville, Florida 32611, USA
13
SLAC National Accelerator Laboratory/Kavli Institute for Particle Astrophysics and Cosmology, Menlo Park, California 94025, USA
14
Department of Physics, Stanford University, Stanford, California 94305, USA
15
Department of Physics, Southern Methodist University, Dallas, Texas 75275, USA
16
Department of Physics & Astronomy, Northwestern University, Evanston, Illinois 60208-3112, USA
17
Department of Physics, South Dakota School of Mines and Technology, Rapid City, South Dakota 57701, USA
18
SNOLAB, Creighton Mine #9, 1039 Regional Road 24, Sudbury, Ontario P3Y 1N2, Canada
19
Department of Physics, Queen
s University, Kingston, Ontario K7L 3N6, Canada
20
School of Physics & Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA
21
D ́
epartement de Physique, Universit ́
e de Montr ́
eal, Montr ́
eal, Qu ́
ebec H3C 3J7, Canada
22
Institute for Astroparticle Physics (IAP), Karlsruhe Institute of Technology (KIT), 76344 Eggenstein-Leopoldshafen, Germany
23
Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany
24
Department of Physics, University of California, Berkeley, California 94720, USA
25
College of Natural and Health Sciences, Zayed University, Dubai, 19282, United Arab Emirates
26
Laurentian University, Department of Physics, 935 Ramsey Lake Road, Sudbury, Ontario P3E 2C6, Canada
27
Department of Physics, University of Colorado Denver, Denver, Colorado 80217, USA
28
Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA
29
Department of Electrical Engineering, University of Colorado Denver, Denver, Colorado 80217, USA
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
0031-9007
=
23
=
131(9)
=
091801(7)
091801-1
© 2023 American Physical Society
30
Department of Physics, University of South Dakota, Vermillion, South Dakota 57069, USA
31
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
32
Department of Physics, Santa Clara University, Santa Clara, California 95053, USA
33
Department of Physics, Duke University, Durham, North Carolina 27708, USA
34
Lawrence Livermore National Laboratory, Livermore, California 94550, USA
(Received 21 March 2023; revised 7 July 2023; accepted 26 July 2023; published 28 August 2023)
We measured the nuclear-recoil ionization yield in silicon with a cryogenic phonon-sensitive gram-scale
detector. Neutrons from a monoenergetic beam scatter off of the silicon nuclei at angles corresponding to
energy depositions from 4 keV down to 100 eV, the lowest energy probed so far. The results show no sign
of an ionization production threshold above 100 eV. These results call for further investigation of the
ionization yield theory and a comprehensive determination of the detector response function at energies
below the keV scale.
DOI:
10.1103/PhysRevLett.131.091801
The identity of dark matter and determination of neutrino
properties are problems at the forefront of physics beyond
the standard model. Rare event searches focused on dark
matter detection
[1
6]
or coherent elastic neutrino-nucleus
scattering (CE
ν
NS)
[7
14]
often detect products generated
by a particle interacting in a target material, requiring a
strong understanding of that material
s response to energy
depositions. Silicon is a commonly used target material.
Particle interactions with the silicon nuclei or electrons
generate free charge carriers, with nuclear recoils generat-
ing fewer charge carriers than electron recoils of the same
energy. The ratio of charge carriers produced by nuclear
and electron recoils, called the ionization yield
Y
, is crucial
to understanding the response of such detectors, and is
believed to be an intrinsic material property. Experimental
measurements of
Y
in silicon
[15
17]
for nuclear recoils
above 4 keV have been consistent with a model developed
by Lindhard
et al.
[18]
. At lower energies, typical in
low-mass dark matter or CE
ν
NS searches, measurements of
Y
indicate a significant deviation from the Lindhard model
[19
21]
. Recent modeling
[22,23]
has focused on under-
standing the origin of these deviations. Furthermore,
measurements of the low-energy yield in another com-
monly used semiconductor, germanium, have been incon-
sistent with each other
[24
26]
. These observations create
the need for more independent ionization yield measure-
ments. Here, we present the result of an ionization yield
measurement in silicon using data taken with a cryogenic
detector
[27,28]
as part of a neutron scattering experiment
called ionization measurement with phonons at cryogenic
temperatures.
We perform the measurements using the tandem accel-
erator at Triangle Universities Nuclear Laboratory. The
accelerator produces a 2.5 MHz pulsed proton beam with
0.5% duty cycle, which is directed onto a 100-nm-thick
lithium fluoride (LiF) target. The target reaction
7
Li
ð
p
;
n
Þ
7
Be produces a neutron beam with a controllable
bimodal energy after collimation
[29,30]
. The ejected
neutrons elastically scatter in the silicon detector, and
neutrons from the higher energy mode are subsequently
detected in liquid scintillator cells located at known
scattering angles corresponding to six recoil energies
between 0.1 and
3
.
9
keV
nr
. While neutrons of the lower
energy mode are also present in the beam, they will not
contribute significantly to our analysis because their energy
falls below the trigger threshold of the liquid scintillator
detectors.
By tuning the proton energy to slightly over the forward
production threshold (1.881 MeV) and selecting the for-
ward-going neutrons, we produce a neutron beam at
55.7 keV with
1
keV spread. This allows us to exploit
the 55.7 keV resonance in the neutron-silicon elastic-
scattering cross section, making the measurement robust
against small drifts in the neutron energy. At the beginning
of the experiment, the neutron energy is measured to match
the expected
56
keV via the time of flight between the
beam pickup monitor (BPM) immediately upstream from
the target and a neutron detector at 0°.
The silicon detector is a
1
-cm
2
-square and 4-mm-thick
SuperCDMS HVeV detector
[27,28]
measuring the total
phonon energy (denoted as
E
t
with unit keV
t
) generated
following particle scattering. The detector is operated at
52 mK in an adiabatic demagnetization refrigerator (ADR).
A voltage bias can be applied across the 4-mm thickness,
producing phonons from the accelerated charge carriers
through the Neganov-Trofimov-Luke effect
[31,32]
. The
total phonon energy
E
t
is
E
t
¼
E
r
þ
n
eh
eV;
where
h
n
eh
Y
E
r
ε
:
ð
1
Þ
The term
E
r
is the recoil energy,
e
the elementary electric
charge,
V
the substrate voltage bias,
n
eh
the number of
electron-hole (eh) pairs generated from this recoil,
h
n
eh
i
its
averaged number for an energy deposit of
E
r
, and
ε
the
average energy per eh in silicon. The choice of
ε
is fully
correlated with the ionization yield, and our results are
based on
ε
¼
3
.
8
eV
[33]
. By operating without voltage
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-2
bias (0 V mode) we measure the recoil energy directly,
whereas by applying a voltage bias (HV mode) our detector
becomes sensitive to the ionization signal. This work
presents an ionization yield measurement based on 11 days
of data taken with a bias of 100 V, with cross-checks
performed in 0 V mode.
We record the silicon detector data continuously at a
sampling frequency of
1
.
5
MS
=
s and use an offline trigger
to identify energy depositions. The detector achieved a
baseline energy resolution of
4
.
5
eV
t
. A trigger threshold
of
50
eV
t
is used to avoid near-threshold effects.
Calibration of the detector up to phonon energies of
120
keV
t
is accomplished with a laser and an
55
Fe source
[28]
. Laser calibration data are taken daily to allow us to
correct for gain variations caused by ADR thermal
cycles.
For the secondary neutron detection, we use a total of 29
liquid scintillator cells filled with Eljen EJ-301 or EJ-309
coupled to Hamamatsu R7724 photoelectron multiplier
tubes (PMT) arranged at angles corresponding to six
different nuclear recoil energies. Twenty-six of the neutron
detectors are mounted on two concentric rings of radii
29.4 cm and 45.2 cm, referenced as ring detectors. The
rings are centered on the beam axis and are first placed at a
distance of 86 cm downstream from the silicon detector for
seven days to perform the measurement at 0.46 and
0
.
22
keV
nr
, and then at 131 cm downstream for four days
for 0.22 and
0
.
1
keV
nr
. The remaining three neutron
detectors, referenced as the lone-wolf (LW) detectors,
are each positioned to measure recoil energies of 0.75,
2.2 and
3
.
9
keV
nr
20
cm away from the silicon detector.
The neutron detectors are calibrated daily against the
Compton edge of
137
Cs gamma rays. The scattered neutrons
deposit a maximum of
50
keV in the liquid scintillator,
equivalent to a 5 keV electron recoil (
5
keV
ee
) assuming a
quenching factor of 10%
[34]
. The variation of this
quenching factor only affects the expected efficiency for
tagging coincidence events and does not have a significant
contribution to the systematic uncertainty.
The PMT and BPM waveforms are digitized simulta-
neously at a sampling frequency of
250
MS
=
s following a
hardware trigger window of
3
keV
ee
to
40
keV
ee
on
each PMT channel. We expect negligible detection of the
low-energy mode of neutrons from the beam
[34]
.To
identify coincidence events between the silicon and the
secondary neutron detectors, we synchronize the clock of
the two data acquisition systems every minute.
To understand the expected recoil energy distribu-
tions for each secondary neutron detector, we simulate
the experiment using
G
eant4
[35]
-10.05.p01 with the
Shielding
physics list. The simulation model includes
the neutron beam collimator, the ADR, the silicon detector,
and the scintillator cells of the secondary neutron detectors.
Neutrons are emitted from a spot of 1-mm radius in a cone
of 4° half angle from the location of the target. Angular
dependence of the neutron kinematics from the emission is
negligible for a collimator with a 3° opening angle
[29]
.
The silicon detector position is determined to within
1
cm by scanning with a tightly collimated
57
Co x-ray
source on a translation stage. The detector position uncer-
tainty translates into an uncertainty on the measured recoil
energy. This uncertainty is negligible for the ring detectors
but not for the LW detectors. To account for this, we
linearly scale the simulated LW recoil energy spectra to the
LW spectra measured in 0 V mode, which gives us a direct
measurement of the recoil energies observed by each LW.
The scale factors, determined by a binned-likelihood
minimization, are reported in Table
I
. We adopt the
uncertainties on these scale factors as a systematic uncer-
tainty on the recoil energy measured by the LWs.
The initial neutron energy distribution in the simulation
is generated from a semianalytical model of the proton
beam interactions in the LiF target. Initial proton energies
are sampled from a Gaussian distribution with a standard
deviation of 2 keV
[36]
and a mean value depending on the
targeted neutron energy. The energy loss in LiF, as
determined by
TRIM
[37]
, can then be converted into a
neutron energy following the kinematics outlined in
Ref.
[30]
. We sample neutron energy distributions for
targeted nominal beam energies ranging from 46 to
60 keV in 1 keV steps.
We use these simulations to evaluate the stability of the
neutron beam energy hourly based on a fit of the spectrum
of neutron energy depositions in the silicon detector
between 4.0 and 8.2 keV. The lower bound is outside
our region of interest for the ionization yield analysis, while
the upper bound marks an energy region where neutron
interactions dominate. Since the discrepancy between the
Lindhard model and existing measurements is only 20% in
this energy range, we use it to convert the simulated recoil
energies to total phonon energies. The neutron beam
simulation matching the data best is chosen to represent
the beam behavior for that hour. This method yields a
spread in the best-fit beam energy of about 3 keV, which we
take to be the systematic uncertainty on the beam energy.
We remove silicon detector events occurring on the tails
of large energy depositions to ensure accuracy in energy
reconstructions. Multiple energy depositions occurring
close in time, referred to as pile-up events, pose challenges
in energy estimations. We allow up to two pileups per event
TABLE I. Nominal recoil energy for the LW detectors and
corrections for detector offset by using the direct measurement
in 0 V mode.
Nominal
E
r
Measured
E
r
Scale factor
0.75 keV
0.89 keV
1
.
18
þ
0
.
03
0
.
15
2.00 keV
2.33 keV
1
.
16
þ
0
.
05
0
.
14
3.87 keV
3.91 keV
1
.
01
þ
0
.
01
0
.
11
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-3
after demonstrating negligible biases in energy estimations
using the
matched-filter-integral
algorithm with built-in
pileup corrections
[28]
. For neutron detector events, we
remove pile-up events and select events below
10
keV
ee
to
reject photon events. Events where most energy is depos-
ited within a single time bin are also removed, as these are
inconsistent with the expected neutron pulse shape
[38]
.
We define a variable
dt
as the time difference between
any silicon detector event and neutron detector events. Pairs
of events in the silicon and the neutron detectors with
j
dt
dt
0
j
<
0
.
9
μ
s are referred to as coincident events, in which
dt
0
is the neutron time of flight between the two detectors.
The
0
.
9
μ
s coincidence window corresponds to 3 standard
deviations of the
dt
timing resolution (see Fig.
1
).
To suppress background from random coincidences, we
require the time of flight between BPM and neutron
detector (denoted as TOF) to be consistent with 55.7 keV
neutrons. This corresponds to a peak in the TOF distribu-
tion (see Fig.
1
). The simulated TOF is fitted to data with a
time offset to account for a small discrepancy in the
geometry. We select events with a window corresponding
to a

2
keV spread in the neutron energy. The energy
spectra of all recoil energies after event selection are shown
in Fig.
2
.
Background events with no correlation between the
neutron detector and the silicon detector form a flat
distribution in
dt
. We estimate them with sidebands from
20 to
50
μ
s before or after the coincidences, as shown in
Fig.
1
. The estimated background spectra are shown as the
blue colored component in Fig.
2
.
We build a detector response model that transforms the
simulated nuclear recoil energy spectra into the total
phonon energy spectra. We parametrize
Y
ð
E
r
Þ
, the
energy-dependent ionization yield, as a piecewise linear
function characterized by the values at the recoil energies
we measure. The function is fixed at low and high energies
to be
Y
ð
E
r
¼
0
keV
Þ¼
0
, and
Y
ð
E
r
¼
10
keV
Þ¼
0
.
3
from the Lindhard model. For each simulated event with
E
r
,
h
n
eh
i
is calculated with Eq.
(1)
. Then
n
eh
is sampled
from a distribution with this mean and a variance [
σ
2
ð
n
eh
Þ
]
characterized by a Fano factor
F
¼
σ
2
ð
n
eh
Þ
=
h
n
eh
i
.
Poission, binomial, and negative binomial distributions
are used to model
F
¼
1
,
F<
1
, and
F>
1
, respectively.
The nuclear-recoil Fano factor is unknown at this energy,
but at larger energies there is evidence that it can be
1
[39]
.
n
eh
is further smeared to account for charge trapping
(12.7% probability) and impact ionization (0.6% proba-
bility) in the detector
[28]
. An energy-dependent Gaussian
distributed detector resolution (see Ref.
[28]
) is applied to
the total phonon energy
E
t
after conversion from
n
eh
with Eq.
(1)
.
We compare the simulated
E
t
spectra to data after
applying their normalization factors, denoted as
n
. Each
simulated
E
t
spectrum has three free parameters:
Y
,
F
, and
n
. For the 220 eV recoil energy point that is measured
FIG. 1. Time difference between the silicon detector and
neutron detectors (
dt
) and time difference between BPM and
neutron detectors (TOF) for
890
eV
nr
scattering. Left:
dt
dis-
tribution. The signal window is marked in orange; the region
corresponding to random coincidences used for background
estimation is marked in blue. Right: TOF distribution. The
selection window is marked in green.
FIG. 2. Comparison of the fit model using the best-fit param-
eters to the data for all recoil energies. For the 100 eV recoil
energy (top panel), the integer number of eh contributions to the
model are shown colored for pair number up to
n
eh
¼
5
.
Individual eh pairs are not shown for other recoil energies
because they are indistinguishable based on our energy resolution
and the number of events. Blue regions are the estimated
background.
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-4
twice, we constrain the two fits to have the same yield and
the same Fano factor. We sample these parameters simul-
taneously with a Markov-Chain Monte Carlo (MCMC)
method using a binned-likelihood loss function imple-
mented in the
B
ayesian
A
nalysis
T
oolkit
[40]
.
To improve the speed of convergence, we perform the fit
in an iterative way. The first fit is for the ring detector
parameters while keeping the yield at the LW detector
energies fixed to the Izraelevitch result
[20]
. We then fit for
the LW detector parameters while keeping the ring param-
eters at the previous result. Finally, the fit is rerun on the
ring parameters with the LW parameters at the previous
result. The results with the best-fit parameters are shown
in Fig.
2
.
The systematic uncertainties are provided below.
(1) Recoil energy uncertainty has two contributions. A

1
.
3%
uncertainty arises from the total phonon energy
scale calibration of the 0 V mode data. For the LW detectors
additional uncertainties come from the scale factors in
Table
I
. (2) Regarding neutron beam energy uncertainty, the
central energy of the neutron beam as measured had a
spread of

3
keV. We vary the beam energy in the
simulation and use the resulting spectra for the TOF cut
and the fit model. (3) The charge trapping and impact
ionization uncertainties are probabilities varied by the
uncertainty from Ref.
[28]
. They are varied conservatively
such that when one probability is increased the other is
decreased. (4) With respect to the time of flight cut
uncertainty, the neutron time of flight is correlated with
their energy. The effect of the TOF selection is evaluated by
choosing a wider (narrower) window selecting events
within

50%
(

20%
) of the simulated TOF distribution
maximum. (5) We also note uncertainty in modeling the
Fano factor. Owing to poor knowledge of the Fano factor,
we perform a fit with the Fano factor fixed to one.
Each systematic uncertainty is evaluated at
1
σ
signifi-
cance by fluctuating the corresponding parameter and
performing the MCMC fit to 100 pseudoexperiments.
These pseudoexperiments are generated by applying our
detector response model to the simulated recoil energies
with the nominal fit results. Each bin in the resulting
distribution is then fluctuated by a Poisson random number,
resulting in a single pseudoexperiment. To be conservative,
we choose the deviation of the central yield value in each
scenario from the original best-fit yield plus the standard
deviation of these fits as our estimate for the given
systematic uncertainty. The systematic uncertainties are
assigned as two-sided asymmetric if the positive and
negative fluctuated samples yield two-sided deviations;
otherwise symmetric two-sided uncertainties are assigned
with the largest deviation. All calculated systematic uncer-
tainties are added in quadrature to obtain the total system-
atic uncertainty. For the LWs, the uncertainty on the recoil
energy dominates because of the large position uncertainty.
The statistical uncertainty dominates for the ring detectors.
The fit results and the uncertainties are provided in
Table
II
and Fig.
3
. The correlations among the fit
parameters are found to be negligible. We provide a
FIG. 3. The measured ionization yields, along with their
statistical and total uncertainties and a fit with a power-law
function. Also shown are data points from previous measure-
ments
[16,17,19
21,41]
. The dashed line shows the Lindhard
model with
k
¼
0
.
146
[42]
.
TABLE II. Measured silicon ionization yield, Fano factor, and signal normalization, with uncertainties. The 220 eV normalization
given without (with) parentheses is for the far (near) position of the ring detectors. The remaining columns provide the statistical (stat)
uncertainty and the systematic uncertainties from the uncertainty on the recoil energy, the neutron beam energy, the charge trapping and
impact ionization probabilities (CT
=
II), the TOF cut, and (potential) deficiencies in modeling the Fano factor.
E
r
[keV
nr
]
Ionization
yield
Y
Fano
factor Normalization Statistical
Recoil
energy
Beam
energy CT
=
II
TOF
Fano factor
mismodeling
0.10
0
.
102
þ
0
.
034
0
.
030
0
.
9
þ
0
.
7
0
.
4
28
þ
6
5
þ
0
.
024
0
.
019
þ
0
.
006
0
.
006
þ
0
.
005
0
.
004

0
.
004

0
.
002

0
.
022
0.22
0
.
108
þ
0
.
009
0
.
010
0
.
5
þ
0
.
2
0
.
1
48
þ
7
7
(
118
þ
9
9
)
þ
0
.
006
0
.
006
þ
0
.
001
0
.
002
þ
0
.
002
0
.
004

0
.
001
þ
0
.
002
0
.
001

0
.
005
0.46
0
.
136
þ
0
.
009
0
.
008
1
.
8
þ
0
.
6
0
.
5
230
þ
14
13
þ
0
.
007
0
.
006
þ
0
.
003
0
.
002

0
.
004

0
.
001

0
.
001

0
.
001
0.89
0
.
127
þ
0
.
031
0
.
015
3
.
7
þ
0
.
8
0
.
9
288
þ
11
12
þ
0
.
006
0
.
006
þ
0
.
028
0
.
006

0
.
008

0
.
006
þ
0
.
001
0
.
002

0
.
007
2.33
0
.
173
þ
0
.
044
0
.
019
7
.
7
þ
3
.
2
2
.
2
377
þ
16
14
þ
0
.
006
0
.
006
þ
0
.
042
0
.
012

0
.
008
þ
0
.
002
0
.
007
þ
0
.
003
0
.
001

0
.
010
3.91
0
.
236
þ
0
.
055
0
.
009
8
.
4
þ
2
.
4
1
.
9
318
þ
12
15
þ
0
.
005
0
.
004
þ
0
.
054
0
.
007
þ
0
.
004
0
.
002
þ
0
.
006
0
.
003

0
.
002

0
.
001
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-5
least-square fit to our results on the ring detectors with an
empirically chosen power-law function
Y
ð
E
r
Þ¼
Y
10
keV
·
ð
E
r
=
10 000
Þ
B
that is constrained to go through the yield
of the Lindhard model at 10 keV (
Y
10
keV
). The resulting
best fit yields
B
¼
0
.
261
þ
0
.
017
0
.
011
,given
Y
10
keV
¼
0
.
302
.
Our results show some tension with an earlier experi-
ment using a photoneutron source
[19]
at 890 eV. There is
also tension at lower energies with the recent result using
silicon neutron capture
[21]
, which may be caused by the
choice to fit to the Sorensen model
[22]
with a finite
ionization threshold. Our results agree with the similar
neutron-scattering setup
[20]
above 2 keV. Our measure-
ment of the ionization yield of nuclear recoils in silicon is
the first reaching down to 100 eV. The previously noted
deviation from the Lindhard model extends down to 100 eV
with no indication of an ionization production threshold.
This latter fact is of great importance to rare event search
experiments in semiconductor detectors.
The SuperCDMS Collaboration gratefully acknowledges
the Triangle Universities Nuclear Laboratory (TUNL)
facility and its staff. 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 PNNL is
operated by the Battelle Memorial Institute for the U.S.
Department of Energy under Contracts No. DE-AC02-
37407CH11359, No. DE-AC02-76SF00515, and No. DE-
AC05-76RL01830, respectively. The TUNL accelerator is
operated and maintained for the U.S. Department of Energy
with Grant No. DE-FG02-97ER41033. This research was
enabled in part by support provided by SciNet
[43]
and the
Digital Research Alliance of Canada
[44]
.
*
Corresponding author.
runzeren2023@u.northwestern.edu
[1] A. H. Abdelhameed, G. Angloher, P. Bauer, A. Bento, E.
Bertoldo
et al.
(CRESST Collaboration),
Phys. Rev. D
100
,
102002 (2019)
.
[2] E. Armengaud, C. Augier, A. Benoit, L. Berge, J. Billard
et al.
(EDELWEISS Collaboration),
Phys. Rev. D
98
,
082004 (2018)
.
[3] R. Agnese
et al.
(SuperCDMS Collaboration),
Phys. Rev. D
95
, 082002 (2017)
.
[4] H. Ma, Z. She, Z. Liu, L. Yang, Q. Yue, Z. Zeng, and T. Xue
(CDEX Collaboration),
J. Phys. Conf. Ser.
1468
, 012070
(2020)
.
[5] O. Abramoff, L. Barak, I. M. Bloch, L. Chaplinsky, M.
Crisler
et al.
(SENSEI Collaboration),
Phys. Rev. Lett.
122
,
161801 (2019)
.
[6] A. Aguilar-Arevalo, D. Amidei, D. Baxter, G. Cancelo,
B. A. Cervantes Vergara
et al.
(DAMIC Collaboration),
Phys. Rev. Lett.
125
, 241803 (2020)
.
[7] D. Akimov
et al.
(COHERENT Collaboration),
Science
357
, 1123 (2017)
.
[8] B. Cabrera-Palmer and D. Reyna, Sandia Report
No. SAND2012-8021, Sandia National Laboratories, 2012.
[9] H. T. Wong,
Nucl. Phys. B, Proc. Suppl.
138
, 333 (2005)
.
[10] V. Belov
et al.
,
J. Instrum.
10
, 12011 (2015)
.
[11] C. Buck
et al.
, J. Phys.
1342
, 012094 (2020).
[12] J. Billard
et al.
,
J. Phys. G
44
, 105101 (2017)
.
[13] G. Agnolet
et al.
,
Nucl. Instrum. Methods Phys. Res., Sect.
A
853
, 53 (2017)
.
[14] A. Aguilar-Arevalo
et al.
,
J. Phys. Conf. Ser.
761
, 012057
(2016)
.
[15] A. R. Sattler,
Phys. Rev.
138
, A1815 (1965)
.
[16] G. Gerbier, E. Lesquoy, J. Rich, M. Spiro, C. Tao
et al.
,
Phys. Rev. D
42
, 3211 (1990)
.
[17] B. L. Dougherty,
Phys. Rev. A
45
, 2104 (1992)
.
[18] J. Lindhard, V. Nielsen, M. Scharff, and P. V. Thomsen, Mat.
Fys. Medd. K. Dan. Vidensk. Selsk
33
, 1 (1963),
https://
gymarkiv.sdu.dk/MFM/kdvs/mfm%2030-39/mfm-33-10
.pdf
.
[19] A. E. Chavarria
et al.
,
Phys. Rev. D
94
, 082007 (2016)
.
[20] F. Izraelevitch
et al.
,
J. Instrum.
12
, 06014 (2017)
.
[21] A. N. Villano, M. Fritts, N. Mast, S. Brown, P. Cushman, K.
Harris, and V. Mandic,
Phys. Rev. D
105
, 083014 (2022)
.
[22] P. Sorensen,
Phys. Rev. D
91
, 083509 (2015)
.
[23] Y. Sarkis, A. Aguilar-Arevalo, and J. C. D
Olivo,
Phys. Rev.
D
101
, 102001 (2020)
.
[24] J. I. Collar, A. R. L. Kavner, and C. M. Lewis,
Phys. Rev. D
103
, 122003 (2021)
.
[25] M. Albakry, I. Alkhatib, D. Amaral, T. Aralis, T. Aramaki, I.
Arnquist, I. A. Langroudy, E. Azadbakht, S. Banik, C.
Bathurst
et al.
,
Phys. Rev. D
105
, 122002 (2022)
.
[26] A. Bonhomme, H. Bonet, C. Buck, J. Hakenmüller, G.
Heusser, T. Hugle, M. Lindner, W. Maneschg, R. Nolte, T.
Rink
et al.
,
Eur. Phys. J. C
82
, 815 (2022)
.
[27] R. K. Romani
et al.
,
Appl. Phys. Lett.
112
, 043501 (2018)
.
[28] R. Ren, C. Bathurst, Y. Y. Chang, R. Chen, C. W. Fink
et al.
,
Phys. Rev. D
104
, 032010 (2021)
.
[29] C. Lee and X.-L. Zhou,
Nucl. Instrum. Methods Phys. Res.,
Sect. B
152
, 1 (1999)
.
[30] A. O. Hanson, R. F. Taschek, and J. H. Williams,
Rev. Mod.
Phys.
21
, 635 (1949)
.
[31] B. Neganov and V. Trofimov, Otkryt. Izobret
146
,53
(1985),
https://inspirehep.net/literature/1416918
.
[32] P. N. Luke,
J. Appl. Phys.
64
, 6858 (1988)
.
[33] R. H. Pehl, F. S. Goulding, D. A. Landis, and M. Lenzlinger,
Nucl. Instrum. Methods
59
, 45 (1968)
.
[34] C. Awe, P. S. Barbeau, J. I. Collar, S. Hedges, and L. Li,
Phys. Rev. C
98
, 045802 (2018)
.
[35] R. Brun, L. Urban, F. Carminati, S. Giani, M. Maire, A.
McPherson, F. Bruyant, and G. Patrick,
G
eant
: Detector
description and simulation tool, Technical Report, CERN,
1993.
[36] TUNL (private correspondence).
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-6
[37] J. F. Ziegler, M. D. Ziegler, and J. P. Biersack,
Nucl.
Instrum. Methods Phys. Res., Sect. B
268
, 1818
(2010)
.
[38] F. T. Kuchnir and F. J. Lynch,
IEEE Trans. Nucl. Sci.
15
,
107 (1968)
.
[39] M. Matheny, A. Roberts, A. Srinivasan, and A. N. Villano,
Phys. Rev. D
106
, 123009 (2022)
.
[40] A. Caldwell, D. Kollar, and K. Kröninger,
Comput. Phys.
Commun.
180
, 2197 (2009)
.
[41] R. Agnese
et al.
(SuperCDMS Collaboration),
Nucl. Ins-
trum. Methods Phys. Res., Sect. A
905
, 71 (2018)
.
[42] J. D. Lewin and P. F. Smith,
Astropart. Phys.
6
, 87 (1996)
.
[43]
www.scinethpc.ca
.
[44]
www.alliancecan.ca
.
PHYSICAL REVIEW LETTERS
131,
091801 (2023)
091801-7