Searching for low-mass dark matter particles with a massive Ge bolometer operated above-ground

TL;DR

This study demonstrates the first sub-GeV dark matter search using a germanium bolometer operated above ground, setting new limits on low-mass dark matter interactions despite high background noise.

Contribution

It presents the first sub-GeV dark matter limits with a germanium detector operated above ground, including novel constraints via the Migdal effect for particles between 45 and 150 MeV/c².

Findings

Established the most stringent above-ground limits on spin-independent interactions above 600 MeV/c².

Set new constraints on spin-dependent interactions below 1.3 GeV/c².

Excluded new parameter space for strongly interacting dark matter particles considering Earth-shielding effects.

Abstract

The EDELWEISS collaboration has performed a search for dark matter particles with masses below the GeV-scale with a 33.4-g germanium cryogenic detector operated in a surface lab. The energy deposits were measured using a neutron-transmutation-doped Ge thermal sensor with a 17.7~eV (RMS) baseline heat energy resolution leading to a 60~eV analysis energy threshold. Despite a moderate lead shielding and the high-background environment, the first sub-GeV spin-independent dark matter limit based on a germanium target has been achieved. The experiment provides the most stringent, nuclear recoil based, above-ground limit on spin-independent interactions above 600~MeV/c$^{2}$. The experiment also provides the most stringent limits on spin-dependent interactions with protons and neutrons below 1.3~GeV/c$^{2}$. Furthermore, the dark matter search results were studied in the context of Strongly Interacting Massive Particles, taking into account Earth-shielding effects, for which new regions of the available parameter space have been excluded. Finally, the dark matter search has also been extended to interactions via the Migdal effect, resulting for the first time in the exclusion of particles with masses between 45 and 150~MeV/c$^{2}$ with spin-independent cross sections ranging from $10^{-29}$ to $10^{-26}$~cm$^2$.

Figures

Equations

$$H(f_i)=h\frac{s^{\ast }(f_i)}{J(f_i)}{e}^{-j2\pi f_i{t}_M},$$

$$h={\left(\sum _i\frac{|s(f_i){|}^2}{J(f_i)}\right)}^{-1},$$

$${\chi }^2(a,t_0)=\sum _i\frac{|\tilde{v}(f_i)-a\tilde{s}(f_i)e^{-j2\pi t_0{f}_i}{|}^2}{J(f_i)},$$

$$\frac{\mathrm{d}R}{\mathrm{d}{E}_{\mathrm{out}}}=\eta (E_{\mathrm{out}}){\int }_0^{\infty }{P}_{\mathrm{OF}}(E_{\mathrm{out}}|E_{\mathrm{in}})\frac{\mathrm{d}R}{\mathrm{d}{E}_{\mathrm{in}}}\mathrm{d}{E}_{\mathrm{in}}.$$

$$P(X^{\ast })$$

$$P(X^{\ast })=\frac{N_s}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{|X^{\ast }{|}^2}{2{\sigma }^2}\right){\left[\mathrm{erf}\left(\frac{|X^{\ast }|}{\sqrt{2}\sigma }\right)\right]}^{(N_s-1)}.$$

$$\begin{array}{rl}{P}_{\mathrm{OF}}(X^{\ast }|E_{\mathrm{in}})& =\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(X^{\ast }-E_{\mathrm{in}}{)}^2}{2{\sigma }^2}\right){\left[\mathrm{erf}\left(\frac{|X^{\ast }|}{\sqrt{2}\sigma }\right)\right]}^{(N_s-1)}\\ & +\frac{(N_s-1)}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{|X^{\ast }{|}^2}{2{\sigma }^2}\right){\left[\mathrm{erf}\left(\frac{|X^{\ast }|}{\sqrt{2}\sigma }\right)\right]}^{(N_s-2)}\frac{1}{2}\left[\mathrm{erf}\left(\frac{|X^{\ast }|-E_{\mathrm{in}}}{\sqrt{2}\sigma }\right)-\mathrm{erf}\left(\frac{-|X^{\ast }|-E_{\mathrm{in}}}{\sqrt{2}\sigma }\right)\right].\end{array}$$

$$P(E_{\mathrm{out}}|E_{\mathrm{in}})=\eta (E_{\mathrm{out}})P_{\mathrm{OF}}(E_{\mathrm{out}}|E_{\mathrm{in}}).$$

$$\frac{\mathrm{d}R}{\mathrm{d}{E}_{\mathrm{out}}}=\eta (E_{\mathrm{out}}){\int }_0^{\infty }{P}_{\mathrm{OF}}(E_{\mathrm{out}}|E_{\mathrm{in}})\frac{\mathrm{d}R}{\mathrm{d}{E}_{\mathrm{in}}}\mathrm{d}{E}_{\mathrm{in}}.$$

Full text

I Introduction

Various cosmological observations indicate that $26\%$ of the energy density of the Universe is in the form of cold, non-baryonic, dark matter (DM) planck2018-vi . Weakly Interacting Massive Particles (WIMPs) are suitable cold DM candidates; they arise in extensions of the Standard Model of Particle Physics, such as Supersymmetry, and are naturally produced in the early Universe with the correct abundance (for reviews see e.g. Refs. Jungman ; Bertone ). WIMPs from the Milky Way’s dark matter halo can be detected directly on Earth, via the keV-scale recoils produced when they elastically scatter off nuclei Goodman:1984dc ; annual . In recent decades, significant advances have been made in the search for WIMPs in the GeV/c2- to TeV/c2-range that is natural for Supersymmetry xenon1t ; lux ; pandax . However, in the light of the absence of signal in that region there is an increasing interest in DM particles in the GeV/c2 and sub-GeV/c2 mass range Essig ; Cheung ; Hooper ; Falkowski ; Petraki ; Zurek ; Bertone:2018krk . These searches require experimental thresholds as low as a few tens of eV, a performance that can be attained by cryogenic detectors cresst ; cdmslite . A particular advantage of such detector technology is that the thermal signal is not affected by the strong quenching effects that tend to severely reduce the amplitude of ionization or scintillation signals at low energy.

