Light Dark Matter Search with a High-Resolution Athermal Phonon Detector Operated Above Ground

We present limits on spin-independent dark matter-nucleon interactions using a $10.6$ $\mathrm{g}$ Si athermal phonon detector with a baseline energy resolution of $\sigma_E=3.86 \pm 0.04$ $(\mathrm{stat.})^{+0.19}_{-0.00}$ $(\mathrm{syst.})$ $\mathrm{eV}$. This exclusion analysis sets the most stringent dark matter-nucleon scattering cross-section limits achieved by a cryogenic detector for dark matter particle masses from $93$ to $140$ $\mathrm{MeV}/c^2$, with a raw exposure of $9.9$ $\mathrm{g}\cdot\mathrm{d}$ acquired at an above-ground facility. This work illustrates the scientific potential of detectors with athermal phonon sensors with eV-scale energy resolution for future dark matter searches.

exposure of 9.9 g d acquired at an above-ground facility. This work illustrates the scientific potential of detectors with athermal phonon sensors with eV-scale energy resolution for future dark matter searches.
Introduction.-Numerous observations have shown that the majority of the Universe is composed of nonluminous matter [1][2][3]. The weakly-interacting massive particle (WIMP) [4] has long been a favored candidate for this dark matter (DM). However, direct detection experiments have ruled out a significant portion of the most compelling WIMP parameter space [5][6][7], which has motivated both theoretical and experimental exploration of alternative DM models [8]. In particular, light dark matter (LDM) with a mass in the keV/c 2 to GeV/c 2 range and coupling to Standard Model particles via a new force mediator provides a well-motivated alternative to the WIMP hypothesis [9][10][11]. While recent LDM searches have focused on DM-electron interactions [12][13][14][15][16], detectors with eV-scale energy thresholds can also be used to study LDM via DM-nucleon interactions.
We present results from a DM search with a new Cryogenic PhotoDetector (CPD) featuring an athermal phonon sensor with a baseline energy resolution of σ E = 3.86 ± 0.04 (stat.) +0. 19 −0.00 (syst.) eV. Although this device was designed for active particle identification in rare event searches [17], such as for neutrinoless doublebeta decay [18,19] and DM, the excellent energy resolution motivated its use as a DM detector itself. As a combined effort of the SuperCDMS and CPD collaborations, a DM search was carried out with 9.9 g d of raw exposure from Sept. 9 th to 10 th 2018. The data were acquired at the SLAC National Accelerator Laboratory in a surface facility of ∼ 100 m in elevation. We discuss data acquisition techniques, device performance, and the results of an exclusion analysis for spin-independent DMnucleon interactions.
Experimental setup.-The CPD substrate is a 1 mm thick Si wafer with a radius of 3.81 cm and a mass of 10.6 g.
It is instrumented on one side with ∼ 1000 Quasiparticle-trap-assisted Electrothermal feedback Transition-edge sensors (QETs) [20,21] distributed over the surface and connected in parallel to a single readout channel. The opposite side of the wafer is unpolished and not instrumented. The distributed channel results in minimal position dependence and fast collection of athermal phonons, which reduces inefficiency due to effects such as athermal phonon down-conversion [22,23]. The eV-scale baseline energy resolution was achieved in part because of the relatively low QET critical temperature of 41.5 mK with a nominal bath temperature of 8 mK.
A collimated 55 Fe source was placed facing the noninstrumented side. The electron capture decay provides Mn K α and K β X-ray lines at 5.9 and 6.5 keV, respectively, for in situ calibration [24]. A 38 µm layer of Al was placed in front of the collimator to attenuate the rate of these photons and provide an additional calibration line at 1.5 keV from Al fluorescence [25].
For the sensor readout, a direct-current superconducting quantum interference device (SQUID) array-based amplifier was used, similar in design to the one described in Ref. [26]. Because of project time constraints and large cosmogenic backgrounds, the DM search was limited to 22 hrs. Data were acquired over this period using a field-programmable gate array (FPGA) triggering algorithm based on the optimal filter (OF) formalism [27,28]. Throughout the exposure, randomly triggered samples of the baseline noise were acquired ("in-run random triggers"), which allowed us to observe any changes in the noise over the course of the search and to calculate and monitor the baseline energy resolution.
The trigger threshold was set at 4.2 σ above the normally-distributed baseline noise level, corresponding to 16.3 eV after calibration. The phonon-pulse template used for the FPGA triggering algorithm was a doubleexponential pulse with a rise time of τ r = 20 µs and a fall time of τ f = 58 µs. The rise time was taken from the expected collection time of athermal phonons, and the fall time was taken from the thermal response time of the QET estimated from a measurement of the complex admittance [20]. Each of these time constants was confirmed by a nonlinear least squares fit to nonsaturated pulses. Although Ref. [17] discusses the existence of extra fall times, their effect on the OF amplitude measured for each event is negligible. Before starting the DM search, a separate, small subset of random triggers was collected. After removing data contaminated by effects such as elevated baselines and phonon pulses, the noise spectrum used by the FPGA algorithm was generated from these random triggers.
For overlapping triggered pulse traces, the triggering algorithm was set to save a trace centered on the pulse with the largest OF amplitude. We note that the FPGA triggering algorithm acted on a trace that was downsampled by a factor of 16, from the digitization rate of 625 kHz to 39 kHz. Additionally, the FPGA triggering algorithm considered only 26.2 ms of the total 52 ms long time trace saved for each triggered pulse trace ("event"). Because of these factors, the energy resolution of the FPGA triggering algorithm is not as good as can be achieved by reconstructing event energies using an offline OF, as described in the following sections.
While the FPGA-based OF was used to trigger the experiment in real time, we ultimately used an offline algorithm to reconstruct event energies, where we again used the OF formalism. For this offline OF, we were able to use a single noise spectrum computed from the in-run random triggers to represent the entire data set because there was negligible time variation of the noise over the 25 1. A zoomed-in portion of an example event within the analysis ROI. The raw pulse (gray) is compared to the offline optimal filter result (blue), the pulse template scaled by the fit result (black dashed), the FPGA filter result (red with dots), and the FPGA trigger threshold (black dotted). The offline and FPGA optimal filters are highly correlated, but not exactly the same, with corresponding energy estimates for this event of 187 eV and 179 eV, respectively. The offset between the optimal filters and the raw pulse is an artifact of the filters, as they were set up to determine the time of the beginning of the pulse, as opposed to the maximum of the pulse.
course of the full exposure. Pulse amplitudes and start times were reconstructed using the same phonon-pulse template as in the FPGA triggering algorithm. Thus, there are two different pulse amplitudes for each eventone from the FPGA triggering algorithm and one from the offline OF. In Fig. 1, we compare the different energy estimators for a representative event.
This detector was optimized for maximum energy sensitivity at low energies and does not have a large enough dynamic range to observe the calibration lines without nonlinear effects from saturation of the QETs. The nonlinearity is minimal within our region-of-interest (ROI), which is below 240 eV. Above the ROI, the fall time of the pulses increases monotonically with energy, which can be explained by effects of local saturation. Localized events can saturate nearby QETs to above the superconducting transition, while QETs far from the event stay within the superconducting transition. Because this is a single-channel device, the saturated and unsaturated QETs are read out in parallel and thus effectively combine into a single phonon pulse with an increased fall time. In order to correct out the saturation effects within the ROI, we follow the calibration method as outlined in Ref. [17]. That is, the energy removed by electrothermal feedback (E ETF ) [20] is saturation-corrected using an exponential model, and the OF-based energy estimators are converted to units of energy via a linear fit to the calibrated E ETF within the ROI. With this method, there is an asymmetric systematic error in the baseline energy Measured energy spectrum in the DM-search ROI for the full exposure after application of the quality cuts. The data have been normalized to events per gram per day per eV and have been corrected for the event-selection efficiency, but not the trigger efficiency. The inset shows the calibrated EETF spectrum up to 7 keV, noting the locations of the different spectral peaks. The known values of the dashed lines are 1.5, 5.9, and 6.5 keV for the Al fluorescence (pink), 55 Fe Kα (blue), and 55 Fe K β (cyan) lines, respectively. The two dotted gray lines between 4 and 5 keV in calibrated EETF are the Si escape peaks [29].
resolution, for which the upper bound corresponds to the value achieved when calibrating E ETF linearly to the Al fluorescence line as opposed to the aforementioned exponential model (the lower bound). In Fig. 2, we show the differential rate spectrum of the calibrated offline OF amplitude, with the inset showing the differential rate spectrum for the calibrated E ETF .
Data selection and efficiency.-We make our final event selection with a minimal number of selection criteria ("cuts") to remove poorly reconstructed events without introducing energy dependence into the selection efficiency. This approach helps to reduce the complexity of the analysis and thus avoid introduction of systematic uncertainties. We apply two data-quality cuts: a prepulse baseline cut and a chi-square cut.
We define the event baseline as the average output in the prepulse section of each event, which is the first 25.6 ms of each trace. Large energy depositions have a long recovery time, which may manifest itself as a sloped baseline for subsequent events. Our trigger has reduced efficiency for any low-energy events occurring on such a baseline. We expected roughly 10% of the events to sit on the tail of a high energy event in part because of the high muon flux at the surface of ∼ 1 muon/cm 2 /min [30]. The baseline cut is performed by binning the data across the search in 400 s long bins and removing from each bin the 10% of events that have the highest baseline.
The chi-square cut is a general cut on our goodness-offit metric, for which we use the low-frequency chi-square χ 2 LF calculated from the offline OF fit [28]. This metric is similar to the χ 2 from the offline OF fit, but we exclude frequencies over f cutoff from the sum. This truncation allows us to remove sensitivity to superfluous degrees of freedom outside of our signal band from the chi-square, thereby reducing both the expected mean and the expected variance of the chi-square distribution. In this analysis, we used f cutoff = 16 kHz because the rise and fall times of our expected pulse shape correspond to frequencies of 8.0 kHz and 2.7 kHz respectively. The pulseshape variation within the DM-search ROI is minimal; this leads to a chi-square distribution that is largely independent of reconstructed event energy within this range. This in turn allows us to set an energy-independent cutoff value for χ 2 LF . Our measured events cannot be used directly to measure the signal efficiency of the chi-square cut because they include some that are not representative of the expected DM signal, e.g. vibrationally-induced events, electronic glitches, pileup events, etc. Therefore, we created a pulse simulation by adding noise from the in-run random triggers to the pulse template, systematically scaling the latter over the range of energies corresponding to the DM-search ROI. We then process and analyze the simulated data in the same way as the DM-search data. In this case, the passage fraction of the chi-square cut, which has an energy-independent value of 98.53 ± 0.01%, represents the cut's efficiency.
We do not apply any other cuts to the DM-search data. The total signal efficiency is thus 88.7% and is independent of energy. A variation of the cut values within reasonable bounds was found to have no significant impact on the experimental sensitivity.
Signal model.-In our DM signal model for spinindependent nuclear-recoil interactions [31], we use the standard astrophysical parameters for the dark matter velocity distribution [32][33][34]: a velocity of the Sun about the galactic center of v 0 = 220 km/s, a mean orbital velocity of the Earth of v E = 232 km/s, a galactic escape velocity of v esc = 544 km/s, and a local DM density of ρ 0 = 0.3 GeV/cm 3 . To take into account the trigger efficiency, we convolve the differential rate with the joint probability density function relating our two energy estimators, including the effects of the applied cuts. The signal model, which includes the estimated trigger efficiency, is given by where E 0 is the true recoil energy, E ′ is the recoil energy measured by the offline OF, E T is the recoil energy measured by the FPGA triggering algorithm, δ is the trigger threshold set on the FPGA triggering algorithm, ε is the efficiency of the two quality cuts and two cuts that are applied to simulated data (as described in the following paragraphs), Θ represents the trigger threshold cut (a Heaviside function), and P (E ′ , E T |E 0 ) is the probability to extract E ′ and E T using the two energy reconstruction algorithms given the true recoil energy E 0 . For the efficiency ε, we have generalized its form to be a function of energy, knowing that the baseline and chi-square cuts themselves are energy independent. The heat quenching factor (the ratio of heat signals produced by nuclear and electron recoils of the same energy that accounts for effects such as displacement damage) has been assumed to be unity for this work. Though measurements of the heat quenching factor have not been made for Si, similar work has been undertaken for Ge, where the heat quenching factor was shown to be very close to unity [35,36]. The model in Eq. (1) was evaluated numerically, taking advantage of our pulse simulation. The pulse simulation includes a software simulation of the FPGA triggering algorithm, which had the same output as the hardware version when run on the DM-search data. With this simulation of the FPGA triggering algorithm, we can use the pulse simulation to determine P (E ′ , E T |E 0 ) directly. Low pulse height events may have their OF energy estimate affected by a shift of the start time estimate, but the simulation automatically takes this effect into account.
We also added two cuts to the simulated data only: a confidence ellipse cut and a trigger time cut. The confidence ellipse cut removed any events with an energy estimator value outside of the 99.7% confidence ellipse, which is defined by the covariance matrix of our two energy estimators for zero-energy events. This cut was implemented to exclude the possible scenario of calculating a finite sensitivity to zero-energy DM, which would be a nonphysical result. The trigger time cut removed events that were not within 29 µs-half of a fall time of a pulseof the true event time, as determined by the energy-scaled pulse template. This cut ensured that the triggering algorithm was able to detect the signal added, as opposed to a large noise fluctuation elsewhere in the trace. These two cuts required knowledge of the true energy of the pulse-they cannot be applied to the data, but can be applied to the signal model-and helped to ensure that our signal modeling was conservative. In adding each of these cuts, we reduced our signal efficiency estimate, which necessarily biased the results in the conservative direction.
Results.-The objective of this DM search was to set conservative limits on the spin-independent interaction of dark matter particles with masses below 1.5 GeV/c 2 . For the lower edge of the limit contour, we use the optimum interval (OI) method [37,38] with unknown background. For the upper edge of the limit contour, we use a modified version of the publicly available verne code [39], which uses a Poisson method to calculate the effects of overburden [40][41][42] on the DM signal. This code has been similarly used in Refs. [43][44][45]. For the overburden assumption, we include the 5 cm of Cu surrounding  3. The 90% C.L. limits on the spin-independent DMnucleon cross section as a function of DM mass for this work (solid red line), compared to results from other experiments [45][46][47][48][49][50][51]. For above-ground experiments with overburden calculations, the previously ruled out parameter space is shown as the gray shaded region, and the new parameter space ruled out from this search is shown as the red shaded region. For the Collar 2018 surface result, which uses a liquid scintillator cell operated at 1 • C, an overburden calculation would be useful for comparison to the upper edges of the various contours for the surface searches. We note that the systematic error in the baseline energy resolution changes the result within the error of the limit's line width, thus we only include the result from the 3.86 eV calibration. the detector, the shielding from the atmosphere, and the shielding from the Earth. Both limit-setting methods assume that the full measured event rate could be due to a DM signal and set the limits at the 90% confidence level (C.L.).
The results of the dark matter search are shown in Fig. 3, compared to other pertinent DM searches in the same parameter space [45][46][47][48][49][50][51]. For DM masses between 93 and 140 MeV/c 2 , these results provide the most stringent limits for nuclear-recoil DM signals using a cryogenic detector. For DM masses between 220 MeV/c 2 and 1.35 GeV/c 2 , they are the most stringent limits achieved in an above-ground facility. For these low DM masses, the large cross sections approach the level at which the Born approximation used in the standard DM signal model begins to fail [52]. However, in the absence of a generally accepted alternative model and to be comparable to other experiments (all of which also use the Born approximation in this regime), we decided to keep it in our signal model as well.
To estimate the systematic error in the limit contour, we compared the results obtained by calculating the signal model using eight different sets of pulse simulations. The variation in the limit was found to be on the order of 10% for DM masses below 200 MeV/c 2 for the lower edge and below 100 MeV/c 2 for the upper edge. Above these DM masses, the variation in each edge decreased to less than 1%. The O(10%) variation observed at the smallest DM masses is attributed to a greater uncertainty in the trigger efficiency for sub-threshold events, as opposed to events that are reconstructed with energies above threshold. In the limit shown in Fig. 3, we have taken the median of the limits calculated for the eight simulations at each DM mass. The 10% variation is not plotted, as it would not be visible in the figure.
In Fig. 4, we show the data spectrum for reconstructed energies below 40 eV and DM signal curves for various DM masses for a single pulse simulation, where the cross sections from the OI limit are used. The approximate location of the optimum interval is apparent for each dark matter mass.
In this search, we see an excess of events for recoil energies below about 100 eV, emerging above the roughly flat rate from Compton scattering of the gamma-ray background. This excess of events could be from an unknown external background or due to detector effects such as crystal cracking [53]. As other experiments have observed excess events in searches for low-mass nuclear-recoiling DM [45][46][47]54], understanding this background is of pivotal importance. Future studies will be devoted to this, and we are actively investigating this excess by operating this detector in an environment with substantially reduced cosmogenic backgrounds.
Conclusion.-Using a detector with σ E = 3.86 ± 0.04 (stat.) +0. 19 −0.00 (syst.) eV baseline energy resolution operated in an above-ground facility with an exposure of 9.9 g d, we probe parameter space for spin-independent interactions of DM with nucleons for dark matter particles with masses above 93 MeV/c 2 . The range from 93 to 140 MeV/c 2 was previously not accessible to cryogenic detectors. These results also set the most stringent limits for above-ground nuclear-recoil signals from dark matter for masses between 220 MeV/c 2 and 1.35 GeV/c 2 . This was achieved using a single readout channel composed of QETs distributed on a Si substrate, with a recoil energy threshold set at 16.3 eV.
The results of this work were accomplished despite the high background rates in our surface facility because of the excellent baseline energy resolution of the detector. We plan to operate this detector in an underground laboratory, where we expect to have a significantly lower Compton scattering background rate. This will allow further study of the excess events observed in the ROI, hopefully providing insight into the origin of the event rate that is limiting the results reported here.
These results also demonstrate the potential of athermal phonon sensors with eV-scale baseline energy resolution for future dark matter searches via DM-nucleon interactions. Because this detector has a large surface area relative to its small volume, it is not optimal for a DM search. The baseline energy resolution of such devices scales with the number of QETs, which itself is proportional to the instrumented area (assuming the same QET design used by the CPD) [55,56]. Thus, a decrease in the instrumented area, with an increase in volume, should lead to improvements in baseline energy resolution. Future work is planned to design detectors of volume ∼ 1 cm 3 , for which it is reasonable to expect roughly an order of magnitude improvement in baseline energy resolution through these geometric considerations alone. With improved baseline energy resolution comes a lower energy threshold, allowing a search for spin-independent DM-nucleon interactions for even lower DM masses and a clear path to surpassing the existing noncryogenic detector constraints on sub-100 MeV/c 2 DM interacting with nucleons.