[table]capposition=above

# A Search for Low-mass Dark Matter via Bremsstrahlung Radiation and the
Migdal Effect in SuperCDMS

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. Alonso  Instituto de Física Teórica UAM/CSIC, Universidad
Autónoma de Madrid, 28049 Madrid, Spain  Instituto de Física Teórica UAM-CSIC,
Campus de Cantoblanco, 28049 Madrid, Spain  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, HBNI, Jatni -
752050, India  C. Bathurst  Department of Physics, University of Florida,
Gainesville, FL 32611, USA  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  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  Corresponding
author<EMAIL_ADDRESS>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  Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de
Madrid, 28049 Madrid, Spain  Instituto de Física Teórica UAM-CSIC, Campus de
Cantoblanco, 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, HBNI,
Jatni - 752050, 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  S. Das  School
of Physical Sciences, National Institute of Science Education and Research,
HBNI, Jatni - 752050, India  F. De Brienne  Département de Physique,
Université de Montréal, Montréal, Québec H3C 3J7, Canada  M. Rios  Instituto
de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, 28049 Madrid,
Spain  Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049
Madrid, Spain  S. Dharani  Institute for Astroparticle Physics (IAP),
Karlsruhe Institute of Technology (KIT), 76344 Eggenstein-Leopoldshafen,
Germany  Institut für Experimentalphysik, Universität Hamburg, 22761 Hamburg,
Germany  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  M. Elwan  Department of Physics, University of
Florida, Gainesville, FL 32611, USA  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  College of Natural
and Health Sciences, Zayed University, Dubai, 19282, United Arab Emirates
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  N. Hassan  Département de Physique, Université de Montréal,
Montréal, Québec H3C 3J7, Canada  B.A. Hines  Department of Physics,
University of Colorado Denver, Denver, CO 80217, 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  Department of Physics, University of Toronto, Toronto,
ON M5S 1A7, Canada  D. Jardin  Corresponding author:
<EMAIL_ADDRESS>Department of Physics, Southern Methodist
University, Dallas, TX 75275, USA  V.K.S. Kashyap  School of Physical
Sciences, National Institute of Science Education and Research, HBNI, Jatni -
752050, 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  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  M. Litke  Department of
Physics, Southern Methodist University, Dallas, TX 75275, USA  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  E. Lopez
Asamar  Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid,
28049 Madrid, Spain  Instituto de Física Teórica UAM-CSIC, Campus de
Cantoblanco, 28049 Madrid, Spain  P. Lukens  Fermi National Accelerator
Laboratory, Batavia, IL 60510, USA  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  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
Eggenstein-Leopoldshafen, 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, HBNI, Jatni - 752050, India  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  S.
Pandey  School of Physics & Astronomy, University of Minnesota, Minneapolis,
MN 55455, 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  L. Perna  Department of Physics, University of
Toronto, Toronto, ON M5S 1A7, 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  D. Sadek  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  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  S.S. Poudel  Pacific Northwest National Laboratory,
Richland, WA 99352, USA  D.J. Sincavage  School of Physics & Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA  P. Sinervo  Department of
Physics, University of Toronto, Toronto, ON M5S 1A7, Canada  J. Street
Department of Physics, South Dakota School of Mines and Technology, Rapid
City, SD 57701, USA  H. Sun  Department of Physics, University of Florida,
Gainesville, FL 32611, USA  G.D. Terry  Department of Physics, University of
South Dakota, Vermillion, SD 57069, 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 Eggenstein-
Leopoldshafen, 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 Eggenstein-Leopoldshafen,
Germany  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  C.-P. Wu  Département de Physique, Université de
Montréal, Montréal, Québec H3C 3J7, Canada  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 Eggenstein-Leopoldshafen, 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  A. Zuniga  Department of Physics, University of
Toronto, Toronto, ON M5S 1A7, Canada

###### Abstract

We present a new analysis of previously published of 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~{}\textrm{MeV}/c^{2}$ at $2.7\times 10^{-30}~{}\textrm{cm}^{2}$ via the
bremsstrahlung channel. The Migdal channel search provides overall
considerably more stringent limits and excludes DM masses down to
$30~{}\textrm{MeV}/c^{2}$ at $5.0\times 10^{-30}~{}\textrm{cm}^{2}$.

SuperCDMS Collaboration

##  I Introduction

An abundance of evidence suggests that most of the Universe is composed of
non-luminous matter [1, 2, 3]. This “dark matter” (DM) may consist of an
undiscovered elementary particle or a set of particles [4]. 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~{}\textrm{GeV}/c^{2}$ through
inelastic detection channels. In canonical direct DM searches, the interaction
between a DM particle and a nucleus is assumed to be an elastic two-body
interaction. For DM particles with masses, $m_{\chi}$, much smaller than that
of the target nucleus $m_{N}$, the recoil energy from an elastic collision is
suppressed by the kinematic term $m_{\chi}^{2}/m_{N}$, resulting in rapidly
diminishing sensitivity when considering lower mass DM candidates. 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 scattering process when $m_{\chi}<<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
[7, 6]. Therefore, these results are directly comparable to existing limits
for DM-nucleon interactions.

In this paper, we present a re-analysis 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 Section VI
describes the limit-setting method. The final results are presented in Section
VII.

##  II SuperCDMS