This paper describes the results obtained by the EDELWEISS collaboration with a 33.4-g Ge detector demonstrating that such a device equipped with a neutron-transmutation-doped Ge (Ge-NTD) sensor NTD can reach the sensitivity to probe the sub-GeV domain. As a proof of the relevance of this technology, it is used in a dedicated EDELWEIS-Surf run devoted to an above-ground search for DM particles. Such a search is bound to be limited by the large background induced by cosmic-ray interactions, but has the advantage of being able to probe models beyond the simple WIMP paradigm by considering relatively large values for the DM-nucleon cross sections. For sufficiently large values, each DM particle will typically interact many times in the atmosphere, Earth and shielding before reaching an underground detector. A DM particle loses energy with each interaction and may therefore arrive at the detector with insufficient energy to be observed above threshold Starkman:1990nj ; Collar1992 ; Collar1993 ; Hasenbalg:1997hs ; Kouvaris:2014lpa ; Kouvaris:2015laa ; Bernabei:2015nia ; Kavanagh:2016pyr . Above-ground searches minimize this ‘stopping’ effect in the Earth and therefore provide good sensitivity to large DM-nucleon cross sections. While most current constraints on strongly interacting massive particles (SIMPs111We use the term ‘SIMP’ here to refer to DM with large scattering cross sections with ordinary matter. However, we note that the term SIMP may also be used for hidden sector DM with strong self-interactions Hochberg:2014dra ; Hochberg:2014kqa .) rely on a reanalysis of public data Albuquerque:2003ei ; Davis:2017noy ; Kavanagh:2017cru , we perform a dedicated search that fully takes into account the detailed detector response, presenting both the smallest and largest cross sections to which the experiment is sensitive.

We also present a search for WIMP-induced nuclear recoils which are accompanied by the ionization of an atomic electron Vergados:2004bm ; Moustakidis:2005gx ; Bernabei:2007jz . The total energy in the nuclear recoil and ionization is typically larger than what can be deposited by elastic nuclear recoils at a given DM mass Ibe:2017yqa . This ‘Migdal effect’ therefore allows one to probe lighter DM particles (for a fixed energy threshold) at the expense of a smaller event rate. Constraints using the Migdal effect have been reported previously for NaI crystal detectors Bernabei:2007jz and liquid Xenon detectors Dolan:2017xbu ; Akerib:2018hck . Here, we present the first Migdal limit using a Germanium target.

The paper is organized as follows. In Sec. 1, we outline the experimental setup. In Sec. III, we give details of the dark matter search, including data processing, detector calibration and data analysis. In Sec. IV, we present limits on weakly- and strongly-interacting DM, for both elastic nuclear recoils and the inelastic Migdal effect. The resulting exclusion regions are summarized in Fig. 6. Finally, we conclude in Sec. V.

II EDELWEISS-Surf

II.1 Detector

The detector prototype consists of a cylindrical high-purity Ge crystal of 20 mm diameter and 20 mm height, corresponding to a total mass of 33.4-g. The thermal sensor design has been optimized for enhanced heat energy resolution. It consists of a Ge-NTD of $2\times 2\times 0.5$ mm3, glued on the top surface of the crystal, weakly thermally coupled to the copper housing of the detector thanks to gold wire bonds connecting its electrodes to two gold pads on a Kapton tape. With a total Ge-NTD electrode surface of 2 mm2 this weak thermal link is about 2.1 nW/K which is sub-dominant with respect to the electron-phonon coupling of 6.7 nW/K ensuring that the detector properly integrates all of the heat signal pyle . The crystal is held by six PTFE clamps (three on each side) in order to ensure the mechanical constraints on all three axes of displacement and to minimize the stress due to PTFE elasticity at low temperatures. Unlike the usual EDELWEISS-III FID800 detectors edwtech , this detector prototype has only one heat channel and no ionization readout. Therefore, a discrimination between nuclear and electron recoils is not possible. However, as there is no electric field applied across the crystal, the detector acts as a true calorimeter measuring the deposited energy of the recoiling particle independently of its type (nuclear or electronic recoil). Quenching effects on the heat energy scale for nuclear recoils in Ge cryogenic detectors have been shown to be very small quenching-edw ; quenching-cdms , and are therefore neglected hereafter.

II.2 Experimental setup

The dark matter search has been performed in the dry dilution cryostat of the Institut de Physique Nucléaire de Lyon installed in a surface building with negligible overburden, see Sec. IV.2. The cryostat is a Hexadry-200 commercially available from Cryoconcept cryoconcept , which has been upgraded to reduce the vibration levels of the mixing chamber by mechanically decoupling the cold head of the pulse tube cryocooler from the dilution unit cryovib . The vibrations at the detector level were further mitigated with the use of a dedicated suspended tower suspendedtower . The latter consists in a 25 cm long elastic pendulum, attached to the 1 K stage by a Kevlar string and a stainless steel spring with an elastic constant of 240 N/m, holding the detector tower situated below the mixing chamber at 10 mK. The detector tower is thermally anchored to an intermediate holding structure, via supple copper braids, which also hosts the connectors for the detector readout. This suspended tower design reduces detector vibrations at the sub-$\mu$g/$\sqrt{\text{Hz}}$ level, with displacements in the order of a few nanometers (RMS) in all three axes, leading to substantial gains in energy resolutions as demonstrated in Ref. suspendedtower . The cold and warm electronics are those described in Ref. edwtech , with a first Bi-FET preamplifier stage at 100 K and a second stage amplifier at 300 K.

The cryostat is surrounded by a 10 cm thick cylindrical lead shield covering a solid angle of $\sim$70% of $4\pi$ around the detector. The materials used for the cryostat construction were not selected for low radioactivity, with the exception of the replacement of the standard glass fiber rods used by Cryoconcept by stainless steel ones, shown to have much less radioactive contamination.

Finally, an 55Fe calibration source was glued on the inner part of the detector’s copper housing and facing the crystal surface opposite to the side on which is glued the Ge-NTD.

III Dark Matter search

The dark matter search started two weeks after the mixing chamber reached its base temperature of 10 mK. During these first two weeks, the thermal response of the detector was studied, and its working point was optimized. The best heat energy resolution achieved was $~{}17.7$ eV (RMS), with the temperature of the suspended tower regulated at 17 mK and the Ge-NTD biased at 1 nA, leading to a steady state resistance of 3.4 M$\Omega$. After the optimization period, it was decided to record six days of data in these conditions from May 22nd until May 27th 2018. It was decided beforehand to blind a 24-hour long data period started at 5 pm on May 26th. The remaining 5-days worth of data were then used to both tune the analysis procedure and selection cuts, and to build a data driven background model, in order to derive dark matter constraints from a blind analysis.

III.1 Data processing

