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

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


I. INTRODUCTION
An abundance of evidence suggests that most of the Universe is composed of nonluminous matter [1][2][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 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 χ , 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 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 χ ≪ 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][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 AE2 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 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.9kg • 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 characteristic 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 truncated at 544 km=s and renormalized [16].The local DM density, ρ χ , is assumed to be 0.3 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][19][20].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 [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.A velocity damping parameter, κ, is defined as: 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 , where the incoming SHM velocity distribution (f) has been transformed by the effects of shielding [22].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 κ 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 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 continuously 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 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 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 hemisphere, the majority of dark matter will come from above, which means this assumption will conservatively underestimate 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: Geological data for the Soudan region is provided in two auxiliary files: The SoudanRegion.csvfile contains the area and fraction of each rock type, and the elemental mass fractions for the chemical composition are described in RockChem.csv.
where r d is the distance from the center of the Earth to the SuperCDMS experiment.The average value of 6471 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 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 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.

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]: 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, where N T is Avogadro's number divided by the atomic mass.The minimum velocity required to induce a recoil of E γ , v min , is given by where E R is the nuclear recoil energy.

Photoelectric absorption cross section uncertainty
The atomic scattering function, f, in Eq. ( 5) is the sum of a real and complex portion, The components, f 1 and f 2 , are defined in Ref. [31] as: 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][34][35][36][37][38][39][40][41][42][43][44][45][46][47], a nominal value and uncertainty band for σ a were derived as a function of energy.Using the resulting values of σ a , jfj 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.calculations of jfj 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 simultaneously 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 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 [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 calculations 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: 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 =ðdE R dvÞ describes the incoming DM flux and scattering rate, 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: 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][53][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, jfj 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 jfj 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.
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, cosmogenic activation of the crystal, radiogenic Compton scattering, 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 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 71 Ge background.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, cosmicray 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 AE 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 distributed 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 surrounding 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][60][61] simulation of 210 Pb surface contamination on a tower of six germanium detectors was performed.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 independent 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 AE 16.6, N SH ¼ 22.3 AE 1.4, and N GC ¼ 23.5 AE 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 uncertainties 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 efficiency models have been extended to 25 keV ee following the procedure in Ref. [57].
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 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 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 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)], 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 AE1σ uncertainty band.production was cut off at eV, according to an independent measurement of the defect energy creation threshold [65].This is implemented with a hyperbolic tangent function, 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, AE1σ, and AE3σ 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] 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.

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 logarithmically 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.V B. 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.three types of terms.The first type is the overall normalization term, (L Poiss ), which allows total number of fitted events to fluctuate around the number of events in the dataset (N).The second type (L P1 χ;bb;sb , L P2 χ;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 P1 eff , L P2 eff ), the resolution model (L P1 res , L P2 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 covariance matrix is calculated from the normalization uncertainty of the individual components and the correlations between them [68].

Variable
This analysis used the MINUIT algorithm [69] via the iminuit [70] Python interface to maximize the log-likelihood 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 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 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: 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 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 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 GeV=c 2 .To test the sensitivity of the limit to yield modelling, we crosschecked 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.
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 bremsstrahlung channel limit to extend below the 0.4 GeV=c 2 reach of the LUX and XENON1T bremsstrahlung analyses.However, both experiments have larger exposures, 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 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 32 and 39 MeV=c 2 .

APPENDIX:
Parameters for the earth and atmosphere are listed in Table IV.The density of the Earth was taken from Ref. [81].

2
Studies based on simulation have shown slightly more bias toward the surface of the crystals due to self-shielding effects.

TABLE I .
Definition of variables in Eq. (1) to calculate the velocity damping parameter, κ.

TABLE II .
[25]ities and heights of the layers in the atmosphere model[25].