The SuperCDMS experiment was operated $\sim$ 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, 9, 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 were
instrumented with interleaved Z-sensitive Ionization and Phonon (iZIP) sensors
[11]. The detectors were biased face-to-face with $\pm~{}2~{}\textrm{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~{}\textrm{keV}$ [12].

Two detectors were, for some periods, operated with a high-voltage bias near
$\sim$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~{}\textrm{eV}_{\textrm{ee}}$ (ER equivalent energy) [9, 10],
but sacrifices all discrimination between NR and ER events.

###  II.1 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~{}\textrm{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.

###  II.2 Event Selection

Since the energy region of interest for this analysis largely overlaps with
Ref. [10], the same data quality selection criteria were used to remove
problematic events such as those arising from low-frequency mechanical noise,
electronic glitches, and poorly reconstructed 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 amplification is not uniform throughout the crystal. To
minimize the number of events with reduced amplification, the same fiducial
volume selection as defined in Ref. [10] was adopted, which rejects events
with low NTL amplification. After applying the selection criteria, the
remaining exposure is $36.9~{}\textrm{kg}\cdot\textrm{d}$ and the analysis
threshold is $70~{}\textrm{eV}_{\textrm{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 characteristic velocity of
$220~{}\textrm{km/s}$ [15]. Particles traveling at velocities greater than the
escape velocity of the galaxy are not gravitationally bound, so the
distribution is truncated at $544~{}\textrm{km/s}$ and re-normalized [16]. The
local DM density is assumed to be
$0.3~{}\textrm{GeV}/(c^{2}\;\textrm{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
$\sim$$10^{-30}~{}\textrm{cm${}^{2}$}$, the Earth and atmosphere can no longer
be neglected. As DM 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 [18].

This analysis accounts for the attenuation effect by damping the DM velocity
distribution, described as “Method B” in Ref. [19]. 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 Section III.1.

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. [19, 20]. Since we concern
ourselves with light dark matter in this paper, we assume a nuclear form
factor of unity. A velocity damping parameter, $\kappa$, is defined as:

$\kappa\equiv\frac{\rho\sigma_{n}}{m_{\chi}\mu_{n}^{2}}\left(\sum_{i}^{\mathrm{elements}}\frac{F_{i}\mu_{i}^{4}A_{i}^{2}}{m_{i}^{2}}\right)d\,,$
(1)

with the variables defined in Table 1 [19]. The calculation of $\kappa$ is
location specific since it depends on the path length, $d$, and the type of
material the particle travels through before reaching the detector.

Table 1: Definition of variables in Equation 1 to calculate the velocity damping parameter, $\kappa$. Variable | Definition
---|---
$\rho$ | material density
$\sigma_{n}$ | DM-nucleon scattering cross section
$\mu_{n}$ | reduced mass (DM particle, nucleon)
$i$ | element index
$F_{i}$ | mass fraction of the element $i$
$\mu_{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

An incoming DM particle with initial velocity $v_{i}$ reaches the detector
with velocity $v_{f}$,

$v_{f}=v_{i}e^{-\kappa}\,,$ (2)

thus the DM velocity distribution at the detector ($f_{d}$) is given by:

$f_{d}\left(v_{f}\right)=e^{2\kappa}\times f\left(v_{f}e^{\kappa}\right)\,,$
(3)

where the incoming SHM velocity distribution ($f)$ has been transformed by the
effects of shielding [19]. The exponential term in front of the distribution
accounts for the normalization of the increased flux caused by the attenuated
velocity distribution.

Implicit in the definition of $\kappa$ 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
overburden, 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 dependence can be found in
Section III.2. We assume that the particles travel in straight line
trajectories even though scattering will affect their trajectories. According
to Ref. [19], this approach still underestimates the number of dark matter
particles with sufficient energy to interact with the detector.

###  III.1 Earth Shielding

Due to the exponential nature of Equation 3, $\kappa$ can be calculated for
each layer of shielding independently using Equation 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 continuously as a function of
altitude. However, for the purpose of simplifying the model, a seven layer
atmosphere model based on the 1976 US Standard Atmosphere was used [21]. 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 2.

Layer | Height range [km] | Maximum density [ g/m3 ]
---|---|---
1 | 0-11 | 1225
2 | 11-20 | 363.91
3 | 20-32 | 88.03
4 | 32-47 | 13.22
5 | 47-51 | 1.43
6 | 51-71 | 0.86
7 | 71-100 | 0.064
Table 2: Densities and maximum heights of the layers in the atmosphere model
[21].

The SuperCDMS experiment at Soudan was located under $714~{}\textrm{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 [22]. 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.111 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.  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
$\kappa$ for each mass and cross section considered.

The dominant component of the shielding is the Earth below the experiment. The
conventional model of the Earth is 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.

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~{}\textrm{cm}$ of polyethylene (C2H4) with a density of
$0.94~{}\textrm{g/cm${}^{3}$}$, followed by $22.5~{}\textrm{cm}$ of Pb and
$3~{}\textrm{cm}$ of Cu [23]. The contribution to the velocity damping
parameter from the shielding around the experiment is 0.2% for a 1
$\textrm{GeV}/c^{2}$ DM particle with a nucleon scattering cross section of
$10^{-30}~{}\textrm{cm}^{2}$ and traveling straight downward, where the
Earth’s shielding contribution is the weakest.

###  III.2 Angular Dependence

The depth of the experiment in the Earth’s crust produces an asymmetric
angular distribution of shielding. 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.

This analysis follows the angular convention in Ref. [24], which defines
$\theta$ as the incident angle between the incoming DM particle’s velocity and
zenith at the experiment. By this definition, $\theta=0$ corresponds to
particles originating from directly below the experiment. The total path
length is calculated as:

$\ell=r_{d}\cos\theta+\sqrt{r_{E}^{2}-\left(r_{d}\sin\theta\right)^{2}}\,,$
(4)

where $r_{d}$ is the distance from the center of the Earth to the SuperCDMS
experiment. The average value of $6471~{}\textrm{km}$ was used for the Earth’s
radius, $r_{E}$. The total path length can be generalized to determine $d$,
the path length through individual layers, used in Equation 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 $\kappa$
that affects the velocity distribution at the detector. The signal model was
integrated over the incoming angle, $\theta$, at twenty discrete points
sampled uniformly in $\cos\theta$ as indicated in Figure 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.

Figure 1: Value of $\kappa$ as a function of incident angle, $\theta$, for a
DM particle with a mass of $1~{}\textrm{GeV}/c^{2}$ and a scattering cross
section of $10^{-30}~{}\textrm{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 $\pi/2$ indicate the location of the
transition between the core and the mantel while the more detailed overburden
model leads to a smooth curve above $\pi/2$.

##  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 radiation, 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.

###  IV.1 Bremsstrahlung Radiation

The differential scattering rate for emitting a photon of energy $E_{\gamma}$
through the bremsstrahlung process has been derived by Kouvaris and Pradler in
Ref. [6]:

$\frac{d\sigma}{dE_{\gamma}}=\frac{4\alpha|f(E_{\gamma})|^{2}}{3\pi
E_{\gamma}}\frac{\mu_{N}^{2}v^{2}\sigma_{N}}{m_{N}^{2}}\sqrt{1-\frac{2E_{\gamma}}{\mu_{N}v^{2}}}\left(1-\frac{E_{\gamma}}{\mu_{N}v^{2}}\right)\,,$
(5)

where $\alpha$ is the fine structure constant, $f$ is the atomic scattering
function discussed in Section IV.1.1, $\mu_{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 $\sigma_{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,

$\frac{dR}{dE_{\gamma}}=N_{T}\frac{\rho_{\chi}}{m_{\chi}}\int_{|\vec{v}|\geq
v_{\rm min}}d^{3}\vec{v}\,vf_{d}(\vec{v})\frac{d\sigma}{dE_{\gamma}}\,,$ (6)

where $N_{T}$ is Avogadro’s number divided by the atomic mass. The minimum
velocity required to induce a recoil of $E_{\gamma}$, $v_{\rm min}$, is given
by

$v_{\rm
min}=\sqrt{\frac{m_{N}E_{R}}{2\mu^{2}_{N}}}+\frac{E_{\gamma}}{\sqrt{2m_{N}E_{R}}},$
(7)

where $E_{R}$ is the nuclear recoil energy.

####  IV.1.1 Photoelectric Absorption Cross Section Uncertainty

The atomic scattering function, $f$, in Equation 5 is the sum of a real and
complex portion,

$|f|^{2}=|f_{1}+if_{2}|^{2}=f_{1}^{2}+f_{2}^{2}\,.$ (8)

The components, $f_{1}$ and $f_{2}$, are defined in Ref. [25] as:

$f_{1}=Z^{*}+\frac{1}{\pi
r_{e}hc}\int_{0}^{\infty}\frac{\epsilon^{2}\sigma_{a}(\epsilon)}{E_{\gamma}^{2}-\epsilon^{2}}d\epsilon\,;\quad
f_{2}=\frac{\sigma_{a}}{2r_{e}\lambda}\,,$ (9)

for the nuclear contribution, where $r_{e}$ is the classical radius of the
electron, $\sigma_{a}$ is the photoelectric absorption cross section,
$\lambda$ is the wavelength of the emitted photon and $\epsilon$ is a variable
of integration. $Z^{*}$ is defined as $Z-(Z/82.5)^{2.37}$, where $Z$ is the
proton number.

Measurements of $\sigma_{a}$ at low energies have a range of values which lead
to significant systematic uncertainties of its value [26]. Both $f_{1}$ and
$f_{2}$ depend on $\sigma_{a}$, so this uncertainty enters into the
calculation of the expected event rate. Based on existing literature [27, 28,
29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41], a nominal value and
uncertainty band for $\sigma_{a}$ were derived as a function of energy. Using
the resulting values of $\sigma_{a}$, $|f|^{2}$ was calculated and compared in
Figure 2, along with the commonly referenced Henke dataset [27]. 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 calculations of $|f|^{2}$,
that conservatively predicts a lower expected signal rate and thus results in
a weaker limit.

Figure 2: The atomic scattering function, $|f|^{2}$, calculated as a function
of energy using different values of the photoelectric absorption cross
section, $\sigma_{a}$, in germanium. The nominal value of $\sigma_{a}$
obtained from literature was used to calculate the solid black curve. The
range of uncertainty on $|f|^{2}$ is indicated by the shaded region.

###  IV.2 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 simultaneously with respect to the nucleus [42]. The displacement of the
nucleus due to a DM scattering event dramatically changes the wavefunction of
the electrons in the surrounding 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 [42]. The formulation of this
effect as applied to DM direct detection has been calculated by Ibe, Nakano,
Shoji, and Suzuki [7]. An alternative approach utilizing the photoelectric
cross section has been developed by Liu, Wu, Chi, and Chen [43]. The Ibe et
al. formulation was then applied to data by Dolan, Kahlhoefer, and McCabe
[44]. 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:

$\frac{\mathrm{d}^{3}R_{\text{ion}}}{\mathrm{d}E_{\mathrm{R}}\,\mathrm{d}E_{e}\,\mathrm{d}v}=\frac{\mathrm{d}^{2}R_{\text{nr}}}{\mathrm{d}E_{\mathrm{R}}\,\mathrm{d}v}\times\sum_{nl}\frac{1}{2\pi}\frac{\mathrm{d}p^{c}_{q_{e}}(nl\to
E_{e})}{\mathrm{d}E_{e}}\,,$ (10)

where $p_{q_{e}}^{c}(nl\to E_{e})$ is the transition probability for an
electron, with momentum $q_{e}=m_{e}v_{\mathrm{nucleus}}$ with respect to a
stationary nucleus, to be ejected from quantum state $nl$ with energy $E_{e}$,
and $\mathrm{d}^{2}R_{\text{nr}}/(\mathrm{d}E_{\mathrm{R}}\,\mathrm{d}v)$
describes the incoming DM flux and scattering rate,

$\frac{\mathrm{d}^{2}R_{\text{nr}}}{\mathrm{d}E_{\mathrm{R}}\,\mathrm{d}v}=\frac{\rho_{\chi}\,\sigma_{N}}{2\,\mu_{N}^{2}\,m_{\chi}}\frac{f_{d}(v)}{v}\,.$
(11)

The factor of 1/2${\pi}$ in Equation 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_{\mathrm{ref}}=10^{-3}c$. Conversion between the reference value and an
arbitrary electron momentum is given by:

$p_{q_{e}}^{c}\left(nl\rightarrow
E_{e}\right)=\left(\frac{q_{e}}{v_{\mathrm{ref}}m_{e}}\right)^{2}p_{v_{\mathrm{ref}}}^{c}\left(nl\rightarrow
E_{e}\right)\,.$ (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 [45, 46, 47]. 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~{}\textrm{GeV}/c^{2}$ mass and $10^{-35}~{}\textrm{cm}^{2}$ 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.

Figure 3: Comparison of the expected DM signal rates in germanium. For the DM
mass chosen, $0.5~{}\textrm{GeV}/c^{2}$, the elastic NR signal is below the
analysis threshold of $70~{}\textrm{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~{}\textrm{keV}$ is where
the $n=2$ electron shell starts contributing to the signal rate.

##  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 252Cf calibration source, cosmogenic
activation of the crystal, radiogenic Compton scattering, and 210Pb 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 $\ll 1$ event for the exposure, and therefore
negligible.

###  V.1 Bulk Background Models

The dominant background originates from activation of the germanium crystal
via the neutron calibration source. When stable 70Ge in the crystal captures a
neutron, it becomes unstable 71Ge that decays via electron capture to 71Ga.
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~{}\textrm{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~{}\textrm{keV}$, and the M-shell captures an electron $1.78\%$
of the time and emits $160~{}\textrm{eV}$ [48]. This background was modeled by
a Gaussian distribution centered on each electron shell peak energy. The
amplitudes 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 71Ge background. Similarly, 68Ge decays to 68Ga which can
decay through electron capture or through beta decay. The 68Ga beta endpoints
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 (3H), which
provides a persistent source of betas due to its half-life of
$12.32~{}\textrm{years}$. The tritium background was modeled by a standard
beta emission spectrum with an endpoint energy of $18.6~{}\textrm{keV}$ and an
unconstrained normalization. The resulting normalization from the fit was 50
$\pm$ 20 events, which is consistent with a previous dedicated analysis [49,
50].

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 68Ga, 65Zn, and 55Fe [50]. Other isotopes
considered are 57Co, 54Mn, and 49Vn, 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 71Ge was treated above.

All of the radioisotopes, created cosmogenically or by source activation,
described in this subsection are distributed nearly uniformly throughout the
detector volume222Studies based on simulation have shown slightly more bias
towards the surface of the crystals due to self-shielding effects.. 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 [51, 10]. 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].

###  V.2 Surface Background Models

Radon daughters plating out on the detectors or surrounding copper housing
were treated differently than the bulk contamination described in Section V.1.
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 [52, 53, 54] simulation of 210Pb
surface contamination on a tower of six germanium detectors was performed. The
simulation allowed for the subsequent alpha decays to implant the long-lived
210Pb and mimic the physical radiocontamination [51]. 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 independent measurement of the
rate of alpha events in each of the detectors [51]. For the energy range used
in this analysis, the normalization to the number of expected events from each
contribution were determined to be $N_{\mathrm{TL}}=158.4\pm 16.6$,
$N_{\mathrm{SH}}=22.3\pm 1.4$, and $N_{\mathrm{GC}}=23.5\pm 5.9$ events, which
contribute to the total number of background events in the likelihood. The
housing copper is not as radiopure 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 uncertainties are taken into account using morphing parameters as
described in Ref. [10].

###  V.3 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~{}\textrm{eV}_{\textrm{ee}}$ and 100% above
$90~{}\textrm{eV}_{\textrm{ee}}$. The data quality selection criteria
efficiency approaches 100% around $150~{}\textrm{eV}_{\mathrm{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~{}\textrm{eV}_{\mathrm{ee}}$ and has almost zero efficiency at
$70~{}\textrm{eV}_{\mathrm{ee}}$.

Figure 4 shows the efficiency curve for each period used in this analysis. The
notable fluctuations in the efficiency curve below
$2~{}\textrm{keV}_{\mathrm{ee}}$ arise from the fiducial volume selection,
which was calculated using simulated data with limited statistics. We take the
efficiency curve with its statistical fluctuations. The efficiency models have
been extended to $25~{}\textrm{keV}_{\textrm{ee}}$ following the procedure in
Ref. [50].

Figure 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\sigma$ uncertainty band.

The uncertainty on the efficiency curve was incorporated 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
accomplished by constructing 1-sigma bands around the median parameterized as
a Gaussian distribution following the same prescription as the surface
background models in Ref. [10].

###  V.4 Resolution Model

The background model and signal models were convolved 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:

$\sigma^{2}(E)=\sigma_{E}^{2}+BE+(AE)^{2}\,,$ (13)

where $\sigma$ is the resolution, $E$ is the total phonon energy, $\sigma_{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 Section II.D of Ref. [10]. At 100
eVee, the resolution of period 1 and 2 is 14 eVee and 16 eVee respectively.
Near the upper end of the analysis range at 25 keVee, the resolution is 193
eVee and 198 eVee for periods 1 and 2, respectively.

###  V.5 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_{\mathrm{NR}})$, is needed. We adopt the
Lindhard model (Equation 14),

$Y_{\text{Lindhard}}(E_{\mathrm{NR}})=k\cdot g(\epsilon)/(1+k\cdot
g(\epsilon)),$ (14)

where $g(\epsilon)=3\epsilon^{0.15}+0.7\epsilon^{0.6}+\epsilon$,
$\epsilon=11.5E_{\mathrm{NR}}({\mathrm{keV}})Z^{-7.3}$, $Z$ is the atomic
number, and $k$ is the free electron energy loss [9, 55].

There is evidence of deviations from the Lindhard model at low energies [56,
57]. Therefore, the ionization production was cut off at 22.7 eV, according to
an independent measurement of the defect energy creation threshold [58]. This
is implemented with a hyperbolic tangent function,

$\displaystyle Y(E_{\mathrm{NR}})=$
$\displaystyle\frac{1}{2}Y_{\text{Lindhard}}(E_{\mathrm{NR}})$ (15)
$\displaystyle\times\left(\tanh\left(\frac{E_{\mathrm{NR}}-22.7\;\mathrm{eV}}{22.32\;\mathrm{eV}}\right)+1\right),$

that 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 $\sigma=0.05$ [10]. The signal
model is calculated for the nominal, $\pm 1\sigma$, and $\pm 3\sigma$ values
of k and intermediate values of $\sigma$ are interpolated.

###  V.6 Analysis Energy Range

The normalization of the models was determined by the fit to data, as
described in Section VI. An example of fitting the data is shown in Figure 5,
where the background models and Migdal signal model, for a WIMP with a mass of
0.5 GeV/c2 and a cross section of 3$\times 10^{-37}\mathrm{cm}^{2}$, have been
scaled by the efficiency and convolved with the resolution model.

Figure 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/c2 and a cross section
of $3\times 10^{-37}$ cm2 (black dashed curve). The data (blue histogram) have
been logarithmically binned and overlaid with the background models (colored
solid curves). The thick black line is the sum of all the models.
Normalization of the surface background model components (TL, SG and GC) are
described in Section V.2. The plot on the bottom shows the residual between
data and the model with the 1 $\sigma$ statistical uncertainty indicated by
the shaded region.

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~{}\textrm{keV}_{\textrm{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 [59, 60] technique because it incorporates knowledge of the
background. It also provides a rigorous and convenient way to account for
systematic uncertainties in signal and background model parameters.

###  VI.1 Likelihood Function

Table 3 contains a list of the variables used in the likelihood function. The
likelihood function is composed of three types of terms. The first type is the
overall normalization term, ($\mathcal{L}_{\textrm{Poiss}}$), which allows the
total number of fitted events to fluctuate around the number of events in the
dataset ($N$). The second type ($\mathcal{L}_{\chi,bb,sb}^{P1}$,
$\mathcal{L}_{\chi,bb,sb}^{P2}$) 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 ($\mathcal{L}_{\textrm{surf}}$), the morphed efficiency
($\mathcal{L}_{\textrm{eff}}^{P1}$, $\mathcal{L}_{\textrm{eff}}^{P2}$), the
resolution model ($\mathcal{L}_{\textrm{res}}^{P1}$,
$\mathcal{L}_{\textrm{res}}^{P2}$), and the yield model
($\mathcal{L}_{\textrm{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 covariance matrix is calculated from
the normalization uncertainty of the individual components and the
correlations between them [61].

Table 3: Definition of variables used in the likelihood function (Equation 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
---|---|---
$P1$ | data period 1 | Identifier
$P2$ | 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
$\nu_{\chi}$ | number of signal events | Free
$\nu_{bb}$ | number of events in $bb$ | Free
$\nu_{sb}$ | number of events in $sb$ | Constrained
$s$ | surface bg morphing parameter | Constrained
$\Xi$ | efficiency morphing parameter | Constrained
$r$ | resolution nuisance parameter | Constrained
$\theta_{Y}$ | yield nuisance parameter | Constrained
$N$ | number of events in data | Constant
$E_{n}$ | energy of event $n$ | Constant
$\mu_{s}$ | expected value of $s$ | Constant
$\mu_{\Xi}$ | expected value of $\Xi$ | Constant
$\mu_{r}$ | expected value of $r$ | Constant
$\sigma_{\Xi}$ | uncertainty of $\Xi$ | Constant
$\boldsymbol{S}$ | surface bg covariance matrix | Constant
$\boldsymbol{R}$ | resolution covariance matrix | Constant
$f_{\chi}$ | signal PDF | Function
$f_{bb}$ | PDF of $bb$ | Function
$\rho_{sb}$ | event density function of $sb$ | Function

This analysis used the MINUIT algorithm [62] via the iminuit [63] Python
interface to maximize the log-likelihood function when evaluating the test
statistic, which will be defined in Equation 17.

The full likelihood function used in this analysis is given by Equation 16,
which is a single likelihood encompassing both data periods.

$\begin{gathered}\ln\mathcal{L}=\\\\[10.0pt]
\ln(\mathcal{L}_{\textrm{Poiss}})+\ln(\mathcal{L}_{\chi,bb,sb}^{P1})+\ln(\mathcal{L}_{\chi,bb,sb}^{P2})+\ln(\mathcal{L}_{\textrm{surf}})\\\
+\ln(\mathcal{L}_{\textrm{eff}}^{P1})+\ln(\mathcal{L}_{\textrm{eff}}^{P2})+\ln(\mathcal{L}_{\textrm{res}}^{P1})+\ln(\mathcal{L}_{\textrm{res}}^{P2})+\ln(\mathcal{L}_{\textrm{yield}})=\\\\[10.0pt]
-\left(\nu_{\chi}+\sum\limits_{bb=1}^{6}\nu_{bb}+\sum\limits_{sb=1}^{3}\nu_{sb}\right)
\\\
+\sum_{n=1}^{N^{P1}}\ln\left[\nu_{\chi}f_{\chi}^{P1}(E_{n})+\sum\limits_{bb=1}^{6}\nu_{bb}f_{bb}^{P1}(E_{n})+\sum\limits_{sb=1}^{3}\rho_{sb}^{P1}(E_{n})\right]\\\
+\sum_{n=1}^{N^{P2}}\ln\left[\nu_{\chi}f_{\chi}^{P2}(E_{n})+\sum\limits_{bb=1}^{6}\nu_{bb}f_{bb}^{P2}(E_{n})+\sum\limits_{sb=1}^{3}\rho_{sb}^{P2}(E_{n})\right]\\\
-\frac{1}{2}\sum\limits_{i,j=1}^{3}\left[(s_{i}-\mu_{s_{i}})^{T}(\boldsymbol{S}_{i,j})^{-1}(s_{j}-\mu_{s_{j}})\right]\\\
-\frac{{(\Xi^{P1}-\mu_{\Xi}^{P1})}^{2}}{2{(\sigma_{\Xi}^{P1})}^{2}}-\frac{{(\Xi^{P2}-\mu_{\Xi}^{P2})}^{2}}{2{(\sigma_{\Xi}^{P2})}^{2}}\\\
-\frac{1}{2}\sum\limits_{i,j=1}^{3}\left[(r_{i}^{P1}-\mu_{r_{i}}^{P1})^{T}(\boldsymbol{R}_{i,j}^{P1})^{-1}(r_{j}^{P1}-\mu_{r_{j}}^{P1})\right]\\\
-\frac{1}{2}\sum\limits_{i,j=1}^{3}\left[(r_{i}^{P2}-\mu_{r_{i}}^{P2})^{T}(\boldsymbol{R}_{i,j}^{P2})^{-1}(r_{j}^{P2}-\mu_{r_{j}}^{P2})\right]\\\
-\frac{1}{2}\theta_{Y}^{2}\,.\end{gathered}$ (16)

###  VI.2 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 due to the shielding discussed in Section
III. According to Equation 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 Equation 1 and the shape of the expected signal to
depend on the DM mass and cross 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:

$q(\nu_{\chi})=\left\\{\begin{array}[]{cr}-2\ln\left(\frac{\mathcal{L}(\nu_{\chi},\hat{\hat{\theta}},m_{\chi},\sigma_{n})}{\mathcal{L}(\hat{\nu}_{\chi},\hat{\theta},m_{\chi},\sigma_{n})}\right)&\nu_{\chi}>\hat{\nu}_{\chi}\\\
0&\nu_{\chi}<\hat{\nu}_{\chi}\end{array}\right.\,,$ (17)

where $\nu_{\chi}$ is the number of signal events, and $\theta$ is the vector
of nuisance parameters. The number of signal events given by the global
likelihood maximum ($\hat{\nu}_{\chi}$) corresponds to the best fit values of
the nuisance parameters ($\hat{\theta}$). The best fit values of the nuisance
parameters when the signal model is fixed is indicated by
$\hat{\hat{\theta}}$. 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 ($\nu_{\chi}$) was increased until the value of $q$ exceeded
1.64 (as shown by Figure 6), which corresponds to the 90% confidence level
based on Wilk’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.

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

##  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: Mean (dotted line) expected sensitivity with 1$\sigma$ (green) and
2$\sigma$ (yellow) bands and the final limit (solid line).

Figure 7 shows the observed limit and projected sensitivity 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 parameter 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 $\textrm{GeV}/c^{2}$. Although the
bremsstrahlung 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 $\textrm{GeV}/c^{2}$. To test the sensitivity of
the limit to yield modelling, we cross-checked our dependence by removing the
Gaussian 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.

Figure 8: Current state of low-mass DM direct detection searches, with results
from SuperCDMS-CPD [64], DAMIC [65], cosmic ray bounds from Collar [66], and
Darkside-50 [67] for traditional elastic interaction searches. The other
curves are published limits from inelastic channel searches: LUX [68],
EDELWEISS [69, 70], XENON1T [71], CDEX [72], CDEX-10 [73], and this result.
The green and yellow bands surrounding the Migdal result indicate the
1$\sigma$ and 2$\sigma$ projected sensitivity ranges, as shown in Figure 7.

The low CDMSlite energy threshold allows the bremsstrahlung channel limit to
extend below the $0.4~{}\textrm{GeV}/c^{2}$ mass reach of the LUX and XENON1T
bremsstrahlung analyses. However, both experiments have larger exposures, and
so place lower cross section limits at higher masses [68, 71].

Comparing the various Migdal channel results, this result also extends to
lower masses than the LUX and XENON1T limits because of the lower threshold,
and is competitive at intermediate masses despite a smaller exposure [68, 71].
CDEX has an energy threshold that is 2–3 times higher than that of CDMSlite,
thus the integrated rate of the Migdal signal is smaller and results in a
slightly less sensitive limit [72]. The EDELWEISS Migdal limit is not as
sensitive in cross section due to the low exposure from operating a
$33.4~{}\textrm{g}$ detector for $24~{}\textrm{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 [69]. Similarly, SuperCDMS-CPD data were collected at a
surface facility using a low threshold detector searching for direct NR events
[64], 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 calculating 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 0.032 and 0.1 $\textrm{GeV}/c^{2}$.

##  VIII Acknowledgements

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

| Atmosphere | Soudan | Earth Crust | Earth Mantle | Earth Outer Core | Earth Inner Core
---|---|---|---|---|---|---
Outer Radius [km] | Table 2 | 0.7135 (depth) | 6.371 x $10^{3}$ | 6.331 x $10^{3}$ | 3.46 x $10^{3}$ | 1.22 x $10^{3}$
Density [g/cm3] | Table 2 | 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 |  |  |
Mass Fraction Ca |  |  | 0.037 |  |  |
Mass Fraction Mn |  |  |  |  |  |
Mass Fraction Fe |  |  | 0.051 |  | 0.8509 | 0.8509
Mass Fraction Ti |  |  |  |  |  |
Mass Fraction S |  |  |  |  | 0.0188 | 0.0188
Mass Fraction N | 0.756 |  |  |  |  |
Mass Fraction Ar | 0.013 |  |  |  |  |
Mass Fraction Ne | 0.00001 |  |  |  |  |
Mass Fraction Ni |  |  |  |  | 0.0520 | 0.0520
Table 4: Model parameters used in the calculation of the velocity distribution
damping. [74, 75, 76, 77, 78]

##  IX Appendix

We describe the rock and chemical composition in SoudanRegion.csv and
RockChem.csv, respectively. The SoudanRegion.csv file contains the area and
fraction of each rock type in eight directions between radii of 100, 500,
1000, 5000, 10000, 20000 and 50000 meters. The elemental mass fractions for
the chemical composition is in RockChem.csv. Parameters for the earth and
atmosphere are listed in Table 4. The density of the Earth was taken from Ref.
[79].

## References

  * Tanabashi  _et al._ [2018] M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Ade  _et al._ [2016] P. A. R. Ade _et al._ (Planck Collaboration), A&A 594, A13 (2016).
  * Clowe  _et al._ [2006] D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, Astrophys. J. 648, L109 (2006), arXiv:astro-ph/0608407 [astro-ph] .
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 [hep-ph] .
  * Battaglieri  _et al._ [2017] M. Battaglieri _et al._ , in  _U.S. Cosmic Visions: New Ideas in Dark Matter College Park, MD, USA, March 23-25, 2017_ (2017) arXiv:1707.04591 [hep-ph] .
  * Kouvaris and Pradler [2017] C. Kouvaris and J. Pradler, Phys. Rev. Lett. 118, 031803 (2017).
  * Ibe  _et al._ [2018] M. Ibe, W. Nakano, Y. Shoji, and K. Suzuki, JHEP 03 (3), 194, arXiv:1707.07258 [hep-ph] .
  * Agnese  _et al._ [2018a] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 120, 061802 (2018a), arXiv:1708.08869 [hep-ex] .
  * Agnese  _et al._ [2018b] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 97, 022002 (2018b), arXiv:1707.01632 [astro-ph.CO] .
  * Agnese  _et al._ [2019a] R. Agnese _et al._ (SuperCDMS), Phys. Rev. D 99, 062001 (2019a), arXiv:1808.09098 [astro-ph.CO] .
  * Agnese  _et al._ [2013] R. Agnese _et al._ (SuperCDMS), Appl. Phys. Lett. 103, 164105 (2013), arXiv:1305.2405 [physics.ins-det] .
  * Agnese  _et al._ [2014] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 112, 10.1103/physrevlett.112.241302 (2014).
  * Neganov and Trofimov [1985] B. S. Neganov and V. N. Trofimov, Otkryt. Izobret. 146, 215 (1985).
  * Luke [1988] P. N. Luke, J. Appl. Phys. 64, 6858 (1988).
  * Kerr and Lynden-Bell [1986] F. J. Kerr and D. Lynden-Bell, Mon. Notices Royal Astron. Soc. 221, 1023 (1986).
  * Smith  _et al._ [2007] M. Smith _et al._ , Mon. Notices Royal Astron. Soc. 379, 755 (2007).
  * Lewin and Smith [1996] J. Lewin and P. Smith, Astropart. Phys. 6, 87 (1996).
  * Emken  _et al._ [2017] T. Emken, C. Kouvaris, and I. M. Shoemaker, Phys. Rev. D 96, 015018 (2017).
  * Emken and Kouvaris [2018a] R. Emken and C. Kouvaris, Phys. Rev. D 97, 115047 (2018a).
  * Starkman  _et al._ [1990] G. D. Starkman, A. Gould, R. Esmailzadeh, and S. Dimopoulos, Phys. Rev. D 41, 3594 (1990).
  * NASA [1976] NASA,  _U.S. Standard Atmosphere, 1976_ , Tech. Rep. NASA-TM-X-74335, NOAA-S/T-76-1562 (National Oceanic and Atmospheric Administration; Department of the Air Force, ; Washington, DC, United States, 1976) 19770009539.
  * Peterson [2020] D. Peterson, personal communication (2020).
  * Akerib  _et al._ [2005] D. S. Akerib _et al._ (CDMS), Phys. Rev. D 72, 052009 (2005), arXiv:astro-ph/0507190 [astro-ph] .
  * Kavanagh [2018] B. J. Kavanagh, Phys. Rev. D 97, 123013 (2018), arXiv:1712.04901 [hep-ph] .
  * Attwood  _et al._ [2009] D. Attwood _et al._ ,  _X-Ray Data Booklet_ , Tech. Rep. LBNLPUB-490 Rev 3 (Lawrence Berkeley National Laboratory, Unversity of California, Berkeley, Ca 94720, 2009).
  * Agnese  _et al._ [2018c] R. Agnese _et al._ (SuperCDMS), Phys. Rev. Lett. 121, 051301 (2018c), [erratum: Phys. Rev. Lett.122,no.6,069901(2019)], arXiv:1804.10697 [hep-ex] .
  * Henke  _et al._ [1993] B. L. Henke, E. M. Gullikson, and J. C. Davis, Atomic Data and Nuclear Data Tables 54, 181 (1993).
  * Chantler [1995] C. T. Chantler, Journal of Physical and Chemical Reference Data 24, 71 (1995).
  * Hunter and Tousey [1964] W. R. Hunter and R. Tousey, Journal de Physique 25, 148 (1964).
  * Baltazar-Rodrigues and Cusatis [2001] J. Baltazar-Rodrigues and C. Cusatis, Nucl. Instrum. Methods. Phys. Res. B: Beam Interactions with Materials and Atoms 179, 325 (2001).
  * Hildebrandt  _et al._ [1973] G. Hildebrandt, J. D. Stephenson, and H. Wagenfeld, Z. Naturforsch. Teil A 28, 588 (1973).
  * Macfarlane  _et al._ [1957] G. G. Macfarlane, T. McLean, J. E. Quarrington, and V. Roberts, Phys. Rev. 108, 1377 (1957).
  * Cardona  _et al._ [1970] M. Cardona, W. Gudat, B. Sonntag, and P. Y. Yu, in  _Proceedings of the Tenth International Conference on the Physics of Semiconductors_ , xiv, Vol. xiv, edited by F. Stern, J. C. Hensel, and S. P. Keller (U.S. Atomic Energy Commission, Division of Technical Information, U.S. Atomic Energy Commission, Division of Technical Information, 1970) Chap. Optical Absorption of Semiconductors from 15 to 170 eV, p. 209.
  * Efimov and Persson [1968] O. Efimov and E. Persson, Fiz. Tverd. Tela, 10, 2234 (1968).
  * Henke  _et al._ [1988] B. L. Henke, J. C. Davis, E. Gullikson, and R. P. ”, A preliminary report on x-ray photoabsorption coefficients and atomic scattering factors for 92 elements in the 10-ev - 10000-ev region (1988).
  * Philipp and Taft [1959] H. R. Philipp and E. A. Taft, Phys. Rev. 113, 1002 (1959).
  * Dash. and Newman [1955] W. C. Dash. and R. Newman, Phys. Rev. 99, 1151 (1955).
  * Potter [1997] R. F. Potter, in  _Handbook of Optical Constants of Solids_, edited by E. D. Palik (Academic Press, Burlington, 1997) pp. 465 – 478\.
  * Marton and Toots [1967] L. Marton and J. Toots, Phys. Rev. 160, 602 (1967).
  * Aspnes and Studna [1983] D. E. Aspnes and A. A. Studna, Phys. Rev. B 27, 985 (1983).
  * Berger  _et al._ [2010] M. J. Berger _et al._ , Xcom: Photon cross section database (version 1.5), https://www.nist.gov/pml/xcom-photon-cross-sections-database (2010), national Institute of Standards and Technology, Gaithersburg, MD.
  * Migdal [1941] A. B. Migdal, J. Phys. 4, 449 (1941).
  * Liu  _et al._ [2020] C. P. Liu, C.-P. Wu, H.-C. Chi, and J.-W. Chen, Phys. Rev. D 102, 121303 (2020), arXiv:2007.10965 [hep-ph] .
  * Dolan  _et al._ [2018] M. J. Dolan, F. Kahlhoefer, and C. McCabe, Phys. Rev. Lett. 121, 101801 (2018), arXiv:1711.09906 [hep-ph] .
  * Kahn  _et al._ [2021] Y. Kahn, G. Krnjaic, and B. Mandava, Phys. Rev. Lett. 127, 081804 (2021).
  * Knapen  _et al._ [2021] S. Knapen, J. Kozaczuk, and T. Lin, Phys. Rev. Lett. 127, 081805 (2021).
  * Liang  _et al._ [2020] Z.-L. Liang, L. Zhang, F. Zheng, and P. Zhang, Phys. Rev. D 102, 043007 (2020), arXiv:1912.13484 [cond-mat.mes-hall] .
  * Pepin [2016] M. D. Pepin,  _Low-Mass Dark Matter Search Results and Radiogenic Backgrounds for the Cryogenic Dark Matter search_ , Ph.D. thesis, Minnesota U. (2016).
  * Fascione [8 08] E. Fascione,  _Measurement of Cosmogenic Production Rates in Ge with CDMSlite_ , Master’s thesis, Queen’s U., Kingston (2018-08).
  * Agnese  _et al._ [2019b] R. Agnese _et al._ (SuperCDMS), Astropart. Phys. 104, 12 (2019b), arXiv:1806.07043 [physics.ins-det] .
  * Barker [2018] D. R. Barker,  _SuperCDMS Background Models for Low-Mass Dark Matter Searches_ , Ph.D. thesis, Minnesota U. (2018).
  * Agostinelli  _et al._ [2003] S. Agostinelli _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 506, 250 (2003).
  * Allison  _et al._ [2016] J. Allison _et al._ , Nucl. Instrum. Methods. Phys. Res. A: Accelerators, Spectrometers, Detectors and Associated Equipment 835, 186 (2016).
  * Allison  _et al._ [2006] J. Allison _et al._ , IEEE Transactions on Nuclear Science 53, 270 (2006).
  * Lindhard [1964] J. Lindhard, Phys. Lett. 12, 126 (1964).
  * Chavarria  _et al._ [2016] A. E. Chavarria _et al._ , Phys. Rev. D 94, 082007 (2016), arXiv:1608.00957 [astro-ph.IM] .
  * Collar  _et al._ [2021] J. I. Collar, A. R. L. Kavner, and C. M. Lewis, Phys. Rev. D 103, 122003 (2021).
  * Agnese  _et al._ [2018d] R. Agnese _et al._ (SuperCDMS), Appl. Phys. Lett. 113, 092101 (2018d), arXiv:1805.09942 [physics.ins-det] .
  * Yellin [2002] S. Yellin, Phys. Rev. D 66, 032005 (2002).
  * Yellin [2007] S. Yellin, Extending the optimum interval method (2007), arXiv:0709.2701 .
  * Cowan [1998] G. Cowan,  _Statistical data analysis_ (Oxford university press, 1998).
  * James and Roos [1975] F. James and M. Roos, Computer Physics Communications 10, 343 (1975).
  * iminuit team [2018] iminuit team, iminuit – a python interface to minuit, https://github.com/scikit-hep/iminuit (2018), accessed: 2018-03-05.
  * Alkhatib  _et al._ [2021] I. Alkhatib _et al._ (SuperCDMS), Light Dark Matter Search with a High-Resolution Athermal Phonon Detector Operated Above Ground (2021).
  * Aguilar-Arevalo  _et al._ [2020] A. Aguilar-Arevalo _et al._ (DAMIC), Phys. Rev. Lett. 125, 241803 (2020), arXiv:2007.15622 [astro-ph.CO] .
  * Collar [2018] J. I. Collar, Phys. Rev. D 98, 023005 (2018).
  * Agnes  _et al._ [2018] P. Agnes _et al._ , Phys. Rev. Lett. 121, 10.1103/physrevlett.121.081307 (2018).
  * Akerib  _et al._ [2019] D. S. Akerib _et al._ (LUX Collaboration), Phys. Rev. Lett. 122, 131301 (2019).
  * Armengaud  _et al._ [2019] E. Armengaud _et al._ (EDELWEISS Collaboration), Phys. Rev. D 99, 082003 (2019).
  * Armengaud  _et al._ [2022] R. Armengaud _et al._ (EDELWEISS Collaboration), Phys. Rev. D 106, 062004 (2022).
  * Aprile  _et al._ [2019] E. Aprile _et al._ (XENON), Phys. Rev. Lett. 123, 241803 (2019), arXiv:1907.12771 [hep-ex] .
  * Liu  _et al._ [2019] Z. Liu _et al._ (CDEX Collaboration), Phys. Rev. Lett. 123, 161301 (2019).
  * Liu  _et al._ [2022] Z. Z. Liu _et al._ (CDEX Collaboration), Phys. Rev. D 105, 052005 (2022).
  * l. Rudnick  _et al._ [2003] R. l. Rudnick, H. D. Holland, and K. K. Turekian, Treatise on Geochemistry , 659 (2003).
  * Emken  _et al._ [2019] T. Emken, R. Essig, C. Kouvaris, and M. Sholapurkar, JCAP 1909 (09), 070, arXiv:1905.06348 [hep-ph] .
  * Emken and Kouvaris [2018b] T. Emken and C. Kouvaris, DaMaSCUS-CRUST: Dark Matter Simulation Code for Underground Scatterings - Crust Edition (2018b), ascl:1803.001 .
  * Möller [2014] D. Möller,  _Chemistry of the Climate System_ (De Gruyter, 2014).
  * Wang  _et al._ [2018] H. S. Wang, C. H. Lineweaver., and T. T. Ireland, Icarus 299, 460 (2018), arXiv:1708.08718 [astro-ph.EP] .
  * Dziewonski and Anderson [1981] A. M. Dziewonski and D. L. Anderson, Physics of the Earth and Planetary Interiors 25, 297 (1981).