The data acquisition used here is the same as from the EDELWEISS-III experiment located at the Laboratoire Souterrain de Modane edwtech . In order to cancel common electronic noise sources and reduce microphonics, the voltage drop across the Ge-NTD is measured differentially and the current across it is modulated from positive to negative values, following a square wave function. A modulation frequency $f_{s}$ of 400 Hz was chosen as it resulted in the best achievable signal-to-noise ratio in the experimental conditions considered here. The data are recorded continuously at this effective sampling frequency $f_{s}$, in so-called stream mode, such that there is no online trigger, unlike the standard EDELWEISS-III data acquisition. Instead, pulse signals are identified offline thanks to a dedicated signal processing pipeline, described in the following, which optimally filters the data based on the frequency dependence of the observed signal and noise. In order to avoid any noise structures at frequencies below the analysis range, from 1 Hz to 200 Hz, a second-order Butterworth numerical filter with a cut-off frequency of 2 Hz has been applied to the entire data stream before any selection of noise traces and triggers. Unless otherwise stated, the remaining part of the data processing is based on this pre-filtered data stream.

III.1.1 Noise PSD estimation

The noise Power Spectral Density (PSD) of the heat channel for an hour of recorded data is recursively determined from a random sample of 1-s time traces uniformly selected throughout the entire data stream for that hour. After an initial selection, based on their RMS dispersion, the removal of time traces containing pulses proceeds recursively by computing their individual frequency-based $\chi^{2}$ with respect to the averaged PSD. The procedure stops once the mean $\chi^{2}$, from both the pre-filtered and raw data stream, from all selected time traces is compatible with the expected value of $N_{s}$ = 400, corresponding to the number of time samples per trace. An average noise PSD is thus determined for each of the 137 hours that comprises the entire data set. These PSDs, corrected for the 2 Hz filter gain, are overlaid on Fig. 1 (left panel). The small dispersion shows that, despite its surface operation, the noise is very stable over the entire 6-days long run.

The PSD plateaus at a value of 4 nV/$\sqrt{\text{Hz}}$ above 20 Hz, with very little electromagnetic pickup contributions. This value is very well explained by the quadratic addition of the Johnson noise of the Ge-NTD (1.8 nV/$\sqrt{\text{Hz}}$) and the current noise of the Bi-FET preamplifier. The slow rise in the noise level below 20 Hz, reaching a value of 6 nV/$\sqrt{\text{Hz}}$ at 2 Hz, is due to internal thermal fluctuation noise from the detector thermModel .

III.1.2 Offline trigger

The offline trigger is based on a match filtered approach trigger designed to maximize the signal-to-noise ratio in the estimation of the signal pulse amplitude. The filter $H(f_{i})$ considered hereafter is derived from the measured frequency dependence of the signal-to-noise ratio and is therefore defined as:

[TABLE]

where $s^{*}(f_{i})$ is the complex conjugate of the signal template (shown as the black dashed line in the left panel of Fig. 1), $J(f_{i})$ is the noise PSD in V2/Hz (black solid lines), $t_{M}$ corresponds to the time position of the pulse template maximum, and $h$ is a normalization constant that preserves the amplitude of the pulse signal such that

[TABLE]

where $f_{i}$ varies between -$f_{s}/2$ to +$f_{s}$/2. An optimally matching filter is determined from each of the recorded 137 hours. Their moduli $|H(f_{i})|$ are also shown as the red solid lines in the left panel of Fig. 1. As one can derive from the latter, only the lowest frequencies (below  50 Hz) are relevant to recover the observed pulse amplitude. The filter from Eq. 1 is applied to the data using the numerical procedure described in Ref. trigger . As an example, the right panel of Fig. 1 shows a 200 eV event prior (grey solid line) and post (red solid line) filtering. Also shown is the best pulse fitting solution following the event processing procedure described in Sec. III.1.3. Candidate events are identified when the filtered data exceed a given threshold level which has been defined in terms of a fixed number $n$ of the baseline energy resolution $\sigma_{OF}$, where $\sigma_{OF}=\sqrt{h}$. The value of $n$ was chosen such that the rate of noise induced triggers is significantly smaller than the rate of physical events, which is about $\sim$1.3 Hz. The dependence on the value of $n$ of the rate of noise induced triggers has been evaluated by simulating a 24-hour long stream of fake data using the observed noise PSD, without injecting any signal pulses, and applying the same triggering procedure as for real data. We chose $n=3$ as it resulted in an optimum between a low energy threshold and a reasonably low expected noise induced trigger rate of 0.15 Hz. With the observed value of $\sigma_{OF}$ of 18 eV, this corresponds to a trigger threshold of 55 eV shown as the blue dotted line in Fig. 1 (right panel).

The triggering procedure searches for candidate events in the filtered data stream, starting with the largest positive deviation from zero. An exclusion interval of $\pm$0.5 s is imposed around each pulse found, and the search is reiterated in the remaining data until no fluctuation greater than $n\sigma_{OF}$ is found. This energy ordering of the pulse finding algorithm affects the energy dependence of the triggering efficiency. For instance, the dead-time associated to the search for low-energy events is effectively greater than that associated to large pulses. A dedicated data-driven procedure, described in Sec. III.3.2, has been developed to fully take into account these effects in the determination of the efficiency as a function of energy. The resulting trigger rate is almost constant over the six days and approximately equal to 0.58 Hz, as shown in the right panel of Fig. 2.

III.1.3 Event processing

Each recorded time trace, corresponding to either triggered events or noise samples, is further processed to estimate its amplitude by minimizing the following $\chi^{2}$ function defined in the frequency domain as:

[TABLE]

where $\tilde{v}(f_{i})$ is the Fourier transform of the observed signal, $a$ is the amplitude of the unit-normalized signal template $s$, and $t_{0}$ corresponds to the starting time of the pulse. For triggered events, the value of $t_{0}$ was allowed to vary within a 20 ms time window centered around the pulse time found by the trigger algorithm. In addition to delivering a slightly more precise estimation of the pulse amplitude $a$, this processing step provides a $\chi^{2}$ value that quantifies the quality of the fit. This is used to reject pulses with shapes that are not consistent with the standard pulse, such as pile-ups and other categories of spurious events further detailed in Sec. III.3.1.

III.2 Detector calibration and stability

III.2.1 Calibration

