Search for low-mass dark matter via bremsstrahlung radiation
and the Migdal effect in SuperCDMS
M. F. Albakry,
1,2
I. Alkhatib,
3
D. Alonso-González,
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,
16
M. Chaudhuri,
11
R. Chen,
17
N. Chott,
18
J. Cooley,
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,
17
C. W. Fink,
16
K. Fouts,
13
M. Fritts,
20
G. Gerbier,
19
R. Germond,
19,2
M. Ghaith,
24
S. R. Golwala,
7
J. Hall,
25,26
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
D. Jardin,
15
,
†
V. K. S. Kashyap,
11
M. H. Kelsey,
10
A. Kubik,
25
N. A. Kurinsky,
13
M. Lee,
10
A. Li,
1,2
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
J. S. Mammo,
30
N. Mast,
20
A. J. Mayer,
2
H. Meyer zu Theenhausen,
22,23
É. Michaud,
21
E. Michielin,
1,2
N. Mirabolfathi,
10
B. Mohanty,
11
J. Nelson,
20
H. Neog,
20
V. Novati,
17
J. L. Orrell,
9
M. D. Osborne,
10
S. M. Oser,
1,2
W. A. Page,
16
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,
16,31
W. Rau,
2
E. Reid,
6
R. Ren,
17
T. Reynolds,
3
A. Roberts,
27
A. E. Robinson,
21
T. Saab,
12
D. Sadek,
12
B. Sadoulet,
16,31
I. Saikia,
15
J. Sander,
30
A. Sattari,
3
B. Schmidt,
17
R. W. Schnee,
18
S. Scorza,
25,26
B. Serfass,
16
S. S. Poudel,
9
D. J. Sincavage,
20
P. Sinervo,
3
J. Street,
18
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,23
S. L. Watkins,
16
O. Wen,
7
Z. Williams,
20
M. J. Wilson,
22
J. Winchell,
10
C.-P. Wu,
21
K. Wykoff,
18
S. Yellin,
14
B. A. Young,
32
T. C. Yu,
13
B. Zatschler,
3
S. Zatschler,
3
A. Zaytsev,
22,23
E. Zhang,
3
L. Zheng,
10
and A. Zuniga
3
(SuperCDMS Collaboration)
1
Department of Physics and Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada
2
TRIUMF, Vancouver, BC V6T 2A3, Canada
3
Department of Physics, University of Toronto, Toronto, ON M5S 1A7, Canada
4
Departamento de Física Teórica, 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, and 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, University of California, Berkeley, California 94720, USA
17
Department of Physics and Astronomy, Northwestern University, Evanston, Illinois 60208-3112, USA
18
Department of Physics, South Dakota School of Mines and Technology, Rapid City,
South Dakota 57701, USA
19
Department of Physics, Queen
’
s University, Kingston, ON K7L 3N6, Canada
20
School of Physics and 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
College of Natural and Health Sciences, Zayed University, Dubai, 19282, United Arab Emirates
25
SNOLAB, Creighton Mine #9, 1039 Regional Road 24, Sudbury, ON P3Y 1N2, Canada
PHYSICAL REVIEW D
107,
112013 (2023)
2470-0010
=
2023
=
107(11)
=
112013(14)
112013-1
Published by the American Physical Society
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, IL 60510, USA
29
Department of Electrical Engineering, University of Colorado Denver, Denver, Colorado 80217, USA
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
(Received 8 March 2022; accepted 12 May 2023; published 30 June 2023)
We present a new analysis of previously published SuperCDMS data using a profile likelihood
framework to search for sub-GeV dark matter (DM) particles through two inelastic scattering
channels: bremsstrahlung radiation and the Migdal effect. By considering these possible inelastic
scattering channels, experimental sensitivity can be extended to DM masses that are undetectable through
the DM-nucleon elastic scattering channel, given the energy threshold of current experiments. We exclude
DM masses down to
220
MeV
=c
2
at
2
.
7
×
10
−
30
cm
2
via the bremsstrahlung channel. The Migdal channel
search provides overall considerably more stringent limits and excludes DM masses down to
30
MeV
=c
2
at
5
.
0
×
10
−
30
cm
2
.
DOI:
10.1103/PhysRevD.107.112013
I. INTRODUCTION
An abundance of evidence suggests that most of the
Universe is composed of nonluminous matter
[1
–
3]
. This
“
dark matter
”
(DM) may consist of an undiscovered
elementary particle or a set of particles
[4]
. A leading
category of hypothesis are Weakly Interacting Massive
Particles (WIMPs). However, since particle DM has not
been detected directly, its exact properties, such as mass
and interaction cross section with standard model particles,
have yet to be determined.
Much effort has been focused on searches for particles
with masses in the GeV
=c
2
to TeV
=c
2
range, where the
favored detection mechanism is rare collisions observed by
terrestrial detectors
[5]
. Some of these approaches can be
extended to reach below
1
GeV
=c
2
through inelastic detec-
tion channels. In canonical direct DM searches, the inter-
action between a DM particle and a nucleus is assumed to be
an elastic two-body interaction. For DM particles with
masses,
m
χ
, much smaller than that of the target nucleus
m
N
, the recoil energy from an elastic collision is suppressed
by the kinematic term
m
2
χ
=m
N
, resulting in rapidly dimin-
ishing sensitivity when considering lower mass DM candi-
dates. This suppression is the result of momentum and
energy conservation with a heavy target nucleus, but it can
be circumvented by involving a third particle in the scatter-
ing process when
m
χ
≪
m
N
. In such inelastic scatterings,
the third particle can receive up to the full energy of the
collision
[6]
. Detection of this higher energy particle
provides sensitivity to DM masses that were not considered
because the energy from the elastic collision was below the
detector threshold.
The inelastic processes considered in this analysis
originate from spin-independent nuclear recoil events that
produce either a photon or an electron
[6,7]
. Therefore,
these results are directly comparable to existing limits for
DM-nucleon interactions.
In this paper, we present a reanalysis of data from the
Super Cryogenic Dark Matter Search (SuperCDMS)
experiment to look for DM scattering inelastically off of
nuclei. Section
II
describes the experiment and data
selection. Section
III
discusses how we account for the
scattering of DM through the atmosphere and Earth before
it reaches the underground experiment. Section
IV
details
the two signal models considered in this analysis. Section
V
specifies the background models included in the likelihood
framework, and Sec.
VI
describes the limit-setting method.
The final results are presented in Sec.
VII
.
II. SUPERCDMS
The SuperCDMS experiment was operated
∼
700
m
underground in the Soudan Underground Laboratory from
2012 to 2015. During this period, 15 germanium crystal
detectors were used to search for DM particle masses from
a few to tens of GeV
=c
2
[8
–
10]
. The 3-inch diameter,
1-inch-thick cylindrical detectors were shielded from
ambient radiation in the experimental cavern by layers
of polyethylene, lead, and a copper cryostat. The crystals
*
Corresponding author.
rcalkins@smu.edu
†
Corresponding author.
daniel.jardin@northwestern.edu
Published by the American Physical Society under the terms of
the
Creative Commons Attribution 4.0 International
license.
Further distribution of this work must maintain attribution to
the author(s) and the published article
’
s title, journal citation,
and DOI. Funded by SCOAP
3
.
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-2
were instrumented with interleaved
Z
-sensitive ionization
and phonon (iZIP) sensors
[11]
. The detectors were biased
face-to-face with
2
V and achieved an electron-recoil
energy threshold of about 860 eV with discrimination
between nuclear recoils (NR) and electron recoils (ER)
down to 2 keV
[12]
.
Two detectors were, for some periods, operated with a
high-voltage bias near
∼
75
V across the crystal
[9,10]
. This
mode of operation, referred to as the CDMS low ionization
threshold experiment (CDMSlite), takes advantage of the
Neganov-Trofimov-Luke (NTL) mechanism
[13,14]
to
amplify small ionization signals. The amplification lowers
the threshold of the experiment below an energy of
100
eV
ee
(ER equivalent energy)
[9,10]
but sacrifices all
discrimination between NR and ER events.
A. CDMSlite data
For this analysis, we consider the data collected by one
of the CDMSlite detectors that was operated from February
2015 to May 2015 and collected 60.9 days of raw livetime
[10]
. The exposure was divided into two segments, period 1
(P1) and period 2 (P2), due to changes in the operating
conditions and parasitic resistance that affected the actual
voltage across the crystal. These data were first analyzed in
Ref.
[10]
. Thus, the search presented in this paper was
conducted on an unblinded dataset.
B. Event selection
Since the energy region of interest,
70
eV
ee
to
30
keV
ee
,
for this analysis largely overlaps with Ref.
[10]
, the same
data quality selection criteria were used to remove prob-
lematic events such as those arising from low-frequency
mechanical noise, electronic glitches, and poorly recon-
structed pulse shapes.
The grounded copper housing surrounding the detector
distorts the electric field near the edges of the crystal. Since
the electric potential is reduced in these regions, the ampli-
fication is not uniform throughout the crystal. To minimize
the number of events with reduced amplification, the same
fiducialvolumeselectionasdefinedinRef.
[10]
wasadopted,
which rejects events with low NTL amplification. After
applying the selection criteria, the remaining exposure is
36
.
9
kg · d, and the analysis threshold is
70
eV
ee
[10]
.
III. DAMPED VELOCITY DISTRIBUTION
Calculation of the DM flux requires knowledge of the
velocity distribution of the incoming DM. The standard
halo model (SHM) was assumed, which is based on a
Maxwell-Boltzmann velocity distribution with a character-
istic velocity of
220
km
=
s
[15]
. Particles traveling at
velocities greater than the escape velocity of the galaxy
are not gravitationally bound, so the distribution is trun-
cated at
544
km
=
s and renormalized
[16]
. The local DM
density,
ρ
χ
, is assumed to be
0
.
3
GeV
=
ð
c
2
cm
3
Þ
[17]
.
The Earth is typically considered to be transparent to
DM. However, for large coupling strengths between DM
and nuclei, on the order of
∼
10
−
30
cm
2
, the Earth and
atmosphere can no longer be neglected
[18
–
20]
.AsDM
particles travel through the Earth and atmosphere, they can
scatter off of atoms and lose energy. In the most extreme
case, the DM can lose enough energy so that it can no
longer create a signal above the detector threshold. This
Earth shielding limits the sensitivity of experiments since
DM with stronger couplings is fully attenuated before
reaching the detector
[21]
.
This analysis accounts for the attenuation effect by
damping the DM velocity distribution, described as
“
Method B
”
in Ref.
[22,23]
. This approach allows for
DM with cross sections in an intermediate region to scatter
and lose some energy yet still reach the detector. It is also
flexible enough to use a more complex shielding model,
which is detailed in Sec.
III A
.
Interactions with normal matter alter the velocity of DM
particles at the detector by shifting the distribution to lower
values. The DM velocity distribution at the detector site is
calculated from the average energy loss via scattering off of
nuclei, as described in Ref.
[22,24]
. Since we concern
ourselves with light dark matter in this paper, we assume a
nuclear form factor of unity. Avelocity damping parameter,
κ
, is defined as:
κ
≡
ρσ
n
m
χ
μ
2
n
X
elements
i
F
i
μ
4
i
A
2
i
m
2
i
d;
ð
1
Þ
with the variables defined in Table
I [22]
. The calculation of
κ
is location specific since it depends on the path length,
d
,
and the type of material the particle travels through before
reaching the detector.
An incoming DM particle with initial velocity
v
i
reaches
the detector with velocity
v
f
,
v
f
¼
v
i
e
−
κ
;
ð
2
Þ
TABLE I. Definition of variables in Eq.
(1)
to calculate the
velocity damping parameter,
κ
.
Variable
Definition
m
χ
Dark matter mass
ρ
Material density
σ
n
DM-nucleon scattering cross section
μ
n
Reduced mass (DM particle, nucleon)
i
Element index
F
i
Mass fraction of the element
i
μ
i
Reduced mass (DM particle, nucleus of element
i
)
A
i
Atomic number of element
i
m
i
Atomic mass of element
i
d
Path length through the shielding layer
SEARCH FOR LOW-MASS DARK MATTER VIA
...
PHYS. REV. D
107,
112013 (2023)
112013-3
thus, the DM velocity distribution at the detector (
f
d
)is
given by:
f
d
ð
v
f
Þ¼
e
2
κ
×
f
ð
v
f
e
κ
Þ
;
ð
3
Þ
where the incoming SHM velocity distribution (
f
) has been
transformed by the effects of shielding
[22]
. The exponen-
tial term in front of the distribution accounts for the
normalization of the increased flux caused by the attenu-
ated velocity distribution.
Implicit in the definition of
κ
is the dependence on the
incident angle of the DM with respect to the detector. For
example, particles originating from directly above the
experiment pass through the atmosphere, the local over-
burden, and the experimental shielding. Meanwhile, a
particle from below the experiment will traverse the internal
structure of the Earth instead of the local overburden.
Therefore, the incoming DM flux must be evaluated for
every angle. More information about the angular depend-
ence can be found in Sec.
III B
. We assume that the
particles travel in straight line trajectories even though
scattering will affect their trajectories. According to
Ref.
[22]
, this approach still underestimates the number
of dark matter particles with sufficient energy to interact
with the detector.
A. Earth shielding
Due to the exponential nature of Eq.
(3)
,
κ
can be
calculated for each layer of shielding independently
using Eq.
(1)
and summed to derive a cumulative value.
Four categories of shielding were considered for the
SuperCDMS experiment at Soudan: the atmosphere, the
overhead rock, the experimental shielding, and the bulk of
the Earth below the experiment.
The density of the Earth
’
s atmosphere decreases con-
tinuously as a function of altitude. However, for the
purpose of simplifying the model, a seven layer atmosphere
model based on the 1976 U.S. Standard Atmosphere was
used
[25]
. The value used for the density of each layer
corresponds to the lowest altitude of that layer, thus
overestimating the amount of shielding. Densities and
altitudes from sea level are listed in Table
II
.
The SuperCDMS experiment at Soudan was located
under 714 m of Ely greenstone and iron ore in northern
Minnesota. We consulted with a geologist at the University
of Minnesota to obtain rock and chemical compositions of
the Soudan region
[26]
. Data were provided for sectors in
eight geographic directions between radii of 100, 500,
1000, 5,000, 10,000, 20,000, and 50,000 meters.
1
The
composition and density of rock depends on the direction.
To simplify the calculation and ensure that the amount of
shielding is not underestimated, we select the direction that
gives the maximum value of
κ
for each mass and cross
section considered.
The dominant component of the shielding is the Earth
below the experiment.The conventionalmodelofthe Earthis
used, consisting of four concentric spheres: crust, mantle,
outer core, and inner core. Details of the parameters defining
the thickness and composition of each layer in the model are
available in the Appendix and in Supplemental Material
[27]
.
The smallest contribution to the damping model comes
from the shielding around the experiment itself, which was
approximated with concentric spheres. The outermost
layer is 50 cm of polyethylene (C
2
H
4
) with a density of
0
.
94
g
=
cm
3
, followed by 22.5 cm of Pb and 3 cm of Cu
[28]
. The contribution to the velocity damping parameter
from the shielding around the experiment is 0.2% for a
1
GeV
=c
2
DM particle with a nucleon scattering cross
section of
10
−
30
cm
2
and traveling straight downward,
where the Earth
’
s shielding contribution is the weakest.
B. Angular dependence
The depth of the experiment in the Earth
’
s crust produces
an asymmetric angular distribution of shielding
[29]
.
Particles originating from below the experiment must
traverse the majority of the Earth
’
s diameter, while particles
originating from directly above are only affected by the
local overburden at the Soudan Mine. We assume an
isotropic distribution of dark matter hitting the Earth.
Since the Soudan Mine is located in the northern hemi-
sphere, the majority of dark matter will come from above,
which means this assumption will conservatively under-
estimate the dark matter flux at the experiment site.
This analysis follows the angular convention in Ref.
[30]
,
which defines
θ
as the incident angle between the incoming
DM particle
’
s velocity and zenith at the experiment. By
this definition,
θ
¼
0
corresponds to particles originating
from directly below the experiment. The total path length is
calculated as:
l
¼
r
d
cos
θ
þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
r
2
E
−
ð
r
d
sin
θ
Þ
2
q
;
ð
4
Þ
TABLE II. Densities and heights of the layers in the atmosphere
model
[25]
.
Layer
Height range [km]
Maximum density [g
=
m
3
]
10
–
11
1225
211
–
20
363.91
320
–
32
88.03
432
–
47
13.22
547
–
51
1.43
651
–
71
0.86
771
–
100
0.064
1
Geological data for the Soudan region is provided in two
auxiliary files: The
SoudanRegion.csv
file contains the area
and fraction of each rock type, and the elemental mass fractions for
the chemical composition are described in
RockChem.csv
.
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-4
where
r
d
is the distance from the center of the Earth to the
SuperCDMSexperiment. Theaveragevalue of6471km was
used for the Earth
’
s radius,
r
E
. The total path length can be
generalized to determine
d
, the path length through indi-
vidual layers, used in Eq.
(1)
.
Since the composition of the Soudan geology is well
understood, the overburden above the experiment was
modeled using the local geometry data, and the shielding
below the experiment was modeled according to the
conventional Earth model. This introduced an angular
dependence on the calculation of
κ
that affects the velocity
distribution at the detector. The signal model was integrated
over the incoming angle,
θ
, at 20 discrete points sampled
uniformly in cos
θ
as indicated in Fig.
1
. We ignore the
relationship between the WIMP wind and the Earth
reference frame and calculate a single velocity distribution,
which is attenuated according to its path through the
shielding.
IV. INELASTIC SCATTERING SIGNALS
In this paper, we report on the search for a signal from
DM interactions through two inelastic scattering channels.
The first search channel is through bremsstrahlung radia-
tion, where a photon is produced during the DM-nucleon
scattering process. The second search channel is induced by
the Migdal effect, where a low-energy NR perturbs the
atomic electron cloud, occasionally emitting electrons and/
or photons.
A. Bremsstrahlung radiation
The differential scattering rate for emitting a photon of
energy
E
γ
through the bremsstrahlung process has been
derived by Kouvaris and Pradler in Ref.
[6]
:
d
σ
dE
γ
¼
4
α
j
f
ð
E
γ
Þj
2
3
π
E
γ
μ
2
N
v
2
σ
N
m
2
N
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1
−
2
E
γ
μ
N
v
2
s
1
−
E
γ
μ
N
v
2
;
ð
5
Þ
where
α
is the fine structure constant,
f
is the atomic
scattering function discussed in Sec.
IVA 1
,
μ
N
is the
reduced mass of the DM-nucleus system,
m
N
is the mass of
the nucleus,
v
is the velocity of the incoming DM particle
relative to the detector, and
σ
N
is the interaction cross
section between the DM particle and the nucleus. A signal
spectrum is obtained by integrating over the velocity
distribution at the detector while accounting for the angular
dependence and the flux,
dR
dE
γ
¼
N
T
ρ
χ
m
χ
Z
j
⃗
v
j
≥
v
min
d
3
⃗
vvf
d
ð
⃗
v
Þ
d
σ
dE
γ
;
ð
6
Þ
where
N
T
is Avogadro
’
s number divided by the atomic
mass. The minimum velocity required to induce a recoil of
E
γ
,
v
min
,isgivenby
v
min
¼
ffiffiffiffiffiffiffiffiffiffiffiffiffi
m
N
E
R
2
μ
2
N
s
þ
E
γ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
m
N
E
R
p
;
ð
7
Þ
where
E
R
is the nuclear recoil energy.
1. Photoelectric absorption cross section uncertainty
The atomic scattering function,
f
, in Eq.
(5)
is the sum of
a real and complex portion,
j
f
j
2
¼j
f
1
þ
if
2
j
2
¼
f
2
1
þ
f
2
2
:
ð
8
Þ
The components,
f
1
and
f
2
, are defined in Ref.
[31]
as:
f
1
¼
Z
þ
1
π
r
e
hc
Z
∞
0
ε
2
σ
a
ð
ε
Þ
E
2
γ
−
ε
2
d
ε
;
f
2
¼
σ
a
2
r
e
λ
;
ð
9
Þ
for the nuclear contribution, where
r
e
is the classical radius
of the electron,
σ
a
is the photoelectric absorption cross
section,
λ
is the wavelength of the emitted photon, and
ε
is a
variable of integration.
Z
is defined as
Z
−
ð
Z=
82
.
5
Þ
2
.
37
,
where
Z
is the proton number.
Measurements of
σ
a
at low energies have a range of
values, which lead to significant systematic uncertainties of
its value
[32]
. Both
f
1
and
f
2
depend on
σ
a
, so this
uncertainty enters into the calculation of the expected event
rate. Based on existing literature
[33
–
47]
, a nominal value
and uncertainty band for
σ
a
were derived as a function of
energy. Using the resulting values of
σ
a
,
j
f
j
2
was calculated
and compared in Fig.
2
, along with the commonly
referenced Henke dataset
[33]
. The range of variation is
on the order of 30% and is most prominent at low energies.
The signal calculation uses the lower bound from all
FIG. 1. Value of
κ
as a function of incident angle,
θ
, for a DM
particle with a mass of
1
GeV
=c
2
and a scattering cross section of
10
−
30
cm
2
. The curve indicates the distribution, and the dots
indicate values at which the curve was sampled to calculate the
signal models. The discontinuity below
π
=
2
indicates the location
of the transition between the core and the mantle, while the more
detailed overburden model leads to a smooth curve above
π
=
2
.
SEARCH FOR LOW-MASS DARK MATTER VIA
...
PHYS. REV. D
107,
112013 (2023)
112013-5
calculations of
j
f
j
2
, which conservatively predicts a lower
expected signal rate and thus results in a weaker limit.
B. Migdal effect
When a NR occurs, there is a delay between the initial
recoil and the response of the surrounding electron cloud,
effectively boosting the entire electron cloud simultane-
ously with respect to the nucleus
[48]
. The displacement of
the nucleus due to a DM scattering event dramatically
changes the wave function of the electrons in the surround-
ing electron cloud. As the electron cloud relaxes back to the
ground state, an electron can transition to a free state (i.e.,
be ejected), a process known as the Migdal effect
[48]
. The
formulation of this effect as applied to DM direct detection
has been calculated by Ibe, Nakano, Shoji, and Suzuki
[7]
.
The Ibe
et al.
formulation was then applied to data by
Dolan, Kahlhoefer, and McCabe
[49]
. An alternative
approach utilizing the photoelectric cross section has been
developed by Liu, Wu, Chi, and Chen
[50]
, and calcu-
lations extending to higher nuclear recoil velocities have
been performed by Cox
et al.
[51]
. In this paper, we adapt
the Ibe
et al.
formalism to be consistent with other results
within the community.
The differential rate for this process can be expressed as:
d
3
R
ion
d
E
R
d
E
e
d
v
¼
d
2
R
nr
d
E
R
d
v
×
X
nl
1
2
π
d
p
c
q
e
ð
nl
→
E
e
Þ
d
E
e
;
ð
10
Þ
where
p
c
q
e
ð
nl
→
E
e
Þ
is the transition probability for an
electron, with momentum
q
e
¼
m
e
v
nucleus
with respect to a
stationary nucleus, to be ejected from quantum state
nl
with
energy
E
e
, and d
2
R
nr
=
ð
d
E
R
d
v
Þ
describes the incoming DM
flux and scattering rate,
d
2
R
nr
d
E
R
d
v
¼
ρ
χ
σ
N
2
μ
2
N
m
χ
f
d
ð
v
Þ
v
:
ð
11
Þ
The factor of
1
=
2
π
in Eq.
(10)
is a normalization constant
and is consistent with the formulation of
[7]
. The transition
probability tables that were provided by Ref.
[7]
have been
evaluated at a reference velocity,
v
ref
¼
10
−
3
c
. Conversion
between the reference value and an arbitrary electron
momentum is given by:
p
c
q
e
ð
nl
→
E
e
Þ¼
q
e
v
ref
m
e
2
p
c
v
ref
ð
nl
→
E
e
Þ
:
ð
12
Þ
These tables were calculated assuming an isolated atom.
In a crystal detector, atoms are not isolated, and the electron
clouds are subject to interactions with atoms in nearby
lattice sites, in particular, the outer electron shells. We
assumed crystal effects are negligible for the inner shells
and exclude the outermost germanium shell (
n
¼
4
) from
the analysis
[52
–
54]
. The inner electron shells are assumed
to be dominated by interactions with the nucleus and are
sufficiently representative of an isolated atom. Excluding
the valence shell causes a minimal decrease in the expected
signal rate because most of the signal originates from the
inner shells; the majority of the valence shell contribution is
below the experimental threshold. Excluding part of the
signal model results in a more conservative estimate of the
expected rate.
Figure
3
shows the expected signal rates for an incident
DM particle with
0
.
5
GeV
=c
2
mass and
3
×
10
−
37
cm
2
FIG. 2. The atomic scattering function,
j
f
j
2
, calculated as a
function of energy using different values of the photoelectric
absorption cross section,
σ
a
, in germanium. The nominal value of
σ
a
obtained from literature was used to calculate the solid black
curve. The range of uncertainty on
j
f
j
2
is indicated by the shaded
region.
FIG. 3. Comparison of the expected DM signal rates in
germanium. For the DM mass chosen,
0
.
5
GeV
=c
2
, the elastic
NR signal is below the analysis threshold of 70 eV. The Migdal
signal extends to higher energies, but at a smaller rate than the NR
signal. The Migdal signal has a sharp cutoff at low energies
because the valence shell was not included in the analysis. The
bremsstrahlung signal extends to the same energy as the Migdal,
but with a smaller expected rate. The bump above 1 keV is where
the
n
¼
2
electron shell starts contributing to the signal rate.
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-6
cross section. This mass was chosen to highlight the sub-
GeV reach of the inelastic channels where the elastic NR
signal is below the energy threshold of the analysis and thus
undetectable. The chosen cross section is small enough that
the Earth
’
s shielding is a negligible effect.
V. BACKGROUND MODELING AND
SYSTEMATICS
In order to perform a profile likelihood analysis, all
background sources must be well understood and modeled.
The background sources considered in this analysis are
neutron activation by the
252
Cf calibration source, cosmo-
genic activation of the crystal, radiogenic Compton scatter-
ing, and
210
Pb surface contamination. Inelastic scattering of
low-energy neutrons was also considered as a background
source, but it was determined via simulation that the
background contribution was
≪
1
event for the exposure,
and therefore negligible.
A. Bulk background models
The dominant background originates from activation of
the germanium crystal via the neutron calibration source.
When stable
70
Ge in the crystal captures a neutron, it
becomes unstable
71
Ge that decays via electron capture to
71
Ga. The electron capture creates a cascade of energy in
the form of x-rays and Auger electrons, where the energy
emitted from the decay is equal to the binding energy of the
shell that captured the electron. The
K
shell has the highest
probability of capture at 87.57% and results in a 10.37 keV
emission. Each successive shell has a lower probability and
emits less energy: The
L
shell captures an electron 10.53%
of the time and emits 1.3 keV, and the
M
shell captures an
electron 1.78% of the time and emits 160 eV
[55]
. This
background was modeled by a Gaussian distribution
centered on each electron shell peak energy. The ampli-
tudes of these peaks were set relative to the
K
shell peak as
determined by the capture probabilities of each shell. There
is one overall normalization parameter for the
71
Ge back-
ground. Similarly,
68
Ge decays to
68
Ga, which can decay
through electron capture or through beta decay. The
68
Ga
beta end points are substantially higher than the energy
range considered in this analysis so its contribution is
neglected since we expect just 0.001 events.
Before the detectors are brought underground, cosmic-
ray spallation can knock nucleons out of the germanium
atoms in the crystal and create radioisotopes. One of the
more problematic cosmogenic isotopes is tritium (
3
H),
which provides a persistent source of betas due to its
half-life of 12.32 years. The tritium background was
modeled by a standard beta emission spectrum with an
end point energy of 18.6 keV and an unconstrained
normalization. The resulting normalization from the fit
was
50
20
events, which is consistent with a previous
dedicated analysis
[56,57]
.
Tritium decays dominate cosmogenic background
rates, but spallation can also leave other unstable nuclei.
The other residual nuclei were considered background
sources if they have a half-life that is long enough that
they will have not decayed away before data taking began
but also short enough that the activity is comparable to
other background rates. The additional isotopes modeled in
this analysis were
68
Ga,
65
Zn, and
55
Fe
[57]
. Other isotopes
considered are
57
Co,
54
Mn, and
49
Vn, but the expected
contribution of each was determined to be
<
1
event for the
given exposure and are neglected. Each of the modeled
isotopes decays via electron capture, like the activated
germanium, but at different energies. Contributions from K,
L, and
M
shells were modeled by Gaussian distributions
with fixed relative amplitudes with a single normalization
parameter, analogous to how
71
Ge was treated above.
All of the radioisotopes, created cosmogenically or by
source activation, described in this subsection are distrib-
uted nearly uniformly throughout the detector volume.
2
Therefore, the efficiency of the physics selection criteria
that were developed for a uniform DM signal could also be
applied to the modeling of the backgrounds originating
from those isotopes.
Gamma rays emitted by long-lived naturally occurring
unstable radioisotopes typically have energies much greater
than those of interest for this analysis, but high-energy
photons can undergo Compton scattering. To first order,
this creates a flat background continuum throughout the
analysis energy region, although the model also includes
low-energy
“
steps
”
at the electron binding energies. These
steps occur when the energy deposited by a scattering
photon has enough energy to overcome the binding energy
of a particular electron shell. As the amount of energy
increases, the number of available electrons to scatter
against increases and so does the corresponding interaction
rate. This Compton background spectrum was modeled as a
flat contribution with an error function at each shell energy,
with a width corresponding to the detector resolution, to
model the steps
[10,58]
. The relative step amplitudes were
determined from an independent fit to simulation data, and
the overall normalization was allowed to float in the
likelihood function
[10]
.
B. Surface background models
Radon daughters plating out on the detectors or sur-
rounding copper housing were treated differently than the
bulk contamination described in Sec.
VA
. Although decays
can implant radon daughters below the surface, they are
predominantly classified as surface events. To model the
expected experimental signature of these surface events, a
Geant4
[59
–
61]
simulation of
210
Pb surface contamination
on a tower of six germanium detectors was performed.
2
Studies based on simulation have shown slightly more bias
toward the surface of the crystals due to self-shielding effects.
SEARCH FOR LOW-MASS DARK MATTER VIA
...
PHYS. REV. D
107,
112013 (2023)
112013-7
The simulation allowed for the subsequent alpha decays to
implant the long-lived
210
Pb and mimic the physical radio
contamination
[58]
. The simulation indicated that surface
events in the detectors originated from three locations with
direct line of sight: the top lid (TL) of the copper housing
that directly faces the top detector in the tower, the sidewall
housing (SH) around the outer radial wall of each detector,
and the surfaces of the germanium crystal (GC) from both
the detector itself and the face of the adjacent detector.
The simulated spectra were normalized using an inde-
pendent measurement of the rate of alpha events in each of
the detectors
[58]
. For the energy range used in this
analysis, the normalization to the number of expected
events from each contribution were determined to be
N
TL
¼
158
.
4
16
.
6
,
N
SH
¼
22
.
3
1
.
4
, and
N
GC
¼
23
.
5
5
.
9
events, which contribute to the total number
of background events in the likelihood. The housing copper
is not as radio pure as the crystals so it dominates the
contribution. Events originating from the bottom tower lid
are shielded by the other detectors in the housing stack.
Correlations between normalization and shape uncertain-
ties are taken into account using morphing parameters as
described in Ref.
[10]
.
C. Efficiency model
All of the background and signal spectra were convolved
with the overall efficiency of the data selection criteria, the
details of which can be found in Ref.
[10]
. The trigger
efficiency is 90% at
70
eV
ee
and 100% above
90
eV
ee
. The
data quality selection criteria efficiency approaches 100%
around
150
eV
ee
. The largest loss of efficiency at low
energies is due to the fiducial volume selection, which
passes roughly 60% of the events above
200
eV
ee
and has
almost zero efficiency at
70
eV
ee
.
Figure
4
shows the efficiency curve for each period
used in this analysis. The notable fluctuations in the
efficiency curve below
2
keV
ee
arise from the fiducial
volume selection, which was calculated using simulated
data. The simulated data was produced using actual noise
pulses, which means that the simulation dataset statistics
are limited by the statistics of the dataset. We take the
efficiency curve with its statistical fluctuations. The effi-
ciency models have been extended to
25
keV
ee
following
the procedure in Ref.
[57]
.
The uncertainty on the efficiency curve was incorpo-
rated into the likelihood function via implementation of
nuisance parameters in the maximum likelihood fit, one
for each data period. The nuisance parameter is not a
simple normalization factor, but a morphing parameter
that allows for correlated variation between shape and
normalization of the efficiency function. This is accom-
plished by constructing one-sigma bands around the
median parametrized as a Gaussian distribution following
the same prescription as the surface background models
in Ref.
[10]
.
D. Resolution model
The background model and signal models were con-
volved with the energy resolution of the detector; thus, a
resolution model was included in the likelihood function as
additional nuisance parameters. The functional form of the
detector energy resolution is
σ
2
ð
E
Þ¼
σ
2
E
þ
BE
þð
AE
Þ
2
;
ð
13
Þ
where
σ
is the resolution,
E
is the total phonon energy,
σ
E
is the baseline noise resolution that originates from the
readout electronics,
B
is a variance that scales with energy,
and
A
accounts for any effects that scale with energy such
as pulse shape variation and position dependence
[9]
. The
parameters used are described in Sec. II. D of Ref.
[10]
.At
100
eV
ee
, the resolution of period 1 and 2 is
14
eV
ee
and
16
eV
ee
, respectively. Near the upper end of the analysis
range at
25
keV
ee
, the 1 sigma resolution is
193
eV
ee
and
198
eV
ee
for periods 1 and 2, respectively.
E. Nuclear recoil ionization yield
The energy spectrum measured in the detector is a
combination of the ER and NR components. In calculating
the expected signal rate, an assumption about the nuclear
recoil ionization yield,
Y
ð
E
NR
Þ
, is needed. We adopt the
Lindhard model [Eq.
(14)
],
Y
Lindhard
ð
E
NR
Þ¼
k
·
g
ð
ε
Þ
=
ð
1
þ
k
·
g
ð
ε
ÞÞ
;
ð
14
Þ
where
g
ð
ε
Þ¼
3
ε
0
.
15
þ
0
.
7
ε
0
.
6
þ
ε
,
ε
¼
11
.
5
E
NR
ð
keV
Þ
Z
−
7
.
3
,
Z
is the atomic number, and
k
is the free electron energy
loss
[9,62]
.
There is evidence of deviations from the Lindhard
model at low energies
[63,64]
. Therefore, the ionization
FIG. 4. Efficiency of the data selection criteria in Ref.
[10]
for
both CDMSlite data periods. The black curves are the median
efficiencies, and the red band indicates the
1
σ
uncertainty band.
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-8
production was cut off at 22.7 eV, according to an
independent measurement of the defect energy creation
threshold
[65]
. This is implemented with a hyperbolic
tangent function,
Y
ð
E
NR
Þ¼
1
2
Y
Lindhard
ð
E
NR
Þ
×
tanh
E
NR
−
22
.
7
eV
22
.
32
eV
þ
1
;
ð
15
Þ
which has a width of 22.32 eV, which was determined by
requiring that the yield nearly vanish (
Y
¼
0
.
001
) near the
band gap energy of 0.74 eV.
Systematic uncertainties on the Lindhard model are
propagated through its uncertainties in the
k
parameter.
For germanium, a nominal value of
k
¼
0
.
157
is assumed,
with a Gaussian uncertainty of
σ
¼
0
.
05
[10]
. The signal
model is calculated for the nominal,
1
σ
, and
3
σ
values
of
k
and intermediate values of
σ
are interpolated.
F. Analysis energy range
The normalization of the models was determined by the
fit to data, as described in Sec.
VI
. An example of fitting
the data is shown in Fig.
5
, where the background
models and Migdal signal model, for a WIMP with a
mass of
0
.
5
GeV
=c
2
and a cross section of
3
×
10
−
37
cm
2
,
have been scaled by the efficiency and convolved with the
resolution model.
The shape of the background models dictated the energy
range used in this analysis. A strong degeneracy between
unconstrained models that are flat or nearly flat resulted in
lack of fit convergence when maximizing the likelihood.
This degeneracy is broken by extending the analysis region
to
25
keV
ee
past the energy region considered in Ref.
[10]
.
In contrast, the surface background rates were constrained
by the independent measurement of alpha rates through
Gaussian constraints in the likelihood function, and the
cosmogenic activation lines were modeled with Gaussian
distributions that are not degenerate with flat background
models.
VI. PROFILE LIKELIHOOD ANALYSIS
An unbinned profile likelihood function was utilized
for this analysis because it provides the ability to quantify
an excess of events above the expected background
and potentially claim a DM discovery. In the absence
of an excess, a likelihood places a more constraining
exclusion limit than the optimum interval
[66,67]
tech-
nique because it incorporates knowledge of the back-
ground. It also provides a rigorous and convenient way to
account for systematic uncertainties in signal and back-
ground model parameters.
A. Likelihood function
Table
III
contains a list of the variables used in the
likelihood function. The likelihood function is composed of
FIG. 5. Example of an energy spectrum from the maximum
likelihood fit for a Migdal signal model for a WIMP with a mass
of
0
.
5
GeV
=c
2
and a cross section of
3
×
10
−
37
cm
2
(black
dashed curve). The data (blue histogram) have been logarithmi-
cally binned and overlaid with the background models (colored
solid curves). The thick black line is the sum of all the models,
signal and background. Normalization of the surface background
model components (TL, SG and GC) are described in Sec.
VB
.
The plot on the bottom shows the residual between data and
the model with the
1
σ
statistical uncertainty indicated by the
shaded region.
TABLE III. Definition of variables used in the likelihood
function (Eq.
(16)
), sorted by their given state in the fit.
Identifiers are used to label variables that are specific to a data
period, iterators are used to distinguish items from a given set,
free/constrained indicates the state in the fitting procedure,
constants are inputs to the likelihood that were calculated in
advance, and the signal and background models are functions of
energy.
Variable
Definition
State
P
1
Data period 1
Identifier
P
2
Data period 2
Identifier
n
Iterator over
N
events
Iterator
bb
Bulk background iterator
Iterator
sb
Surface background iterator
Iterator
i
,
j
General nuisance iterators
Iterator
ν
χ
Number of signal events
Free
ν
bb
Number of events in
bb
Free
ν
sb
Number of events in
sb
Constrained
(Table continued)
SEARCH FOR LOW-MASS DARK MATTER VIA
...
PHYS. REV. D
107,
112013 (2023)
112013-9
three types of terms. The first type is the overall normali-
zation term, (
L
Poiss
), which allows the total number of fitted
events to fluctuate around the number of events in the
dataset (
N
). The second type (
L
P
1
χ
;bb;sb
,
L
P
2
χ
;bb;sb
) is the core
of the likelihood that uses the signal and background PDFs
to determine the most likely signal and background rates.
The third type constrains the nuisance parameters using
auxiliary measurements. In this analysis, there are six
terms: the morphed surface backgrounds (
L
surf
), the
morphed efficiency (
L
P
1
eff
,
L
P
2
eff
), the resolution model
(
L
P
1
res
,
L
P
2
res
), and the yield model (
L
yield
). The constraint
terms are either a univariate Gaussian PDF, in the case
of the efficiency and yield, or multivariate Gaussian
distributions for the surface backgrounds and resolution.
In order to implement a multivariate constraint, a covari-
ance matrix is calculated from the normalization uncer-
tainty of the individual components and the correlations
between them
[68]
.
This analysis used the MINUIT algorithm
[69]
via the
iminuit
[70]
Python interface to maximize the log-like-
lihood function when evaluating the test statistic, which
will be defined in Eq.
(17)
.
The full likelihood function used in this analysis is given
by Eq.
(16)
, which is a single likelihood encompassing both
data periods.
B. Limit calculation
In a typical DM search, the number of signal events is
directly proportional to the interaction cross section.
Therefore, a constraint on the normalization of the signal
model can be directly converted to a single cross section
because there is a one-to-one mapping between cross
section and number of signal events. In these searches,
the shape of the signal model is determined by the mass of
the DM particle since the velocity distribution is considered
unchanged from space compared to the location of the
experiment on Earth.
This assumption does not apply for any analysis
that involves DM models with potentially large cross
sections
ln
L
¼
ln
ð
L
Poiss
Þþ
ln
ð
L
P
1
χ
;bb;sb
Þþ
ln
ð
L
P
2
χ
;bb;sb
Þþ
ln
ð
L
surf
Þþ
ln
ð
L
P
1
eff
Þþ
ln
ð
L
P
2
eff
Þþ
ln
ð
L
P
1
res
Þþ
ln
ð
L
P
2
res
Þþ
ln
ð
L
yield
Þ
¼
−
ν
χ
þ
X
6
bb
¼
1
ν
bb
þ
X
3
sb
¼
1
ν
sb
þ
X
N
P
1
n
¼
1
ln
ν
χ
f
P
1
χ
ð
E
n
Þþ
X
6
bb
¼
1
ν
bb
f
P
1
bb
ð
E
n
Þþ
X
3
sb
¼
1
ρ
P
1
sb
ð
E
n
Þ
þ
X
N
P
2
n
¼
1
ln
ν
χ
f
P
2
χ
ð
E
n
Þþ
X
6
bb
¼
1
ν
bb
f
P
2
bb
ð
E
n
Þþ
X
3
sb
¼
1
ρ
P
2
sb
ð
E
n
Þ
−
1
2
X
3
i;j
¼
1
h
ð
s
i
−
μ
s
i
Þ
T
ð
S
i;j
Þ
−
1
ð
s
j
−
μ
s
j
Þ
i
−
ð
Ξ
P
1
−
μ
P
1
Ξ
Þ
2
2
ð
σ
P
1
Ξ
Þ
2
−
ð
Ξ
P
2
−
μ
P
2
Ξ
Þ
2
2
ð
σ
P
2
Ξ
Þ
2
−
1
2
X
3
i;j
¼
1
h
ð
r
P
1
i
−
μ
P
1
r
i
Þ
T
ð
R
P
1
i;j
Þ
−
1
ð
r
P
1
j
−
μ
P
1
r
j
Þ
i
−
1
2
X
3
i;j
¼
1
h
ð
r
P
2
i
−
μ
P
2
r
i
Þ
T
ð
R
P
2
i;j
Þ
−
1
ð
r
P
2
j
−
μ
P
2
r
j
Þ
i
−
1
2
θ
2
Y
;
ð
16
Þ
due to the shielding discussed in Sec.
III
. According
to Eq.
(3)
, the effects of the shielding shift the
velocity distribution of DM particles at the detector site
to lower values. The moderated velocity distribution
affects the expected DM-nucleon scattering rate and
the shape of the recoil spectrum. This results in the
shielding parameter, defined by Eq.
(1)
and the shape of
theexpectedsignaltodependontheDMmassandcross
section.
In order to test DM hypotheses in the mass and cross
section parameter space, a test statistic based on the profile
likelihood ratio was defined as:
TABLE III.
(Continued)
Variable
Definition
State
s
Surface bg morphing parameter
Constrained
Ξ
Efficiency morphing parameter
Constrained
r
Resolution nuisance parameter
Constrained
θ
Y
Yield nuisance parameter
Constrained
N
number of events in data
Constant
E
n
Energy of event
n
Constant
μ
s
Expected value of
s
Constant
μ
Ξ
Expected value of
Ξ
Constant
μ
r
Expected value of
r
Constant
σ
Ξ
Uncertainty of
Ξ
Constant
S
Surface bg covariance matrix
Constant
R
Resolution covariance matrix
Constant
f
χ
Signal PDF
Function
f
bb
PDF of
bb
Function
ρ
sb
Event density function of
sb
Function
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-10
q
ð
ν
χ
Þ¼
8
>
>
<
>
>
:
−
2
ln
L
ð
ν
χ
;
ˆ
ˆ
θ
;m
χ
;
σ
n
Þ
L
ð
ˆ
ν
χ
;
ˆ
θ
;m
χ
;
σ
n
Þ
!
ν
χ
>
ˆ
ν
χ
0
ν
χ
<
ˆ
ν
χ
;
ð
17
Þ
where
ν
χ
is the number of signal events, and
θ
is the vector
of nuisance parameters. The number of signal events given
by the global likelihood maximum (
ˆ
ν
χ
) corresponds to the
best fit values of the nuisance parameters (
ˆ
θ
). The best fit
values of the nuisance parameters when the signal model is
fixed is indicated by
ˆ
ˆ
θ
. In order to calculate an upper limit,
q
was set to zero when the number of signal events being
tested was lower than the best fit number of signal events.
Equation
(17)
was evaluated over a grid of mass and
cross section points. At each point being tested, the signal
shape was held constant, and the number of signal events
(
ν
χ
) was increased until the value of
q
exceeded 1.64 (as
shown by Fig.
6
), which corresponds to the 90% confidence
level based on Wilks
’
s theorem. If the upper limit was
greater than the predicted number of signal events, then the
signal hypothesis is consistent with the data and that
combination of DM mass and cross section could not be
excluded.
VII. RESULTS
No significant excess of events was observed above the
expected background rate based on a background only fit to
the data, so the procedure outlined above was used to
calculate a set of bands of DM mass and DM-nucleon cross
section that are excluded at the 90% confidence level. This
process was repeated for both inelastic scattering channels
considered in this analysis. The region outside the bands
could not be excluded because more extreme cross sections
result in too few signal events; small cross sections produce
a small signal rate, and large cross sections cause stronger
attenuation of the signal.
Figure
7
shows the observed limit and projected sensi-
tivity of the Migdal search estimated using pseudo datasets.
Figure
8
shows the exclusion regions and the current state
of low-mass DM direct detection searches, including
bremsstrahlung and Migdal channel searches. New param-
eter space that was not previously tested by other DM
searches is excluded via the Migdal channel for DM masses
between 0.032 and
0
.
1
GeV
=c
2
. Although the bremsstrah-
lung result presented in this paper does not exclude any new
parameter space, it is the most sensitive search using this
channel for DM masses between 0.22 and
0
.
4
GeV
=c
2
.To
test the sensitivity of the limit to yield modelling, we cross-
checked our dependence by removing the Gaussian
FIG. 6. Example of a profile likelihood ratio (PLR) scan of
q
for
a given mass and cross section point, with the number of expected
signal events given by N_expected. The upper limit on the
number of signal events is determined by the position where the
ratio crosses 1.64.
FIG. 7. Mean (dotted line) expected sensitivity with
1
σ
(green)
and
2
σ
(yellow) bands and the final limit (solid line).
FIG. 8. Current state of low-mass DM direct detection searches,
with results from SuperCDMS-CPD
[71]
, EDELWEISS
[72]
, and
cosmic ray bounds from Collar
[73]
for traditional elastic
interaction searches. The other curves are published limits from
inelastic channel searches: LUX
[74]
, EDELWEISS
[72,75]
,
XENON1T
[76]
, CDEX-10
[77]
, DarkSide-50
[78]
, and this
result. CRESST results currently in review
[79]
. The green and
yellow bands surrounding the Migdal result indicate the
1
σ
and
2
σ
projected sensitivity ranges, as shown in Fig.
7
.
SEARCH FOR LOW-MASS DARK MATTER VIA
...
PHYS. REV. D
107,
112013 (2023)
112013-11
uncertainty constraint term in the likelihood and allowed
our model to float unconstrained in the fit, corresponding to
no prior knowledge of yield. The impact on the final limit
was negligible.
The low CDMSlite energy threshold allows the brems-
strahlung channel limit to extend below the
0
.
4
GeV
=c
2
mass reach of the LUX and XENON1T bremsstrahlung
analyses. However, both experiments have larger expo-
sures, and so place lower cross section limits at higher
masses
[74,76]
.
Comparing the various Migdal channel results, this result
also extends to slightly lower masses than the LUX and
XENON1T limits because of the lower threshold
[74,76]
.
CDEX has an energy threshold that is to the three times
higher than that of CDMSlite; thus, the integrated rate of
the Migdal signal is smaller and results in a slightly less
sensitive limit
[80]
. The EDELWEISS Migdal limit is not
as sensitive in cross section due to the low exposure from
operating a 33.4 g detector for 24 hours. However, the
shallow depth at which the EDELWEISS dataset was
acquired allows the exclusion region to extend to higher
cross sections than this analysis because these particles
would lose their energy from scattering in the Earth before
reaching the SuperCDMS experiment
[72]
. Similarly,
SuperCDMS-CPD data were collected at a surface facility
using a low threshold detector searching for direct NR
events
[71]
, making these data a prime candidate for
repeating a bremsstrahlung or Migdal analysis in the future.
In summary, this analysis of CDMSlite data accounts for
the shielding of strongly interacting DM particles by the
Earth and atmosphere. This was implemented by calculat-
ing a damping parameter for the DM velocity distribution
and includes angular dependence of the incident DM.
Using a profile likelihood framework, the bremsstrahlung
channel was not found to probe any new parameter space,
but the Migdal effect channel excludes new parameter
space between 32 and
39
MeV
=c
2
.
ACKNOWLEDGMENTS
The authors would like to thank Matthew Dolan, Timon
Emken, Masahiro Ibe, Felix Kahlhoefer, Chris McCabe,
Wakutaka Nakano, Jayden Newstead, Yutaro Shoji, and
Kazumine Suzuki for their useful discussions. The
SuperCDMS collaboration gratefully acknowledges tech-
nical assistance from the staff of the Soudan Underground
Laboratory and the Minnesota Department of Natural
Resources. The CDMSlite and iZIP detectors were fab-
ricated 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. Fermilab 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 DE-AC02-37407CH11359, DE-
AC02-76SF00515, and DE-AC05-76RL01830, respec-
tively. E. L. A. acknowledges the support from Spanish
grant CA3/RSUE/2021-00827, funded by Ministerio de
Universidades, Plan de Recuperacion, Transformacion y
Resiliencia, and Universidad Autonoma de Madrid.
APPENDIX:
Parameters for the earth and atmosphere are listed in
Table
IV
. The density of the Earth was taken from Ref.
[81]
.
TABLE IV. Model parameters used in the calculation of the velocity distribution damping.
[82
–
86]
.
Atmosphere
Soudan
Earth Crust Earth Mantle Earth Outer Core Earth Inner Core
Outer Radius [km]
Table
II
0.7135 (depth)
6
.
371
×
10
3
6
.
331
×
10
3
3
.
46
×
10
3
1
.
22
×
10
3
Density [g
=
cm
3
]
Table
II
see aux data
3.1
5.514
11
12.6
Mass Fraction H
Mass Fraction C
0.0002
Mass Fraction O
0.231
0.476
0.503
0.0273
0.0273
Mass Fraction Na
0.029
Mass Fraction Mg
0.015
0.256
Mass Fraction Al
0.083
Mass Fraction Si
0.283
0.241
0.0509
0.0509
Mass Fraction P
Mass Fraction K
0.027
(Table continued)
M. F. ALBAKRY
et al.
PHYS. REV. D
107,
112013 (2023)
112013-12