The energy calibration of the reconstructed amplitudes $a$ was ensured by the use of a low-energy X-ray 55Fe source irradiating the bottom side of the Ge crystal, opposite to the Ge-NTD heat sensor, inducing an interaction rate of $\sim$0.3 Hz. The 55Fe source produces two lines corresponding to the $K_{\alpha}$ and $K_{\beta}$ lines of Mn at 5.90 keV and 6.49 keV, respectively. They are clearly visible on the left panel of Fig. 2, showing the calibrated energy as a function of time over the 137-hour acquisition period. The energy resolution of these peaks is 34 eV (RMS). There is sufficient statistics to calibrate each hour separately and to follow precisely any time evolution of the detector gain, defined as the voltage sensitivity per unit of deposited energy, and resolution. The red solid line on the same figure shows the variation of the heat sensor gain (in $\mu$V/keV) as a function of time deduced from the peak position. The gain fluctuations are remarkably small during the entire time period, with a mean value of 1.56 $\mu$V/keV and a maximal dispersion of $\sim 2$%.

III.2.2 Baseline resolution stability

As mentioned is Sec. III.1.3, the time traces used for the PSD determination, also called noise events, were also fitted with Eq. 3. However, in this case $t_{0}$ was arbitrarily fixed to the center of the 1-s time window, in order to estimate the amplitude distribution associated with noise events. The amplitudes are distributed according to a gaussian distribution with a standard deviation (RMS) $\sigma$ that corresponds to the baseline energy resolution. This quantity, measured hour by hour, is compared to the prediction of the optimal filter resolution ($\sigma_{OF}$) on the right panel of Fig. 2. The good agreement between these two estimates of the energy resolution validates the assumption that the noise in each frequency bin is independent. The average baseline energy resolution is 18 eV (RMS), with a $\sim 3$% overall decrease in $\sigma_{OF}$ over the six days of the search. The average value of $\sigma$ during the blinded data period is 17.7 eV (RMS). The stability of the noise is instrumental at maintaining an almost constant trigger rate, shown as the red dots in Fig. 2 (right panel), over the entire search period.

III.3 Data analysis

III.3.1 Selection cuts

Because of its above-ground operation and the moderate lead shield, the detector was exposed to a rather intense rate of high energy interactions in the bulk of the crystal as well as in the holding PTFE clamps and the Ge-NTD heat sensor. The pulse shape of events produced in these cases are very different. In total, the following four types of events were observed:

• •

Normal events: these correspond to the vast majority of observed events for which the incoming particle has interacted in the Ge crystal. They are characterized by a rise time of 6 ms and a decay time of 16 ms, as expected from our thermal model calculations thermModel .

• •

Fast events: these are induced by incoming particles impinging the Ge-NTD heat sensor which has a non-negligible volume of 2 mm3, and from internal radioactivity of the Ge-NTD. They are characterized by a rise time of 2 ms, faster than normal events as expected, and a decay time constant of 41 ms. This longer decaying time with respect to normal events could be explained by athermal phonons heating up the Au wire bonds from the weak thermal leak. This hypothesis is further strengthened by the fact that the amplitudes of these fast events are much smaller than anticipated from our theoretical calculations thermModel .

• •

Slow events: these could be produced by muon interactions in the holding clamps where $\mathcal{O}(1)$ MeV energy is deposited. Despite the very weak thermal coupling of these PTFE clamps to the Ge crystal, such high energy deposition can still produce some non-negligible rise in temperature of the crystal. These events are characterized by a mean rise time of 26 ms and a decay time constant of 115 ms, hence much slower than normal events. This population was later observed to be significantly reduced in detectors where the PTFE clamps have been replaced by a combination of sapphire balls with chrysocale sticks.

• •

*Spike * events: the data acquisition system used in this experiment (a former version of the one used by EDELWEISS-III edwtech ) suffers from random synchronization losses that result in occasional octal jumps in the data flow. These are well modelled by delta functions that can easily be discriminated from physical events.

To discriminate normal events from the other populations of events, including also pile-ups, cuts are performed on the goodness of fit parameter quantified by the $\chi^{2}_{\text{Normal}}$ value from Eq. 3, where the subscript ”Normal” refers to the use of the standard pulse template $s$ in that equation. Additional pulse shape related selection cuts have been designed based on the $\chi^{2}$ differences $\Delta\chi^{2}_{k}=\chi^{2}_{\text{Normal}}-\chi^{2}_{k}$, where $\chi^{2}_{k}$ corresponds to the value calculated by replacing the standard pulse template for normal events $s$ with the templates of the non-standard event population discussed above, i.e. $k$ stands for Fast, Slow and Spike. All cuts were optimized on the non-blinded data set prior to un-blinding.

Figure 3 shows the event distribution in the $\chi^{2}_{\text{Normal}}$ versus energy plane before and after applying the $\Delta\chi^{2}$ cuts (black and grey points, respectively). The colored bands, estimated from the dedicated pulse simulation discussed in Sec. III.3.2, correspond to the 90% C.L. confidence bands for the four different types of events. Normal events are expected to have an averaged $\chi^{2}_{\text{Normal}}$ value centered around $N_{s}=400$ while the other event categories exhibit a quadratic divergence with increasing energy. As expected, the population of events that is the most alike to the normal ones, and therefore the hardest to reject efficiently, are the fast events. After the combination of the various $\Delta\chi^{2}$ cuts and selecting events with $350<\chi^{2}_{\text{Normal}}<450$ we found that the survival probabilities of the three types of spurious events reach $10^{-2}$ above 60 eV for Slow and Spike events, and above 400 eV for Fast events. The additional event population observed centered at 100 eV and $\chi^{2}_{\text{Normal}}\sim$ 2000 corresponds to noise events where the first half of the pre-filtered time trace is affected by the tail of a previous pulse (mostly from 55Fe), occurring a few hundreds of ms before the start of its 1-s time window. Note that this particular population of pile-up events, as all the other types of pile-up events that appear in the gaps between the colored bands in the left panel of Fig. 3, is fully taken into account in the pulse simulation procedure described in Sec. III.3.2, and is very well rejected by the standard $\chi^{2}_{\text{Normal}}$ cut.

III.3.2 Efficiency

The trigger efficiency was evaluated using a dedicated pulse simulation where pulses are generated at random times throughout the entire real data streams. This procedure samples rigorously all possible baseline fluctuations, including those induced by tails of high-energy events, and other non-stationary behaviour, hence avoiding any possible selection-induced biases. It naturally accounts for live-time losses due to the physical event rate and comprehensively takes into account any systematic uncertainties or biases related to the processing pipeline. The following efficiency estimates, as well as the various simulated DM signals discussed in Sec. IV, were obtained by generating a total of 106 simulated events distributed over 1000 Monte Carlo iterations of injecting 1000 simulated pulses in the 24-hour long blinded data stream. This way, the simulated event rate is about two orders of magnitude lower than the real physical event rate, hence inducing negligible additional live-time losses. A full simulation of the DM signal is performed for each mass value in order to evaluate in each case the combined effect of the noise observed in the actual data stream, the trigger selection and the analysis cuts. To better illustrate and understand the detector performance, it is however useful to calculate an efficiency taken as the fraction of simulated events that survive all of these selection criteria, from a population of simulated events with energies uniformly distributed between 0 and 2.5 keV. These resulting efficiencies at different stages of the trigger and analysis are shown in the right panel of Fig. 3. The black and grey points correspond to the variation of the combined trigger and live-time efficiency as a function of the reconstructed and input energy, respectively. The sharp rise around 55 eV of the efficiency curves expressed as a function of the reconstructed energy is due to the fact that the amplitude estimates from the triggering algorithm and the pulse processing are very similar. The rise of the efficiency as a function of the true input energy is much smoother and is well described by an error function centered at 55 eV with a dispersion of 18 eV, as anticipated from the observed baseline energy resolution and the considered 3$\sigma$ trigger level. As discussed in Sec. III.1.2, the trigger efficiency rises smoothly with energy, because of the energy ordering of the algorithm. Hence, between 0.1 and 2.0 keV, the trigger and live-time efficiency rises from 30% to almost 40%. These large dead-time losses can be expected given the size of the $\pm 0.5$ s exclusion time window compared to the observed trigger rate of 0.58 Hz.

It has been tested that the dead time could be reduced by using shorter time traces, but at the cost of a slight deterioration of the energy resolution. As it was anticipated that this high-background DM search would not be limited by statistics, the emphasis was put on low thresholds and therefore energy resolution.

The right panel of Fig. 3 also shows the effect of the $\Delta\chi^{2}$ on the analysis efficiency. Those cuts have a significant effect on the signal below 200 eV, where all three background populations compete with the signal. Finally, the cut on the $\chi^{2}_{\text{Normal}}$ value has a more uniform effect as a function of energy, as shown as the red points on the same figure.

III.3.3 Observed energy spectrum

The energy spectrum recorded within the blinded day of DM search data is shown in Fig. 4. The spectrum is dominated at high energy by the calibration source lines at 5.90 keV and 6.49 keV. A continuous and relatively flat background of $8\times 10^{3}$ evt/kg/keV/day is observed between 500 eV and 8 keV222The quoted background rates are corrected for the final efficiency as a function of energy (see red curve in Fig. 3).. The data below 500 eV, better illustrated in the inset figure, can be described approximately by an exponentially decreasing spectrum with a characteristic slope of 25 eV reaching $10^{5}$ evt/kg/keV/day at 200 eV and culminating at $2\times 10^{7}$ evt/kg/keV/day at the 60 eV analysis energy threshold. While the flat component of $8\times 10^{3}$ evt/kg/keV/day corresponds to the expected gamma-ray background for a detector operated in a non-low-radioactivity cryostat with a moderate lead shield, the exponential rise at the lowest energies is not yet fully explained. Dedicated studies are ongoing to better assess their origin, such as adding electrodes to read out the ionization yield associated to these interactions, and modifying the detector holders.

Finally, the inset figure also shows the contribution from the noise induced triggers (red histogram) that has been evaluated using the noise-only simulated data stream discussed in Sec. III.1.2. This event population explains very well the observed spectrum below 70 eV and it has a negligible contribution to the total observed event rate above 80 eV. It should be noted that in order to derive conservative constraints on DM interactions, it was decided prior to un-blinding that this well anticipated and modeled background will not be subtracted.

IV Results

This section presents the limits on spin-independent interactions of DM particles with nucleons derived from the experimental spectrum of Fig. 4. In Sec. IV.1 the data are interpreted in the standard WIMP context where the signal is the kinetic energy of the nuclear recoil induced in the collision. In Sec. IV.2 we extend its interpretation to the case of SIMPs where Earth-shielding effects are taken into account. In Sec. IV.3 we consider the so-called Migdal effect where the DM-nucleus interaction is inelastic and simultaneously produces a nuclear and an electronic recoil, opening up access to an unexplored domain of cross sections for DM particles with masses well below 100 MeV/c2. Finally, these results are interpreted in terms of spin-dependent interactions of DM particles in Sec. IV.4.

IV.1 Weakly Interacting Spin-independent Dark Matter

The expected signal for the standard elastic DM interactions in the detector was simulated using the comprehensive pulse simulation, described in Sec. III.3.2, that takes into account all systematic effects related to the detector response and data analysis. A total of $10^{6}$ fake DM-induced events were simulated for each DM mass considered. Their amplitudes were drawn from the theoretical distribution of induced nuclear recoils calculated assuming standard spin-independent interactions and using the standard astrophysical parameters for the DM velocity distribution standardassumptions : a Maxwellian velocity distribution with an asymptotic velocity of $v_{0}=$  220 km/s, and a galactic escape velocity of $v_{\text{esc}}=$ 544 km/s. The local DM density at Earth position is assumed to be $\rho_{0}=$ 0.3 GeV/c2/cm3 and the lab velocity with respect to the DM halo is $V_{\text{lab}}=$ 232 km/s. Also, the loss of coherence at high momentum is taken into account using the standard Helm form factor.

The objective of the search was to establish upper limits on the DM-nucleon cross section for each considered DM particle mass, using Poisson statistics and assuming that all events observed in a given energy interval are signal candidates. Prior to un-blinding the DM search data, the optimal energy intervals for such a purpose were determined by maximizing for each DM particle mass the signal-to-noise ratio between the simulated signals and the background model extracted from the 5 days of non-blinded data333In the presence of a signal, this procedure overestimates the background and yields conservative bounds..

After the un-blinding, the numbers of counts in these intervals were extracted from the data shown in Fig. 4. In the left panel of Fig. 5, the derived 90% C.L. upper limits of signals for WIMPs with masses of 0.7, 2 and 10 GeV/c2 are compared to the experimental data. The background model, i.e. the average spectrum recorded in the 5 non-blinded days, is also shown.

Finally, the spin-independent WIMP-nucleon cross sections excluded at 90% C.L. as a function of the WIMP mass are shown as the solid red curve in Fig. 6, and compared with the other experimental results from Refs. xenon1t ; lux ; cdmslite ; cresst ; edwiii ; darkside ; nucleus ; newsg ; xqc ; Dolan:2017xbu ; Akerib:2018hck , and the so-called neutrino floor neutrinoBillard . The EDELWEISS-Surf result is the most stringent, nuclear recoil based, above-ground limit on spin-independent interactions above 600 MeV/c2.

IV.2 Strongly Interacting Dark Matter

Thanks to its above-ground operation, the present DM search can probe SIMPs that would escape detection in underground experiments as the DM particles would be stopped in the rock overburden before reaching the detectors. We therefore extend the data interpretation of Sec. IV.1 to determine for each mass a range of excluded cross sections that take into account the absorption of the SIMP flux in the material above the detector. Overly conservative SIMP limits can be obtained by including in the analysis only those SIMPs that reach the detector without scattering Hooper:2018bfw . Here, more stringent limits are obtained by fully taking into account the effect of scattering on the velocity distribution of the SIMP flux reaching the detector. This velocity distribution is calculated using the publicly-available verne code verne , introduced in Ref. Kavanagh:2017cru . It assumes continuous energy losses and straight-line trajectories for SIMPs travelling through the atmosphere, overburden and detector shielding Starkman:1990nj . Ref. Emken:2018run showed that this simplified formalism leads to constraints similar to more complete but computationally expensive Monte Carlo simulations (e.g. Emken:2017qmp ; Mahdawi:2017cxz ; Mahdawi:2017utm ; Mahdawi:2018euy ).

The SIMP flux calculation takes into account the variation of the direction of the mean DM flux Mayet:2016zxu over the course of the blinded EDELWEISS-Surf exposure (24 hours, starting 17h00 on 26th May 2018). It also accounts for the atmospheric overburden above the detector as well as shielding provided by the material in the building where the detector is located and the lead, steel and copper components in its vicinity. The dominant sources of stopping for downward-travelling particles are 20 cm (40 cm) of concrete from the ceiling (walls) and 10 cm of lead shielding which surrounds the detector in all directions, apart from an opening of around $50^{\circ}$ above the detector. Upward-travelling particles are almost entirely stopped by the Earth.

The escape velocity from the surface of the Earth is around 11 km/s. For DM particles at such low speeds, effects such as gravitational capture Gould:1987ww ; Mack:2007xj and gravitational focusing Kouvaris:2015xga may become important. These effects are not incorporated in the flux calculation. Instead, the DM velocity distribution is conservatively set to zero below $v_{\mathrm{cut}}=20$ km/s when calculating SIMP bounds.

Because of the very large values of cross sections involved and consequently large attenuation of the flux, the simulation of the SIMP signals corresponding to the upper bound of the excluded cross section contour requires samples many orders of magnitude larger than those required in the simple WIMP analysis of Sec. IV.1. As scaling up the simulated sample size from $10^{6}$ to $>10^{10}$ was not technically feasible for computational reasons, we developed an analytic model for the detector response based on the simulation of $10^{7}$ events with input energies ranging from 0 to 2.5 keV (see Appendix A). This model describes the probability $P_{\mathrm{OF}}(E_{\mathrm{out}}|E_{\mathrm{in}})$ of reconstructing an energy $E_{\mathrm{out}}$ given an initial energy $E_{\mathrm{in}}$ when applying the optimal filter algorithm of Sec. III.1.3. The observed spectrum of events $\frac{\mathrm{d}R}{\mathrm{d}E_{\mathrm{out}}}$ is thus given by:

[TABLE]

The measured efficiency as a function of output energy is $\eta(E_{\mathrm{out}})$, as shown by the red curve in the right panel of Fig. 3. The calculation of $P_{\mathrm{OF}}$ and the comparison of the analytic detector response with results of the pulse simulations is discussed in Appendix A.

Using the signal calculated in these simulations, the same statistical procedure described in Sec. IV.1 is applied to derive the 90% C.L. upper bounds on the excluded cross section interval as a function of SIMP mass, resulting in the red contours shown in Fig. 6. The upper bound reported in this work improves upon the high-cross section reach of the CRESST 2017 surface run nucleus (thin blue), as reported in Refs. Davis:2017noy ; Kavanagh:2017cru ; Emken:2018run . This improvement is driven in part by the longer exposure of the EDELWEISS-Surf run, which covers a full day. This includes periods when the mean direction of the DM flux (set by the Sun’s velocity) is perpendicular to the Earth’s surface, reducing the attenuation effect of the Earth and atmosphere.

Sensitivities to 100 MeV/c2-scale SIMPs with spin-independent cross section of the order of 10*-30* cm2 have been reported in Collar:2018ydf . These were derived based on the dark counts observed in a photo-multiplier tube coupled to a liquid scintillator cell operated in a shallow site with 6 meter water equivalent overburden. However these resulting constraints did not take into account Earth-shielding effects, which are particularly relevant for such low-mass dark matter particles with cross section values above $10^{-31}$ cm2, hence preventing them to be compared with the ones presented in this work.

IV.3 Migdal Search

As discussed in Sec. II.1, the detector acts as a true calorimeter with equal sensitivity to the energy deposited by nuclear and electronic recoils. In this section, we consider the case where the WIMP or SIMP interaction with the target atoms induces simultaneously a nuclear recoil and the ionization of an electron. The final state comprises the two types of recoils. This is true in the case of the so-called Migdal effect, whose rate was calculated numerically in Ref. Ibe:2017yqa . The calculations therein were performed for the case of isolated atoms and does not fully consider the band structure of the germanium semi-conductor which is particularly important for the valence electrons ($n=4$).

In the absence of detailed calculations, we chose to neglect the contribution from the ionization of electrons in the $n=4$ valence shell, as well as the much smaller contribution from excitation of electrons into the valence shell. The inner electrons $n\leq 2$ are too tightly bound to give an appreciable signal. Therefore, the only contribution considered here is that from the ionization of M-shell ($n=3$) electrons. Once ionized, electrons are not free (as in Ref. Ibe:2017yqa ), but instead populate the conduction band. However, the band gap in germanium is typically much smaller than the M-shell ionization energies ($\sim 0.74\,\mathrm{eV}$ compared to $35\text{--}170\,\mathrm{eV}$, respectively) and here we will also neglect this small correction. Electrons, radiation and nuclear recoils in that energy range have very short absorption lengths in germanium, and it can be considered that the detector will completely collect the energy from all contributions.

The same standard spin-independent DM-nucleus interactions are assumed as in the previous sections, and notably the ($\propto A^{2}$) dependence of the cross section arising from the coherent coupling to the whole nucleus. The Earth-shielding effects are taken into account as in Sec. IV.2. While the observable signal arises from inelastic Migdal events, the predominant stopping power comes from the elastic DM-nucleus collisions, and therefore only these are taken into account in the calculation of Earth-shielding effects.

The right panel of Fig. 5 shows the expected Migdal spectra for WIMPs with masses of 0.05, 0.1 and 1.0 GeV/c2 corresponding to the respective interaction cross section values excluded at 90% C.L. The comparison with the standard signal in the left panel for WIMPs with considerably larger masses clearly shows how the Migdal effect greatly enhances the sensitivity to very light WIMPs, albeit for significantly larger cross section values.

The resulting 90% C.L. exclusion region is shown in Fig. 6, bounded from below by the thick dashed red line. The constraints weaken rapidly for WIMPs with masses below 100 MeV/c2. For lighter WIMPs, the peak in the expected Migdal signal spectrum overlaps in large part with the observed rapid rise of the data below 100 eV. Nevertheless, limits are obtained for WIMP masses as low as 45 MeV/c2. Below this, the large cross sections required to give an observable signal lead to substantial stopping effects in the atmosphere and shielding, meaning that no constraints can be obtained anymore.

For comparison, we also show in Fig. 6 the Migdal limit obtained by the LUX collaboration Akerib:2018hck (purple dashed line) but we do not show the unofficial limits based on XENON1T as presented in Ref. Dolan:2017xbu .

IV.4 Spin-dependent Interactions

In this last section we present our result in the context of a spin-dependent coupling between the dark matter particle and the nucleus. The isotope ${}^{73}\mathrm{Ge}$ (with a natural abundance of 7.73%) is the only stable germanium isotope with non-zero nuclear spin ($J=9/2$), hence providing sensitivity to spin-dependent (SD) WIMP-nucleus interactions. For the spin expectation values, we assume $\langle S_{p}\rangle=0.031$ and $\langle S_{n}\rangle=0.439$ Klos:2013rwa , with the larger $\langle S_{n}\rangle$ value is coming from a single unpaired neutron. For the spin-structure functions, we take the average of the values reported in Table V of Ref. Klos:2013rwa which includes both 1- and 2-body currents.

We account for Earth-scattering effects, including only the SD scattering of the dark matter particles with nitrogen-14 in the atmosphere, which is the dominant source of stopping for surface-based SD searches Hooper:2018bfw . Nitrogen-14 has nuclear spin $J=1$ and we assume equal proton and neutron spin contents $\langle S_{p}\rangle=\langle S_{n}\rangle=1$. The results, presented as 90% C.L. excluded regions and limits, are shown in Fig. 7 for DM-neutron and DM-proton SD couplings on the left and right panels respectively. Our results, show as red lines and contours, are compared to those of other direct detection experiments shown as solid lines: LUX Akerib:2017kat (purple), XENON1T Aprile:2019dbj (green), PICO-60-II Amole:2019fdf (brown), CDMSLite Agnese:2017jvy (pink), and PANDAX-II Fu:2016ega (blue). The other shaded contours correspond to the SIMP analyses from the XQC rocket xqc ; Mahdawi:2018euy ; Hooper:2018bfw (black solid line), the RRS baloon experiment Hooper:2018bfw , and the CMB Gluscevic:2017ywp (grey contour with dashed line). As for the SI case, see Sec. IV.2, the above-ground dark matter searches reported in Collar:2018ydf ; Abdelhameed:2019szb do not take into account Earth-shielding effects which are particularly relevant for cross section values above $10^{-31}$ cm2, and are therefore not shown in Fig. 7.

The EDELWEISS-Surf results are exploring new parameter space of the SD couplings, on both neutrons and protons, for dark matter masses ranging from 500 MeV/c2 to 1.3 GeV/c2 considering the standard DM induced nuclear recoil signal. The Migdal effect, producing both a nuclear recoil and an ionised electron, extends our lower DM mass bound down to 400 MeV/c2 only in the SD-neutron case because of the dominating Earth-shielding effects happening at such high cross sections. In the case of WIMP-proton SD interactions, we obtain no limit from the Migdal effect due to the strong Earth-shielding effects at such high cross sections. Indeed, the combination of the small proton spin content of ${}^{73}\mathrm{Ge}$ and the low ionization probability means that the scattering rate in this case is extremely small, hence requiring a large cross section to be detectable (around $10^{-24}\,\mathrm{cm}^{2}$ at $1\,\mathrm{GeV}/c^{2}$). However, at such large cross sections, atmospheric stopping is too strong for a DM particle to reach our detector. We also note that at low masses ($\mathcal{O}(500\,\mathrm{MeV}/c^{2})$ both the Migdal WIMP-neutron and the standard WIMP-proton constraints are affected by Earth-stopping effects at both the upper and lower parts of the excluded regions.

V Conclusion

The EDELWEISS collaboration has searched for dark matter particles with masses between 45 MeV/c2 and 10 GeV/c2 using a 33.4-g germanium detector operated in a surface lab, and thus relevant in the search for Strongly Interacting Massive Particles (SIMPs). The energy deposits were measured using a Ge-NTD thermal sensor with a 17.7 eV (RMS) heat energy resolution leading to a 60 eV analysis threshold. This performance, combined with the nearly completely stationary behavior of the detector, led to the achievement of the first limit for the spin-independent interaction of sub-GeV WIMPs based on a germanium target. The experiment provides the most stringent, nuclear recoil based, above-ground limit on spin-independent interactions above 600 MeV/c2.

The search results were also interpreted in the context of SIMPs, taking into account the screening effect of the atmosphere and material surrounding the detector. The lower part of the excluded region for SIMPs corresponds to the previously mentioned WIMP limit and represent the most stringent constraint for masses above 600 MeV/c2. The upper part of the excluded region is limited by Earth-shielding effects: it probes the largest SIMP-nucleon cross sections of any direct detection experiment, excluding a value of $10^{-27}\,\mathrm{cm}^{2}$ for a 1 GeV/c2 WIMP444Note that CDMS-II provides the strongest SIMP constraints for DM masses above around $10^{8}$ GeV/c2 Albuquerque:2003ei ; Kavanagh:2017cru .. There are a number of complementary constraints on SIMP DM, including searches for DM annihilation to neutrinos Albuquerque:2010bt ; anomalous heating of the Earth Mack:2007xj ; Mack:2012ju ; heating of Galactic gas clouds Bhoonah:2018wmw ; Bhoonah:2018gjb ; and DM-cosmic ray interactions Cyburt:2002uw . These typically require additional assumptions about the properties of the DM particle, while the current constraints depend only on its large scattering cross section with nuclei and its interactions near the Earth.

The dark matter search has also been extended to interactions via the Migdal effect, resulting in the exclusion for the first time of particles with masses between 45 and 150 MeV/c2 with cross sections ranging from $10^{-29}$ to $10^{-26}$ cm2. These limits also take fully into account the modeling of Earth-shielding effects essential for obtaining accurate constraints for such large cross section values.

Finally, interpreted in terms of spin-dependent interactions with protons or neutrons, these results exclude new regions of the parameter space below masses of 1.3 GeV/c2. In this case, atmospheric stopping is significant and has a strong effect on the derived exclusion limits, in particular for WIMP-proton interactions and the Migdal effect.

The next steps for the EDELWEISS collaboration are twofold. Such high-performance detectors are currently running in the low-background environment of its underground facility edwtech at the Laboratoire Souterrain de Modane to better understand the origin of the observed events at low energy. Also, the detectors are equipped with electrodes in order to amplify the signal using the Neganov-Trofimov-Luke amplification Neganov ; Luke , a particularly enticing prospect for searches using the electron recoils produced by the Migdal effect. The level of performance achieved in this work also opens the possibility of a first experimental measurement of the Migdal effect using a neutron calibration source.

Eventually, this level of detector performance, achieved in an above-ground laboratory with a 30 g-scale massive bolometer, is also very promising in the context of a low-energy and high precision measurement of the Coherent Elastic Neutrino-Nucleus Scattering process Billard:2016giu .

Acknowledgements.

We thank Felix Kahlhoefer for helpful discussions about the Migdal effect in semi-conductor materials. BJK acknowledges funding from the Netherlands Organization for Scientific Research (NWO) through the VIDI research program “Probing the Genesis of dark matter” (680-47-532). Part of this work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. The EDELWEISS project is supported in part by the German Helmholtz Alliance for Astroparticle Physics (HAP), by the French Agence Nationale pour la Recherche (ANR) and the LabEx Lyon Institute of Origins (ANR-10-LABX-0066) of the Université de Lyon within the program “Investissements d’Avenir” (ANR-11-IDEX-00007), by the P2IO LabEx (ANR-10-LABX-0038) in the framework “Investissements d’Avenir” (ANR-11-IDEX-0003-01) managed by the ANR (France), and the Russian Foundation for Basic Research (grant No. 18-02-00159).

Appendix A Analytic Detector Response

This appendix describes the analytic model for the detector response which is used to set limits in the SIMP and Migdal effect analyses.

As described in Sec. III.1.3, the optimal filter (OF) algorithm is used to search for candidate events in the data stream and recover the pulse amplitude. The OF algorithm locks on the largest fluctuation away from zero in a given search window. Note that this could be a downward fluctuation, such that the reconstructed pulse amplitude may be negative. The sampling frequency is $400\,\mathrm{Hz}$ while the OF search window has a width of $20\,\mathrm{ms}$. If each sample were independent, this would correspond to a total of 8 samples. However, correlations in the noise mean that the effective number of independent samples $N_{s}$ in each search window is smaller.

We now calculate the probability distribution for the largest fluctuation (i.e. deviation from zero) among $N_{s}$ independent random variables $X_{i}$. We denote the probability distribution of each random variable $P_{i}(X_{i})$ and the largest fluctuation $X^{*}$. With this, we obtain:

[TABLE]

This is the probability that one of the random variables lies in the interval $X^{*}\rightarrow X^{*}+\mathrm{d}X^{*}$, while the remaining random variables $X_{k\neq i}$ have absolute values smaller than $X^{*}$.

We assume Gaussian noise and fix the resolution to $\sigma=18\,\mathrm{eV}$ (see right panel of Fig. 2). If there are no signal events in the OF search window, then the distribution of the largest fluctuation is given by:

[TABLE]

Identifying $X^{*}$ with $E_{\mathrm{out}}$, we can compare this distribution with the reconstructed energy of noise-only events, for which we find that $N_{s}=3$ provides a very good fit.

If a signal event with energy $E_{\mathrm{in}}$ appears in one of these $N_{s}=3$ independent samples, the resulting random variable will still be Gaussian-distributed with width $\sigma$ but now centred on $E_{\mathrm{in}}$ (the remaining two random variables are Gaussian-distributed with mean zero). Again, applying Eq. 5, we obtain the probability distribution for the largest fluctuation:

[TABLE]

At small input energies, this expression accounts for the fact that the largest fluctuation in the OF search window may not correspond to the signal pulse. At large input energies, we recover a standard Gaussian distribution for the output energy.

Finally, we apply the trigger efficiency $\eta(E_{\mathrm{out}})$ as a function of the reconstructed energy, given by the red curve in the right panel of Fig. 3. The final detector response is then:

[TABLE]

The observed spectrum of events is then obtained by a convolution:

[TABLE]

In Fig. 8, we compare the probability distribution for the output energy (at fixed input energy) estimated using pulse simulations (left panel) and using the analytic prescription (right panel). We use more than $10^{7}$ pulse simulations, which were performed in order to calculate the standard WIMP limits presented in Sec. IV.1. We include only those simulated events passing all cuts and present the results in bins of 2.5 eV in both $E_{\mathrm{in}}$ and $E_{\mathrm{out}}$. The analytic response function captures the main features of detector response seen in the pulse simulations.

We have also performed a more quantitative comparison. Binning the events in bins of 1 eV in input energy, we performed a standard Kolmogorov-Smirnov (KS) test Kolmogorov ; Smirnov in each bin. We compared the output energy distribution of pulse simulations with the distribution expected from the analytic model. We found that the two distributions are consistent, with $p$-values ranging from $10^{-3}$ to 1 over the 2000 bins as expected.

Finally, we have verified that the standard WIMP limits calculated using this analytic model are in good agreement with those obtained using the pulse simulations (as described in Sec. IV.1). We find that the two methods agree at the 10% level.