Constraints on low-mass, relic dark matter candidates from a surface-operated SuperCDMS single-charge sensitive detector

This article presents an analysis and the resulting limits on light dark matter inelastically scattering off of electrons, and on dark photon and axion-like particle absorption, using a second-generation SuperCDMS high-voltage eV-resolution detector. The 0.93 gram Si detector achieved a 3 eV phonon energy resolution; for a detector bias of 100 V, this corresponds to a charge resolution of 3% of a single electron-hole pair. The energy spectrum is reported from a blind analysis with 1.2 gram-days of exposure acquired in an above-ground laboratory. With charge carrier trapping and impact ionization effects incorporated into the dark matter signal models, the dark matter-electron cross section $\bar{\sigma}_{e}$ is constrained for dark matter masses from 0.5--$10^{4} $MeV$/c^{2}$; in the mass range from 1.2--50 eV$/c^{2}$ the dark photon kinetic mixing parameter $\varepsilon$ and the axioelectric coupling constant $g_{ae}$ are constrained. The minimum 90% confidence-level upper limits within the above mentioned mass ranges are $\bar{\sigma}_{e}\,=\,8.7\times10^{-34}$ cm$^{2}$, $\varepsilon\,=\,3.3\times10^{-14}$, and $g_{ae}\,=\,1.0\times10^{-9}$.

In a prior work [12], we undertook an above-ground search with a first-generation Si HVeV detector. Contemporaneously, the SENSEI Collaboration reported an underground search with Skipper CCDs [23]. Both works excluded new parameter space for light DM scattering and dark photon absorption in similar mass ranges, but did not report on the axioelectric coupling, which is most strongly constrained by astronomical observations [24][25][26][27]. In this work, we analyze a slightly larger aboveground exposure of 1.2 gram-days of a second-generation Si HVeV detector with the same dimensions but modified sensor design compared to [12], leading to an improved phonon energy resolution as good as σ E = 3 eV at the single-e − h + -pair energy (3 % charge resolution for a 100 V bias). Using signal models that include the contributions from charge carrier trapping and impact ionization, we report the constraints obtained from a blind analysis on χ scattering for DM masses from 0.5-10 4 MeV/c 2 , as well as V and ALP absorption for masses from 1.2-50 eV/c 2 .

II. EXPERIMENTAL SETUP
The data were acquired in a surface laboratory at Northwestern University (Evanston, IL), with the overburden and environmental radioactivity of a typical steelconcrete building. The detector is made of a 0.93 gram * yychang@caltech.edu † mwilson@physics.utoronto.ca high-purity Si crystal (1 × 1 × 0.4 cm 3 ). We clamped the detector between two printed circuit boards for thermal sinking and electrical connections. To reject correlated environmental noise, we installed an anticoincidence (veto) detector adjacent to the HVeV detector in the same light-tight copper housing mounted to the cold finger of an Adiabatic Demagnetization Refrigerator (ADR). More information about the detector setup and the infrared radiation shield is available in Ref. [14]. SuperCDMS HVeV detectors measure phonons created by particle interactions in the Si crystals with two distributed channels of Quasiparticle-trap-assisted Electrothermal-feedback TESs 1 (QETs) [14]. The QETs fabricated on this device have a superconducting transition temperature T c ≈ 65 mK. One QET channel is a square with a sensitive area of 0.5 cm 2 , and the other is a surrounding frame of equal area. Both are on the detector's top surface. On the bottom surface, an aluminum grid with 5 % surface coverage provides a uniform electric field between the two surfaces. The veto detector consists of a single TES on a thin Si wafer that is identical to the TES described in Ref. [29] but with T c ≈ 52 mK.
We cycled the ADR daily from 4 K to the base temperature and then regulated it at 50-52 mK during data taking to obtain a 10-12 hour/day hold time [14]. To induce Neganov-Trofimov-Luke (NTL) amplification [30,31], the aluminum grid was biased at V bias = 100 V while the QETs and detector housing were held at ground potential. At the start of each daily cycle, we set the operating point of each QET to ∼ 300 mΩ (about 45 % of its normal-state resistance) and recalibrated the detector using a 635-nm laser that was fiber-coupled from room temperature. Each QET was read out with a DC superconducting quantum interference device (SQUID) at 1 K operated in a flux-locked feedback loop, and the signals were digitized continuously at 1.51 MHz. The laser intensity was adjusted to achieve a mean number of e − h + pairs per pulse between 1 and 4, which produced enough events for calibration up to seven e − h + pairs per event. We also took a dedicated laser dataset in which we varied the crystal temperature but held the QETs at their nominal operating point; we used this dataset to reconstruct the temperature dependence of the QET responsivity.

III. DATA COLLECTION AND EVENT RECONSTRUCTION
A raw exposure of 3.0 gram-days was collected over 7 days during April-May of 2019. By partitioning the continuous-acquisition data stream into 10-second long intervals, we performed a three-stage blind analysis. The first second of each interval, i.e. 10 % of the data, was used to develop the analysis pipeline but was not included in the final spectrum. Data from seconds 2-3 of each interval were unblinded to verify that the analysis pipeline was indeed invariant under the presence of a larger statistical sample. Given that the initial unblinding satisfied this condition, we unblinded the remaining data and analyzed seconds 2-10 from each data partition, i.e. 90 % of the data defined as the DM-search data, to extract the final results.
To identify physics events, we triggered on pulses within the continuous-acquisition data stream offline. To issue triggers, we first applied a matched filter to the data stream using an exponential pulse template (20 µs rise time and 80 µs fall time) and then set a trigger threshold equivalent to ∼ 0.4 e − h + pairs for event identification. The event trigger time is the time at which the triggered pulse is at its maximum. After verifying that the two QET channels on the HVeV detector have equal gain, this trigger scheme was applied to the sum of the two channels' data streams and separately to the veto detector.
Pulse energy and time were reconstructed using an optimal filter (OF) algorithm [32,33]. The OF algorithm requires a pulse template and the noise power spectral density (PSD). We constructed the pulse template for the OF algorithm from the laser-calibration event pulses with high-frequency noise filtered out. We measured the noise PSD on an hourly basis to account for variations of the environmental noise, using the first 100 seconds of each one-hour data partition with triggered pulses removed. The pulse amplitude and time that minimize the frequency-domain χ 2 were determined within a time window of [− 678 µs, + 2034 µs] centered on the trigger. These amplitude and time quantities of the OF algorithm were also used to compute a time-domain χ 2 which was used later in the analysis.
Temperature fluctuations at the detector and small variations in the HV bias resulted in a small variation (< 1 %) in the detector gain. We used the quantized e − h + -pair peaks in the aforementioned temperaturecontrolled sample spectrum, as well as the daily lasercalibration spectra, to linearly correct for the temperature and voltage dependencies. We then corrected for nonlinearities in the detector response with a secondorder polynomial.
To calibrate the OF pulse amplitudes to energies we rescaled the e − h + -pair peaks assuming where n denotes the n-photon absorption peak, E γ = 1.95 eV is the laser photon energy, and e is the absolute value of the electron charge. Figure 1 (top panel) shows the resulting spectra from the DM-search and laser-calibration data. Both spectra show the data measured with a detector bias of 100 V after applying the live-time and dataquality cuts. The peak seen at ∼ 50 eV in the DM-search data is due to non-quantized events restricted to the outer QET channel [14]. Light gray-shaded regions on the left-and right-hand sides mark the energy ranges outside the region of interest; vertical lines correspond to the phonon energy En of the n-photon absorption peak (Eq. 1). The black curve is an example of a signal produced by electron-recoiling dark matter particles with a mass of 1 GeV/c 2 and form factor FDM ∝ 1/q 2 . This model assumes a Fano factor of F = 0.155, an impact ionization (II) probability of 2 %, and a charge trapping (CT) probability that varies from 0-15 % shown by the hatched region. The curve is scaled to the dark matterelectron cross sectionσe that sets the limit at the 2 nd e − h +pair peak. The bottom panel shows the binned efficiency data (Ei) (grey solid line), where the corresponding shaded region indicates the 1 σ statistical uncertainty in each bin. The red dashed curve is the efficiency curve, and the corresponding shaded region is the conservative efficiency uncertainty envelope, which accounts for the statistical and systematic uncertainties.

IV. DATA SELECTION
To ensure accurate event reconstruction, individual live-time intervals from the DM-search and lasercalibration data were discarded (cut) based on various criteria: (1) the ADR temperature to exclude data outside the range of the temperature calibration; (2) the pre-pulse baseline averaged over one-second intervals to reject periods of time when the detector was still recovering from a preceding high energy deposition; and (3) the trigger rate to remove bursts of non-DM triggers. The trigger-rate cut was not applied to the laser-calibration data. The data remaining after these live-time cuts de-fine the science exposure for this analysis, and yielded a DM-search exposure of 1.2 gram-days.
To reject poorly reconstructed events in the DM-search exposure, a set of event-by-event data-quality cuts were applied based on: (1) the difference between the OFdetermined pulse time and the event trigger time to reject noise triggers and pulses affected by pile-up events; (2) the event-by-event average pre-pulse baseline to ensure the detector is at a steady working condition before an event occurs; and (3) the energy-dependent frequencyand time-domain χ 2 quantities to remove pile-up events and baseline excursions that are unlikely to have been caused by DM-triggers. To define the cuts, we determined the nominal distributions of each parameter using the laser-calibration data and discarded events in the DM-search exposure exhibiting an excursion of > 3 σ in any of these parameters. Lastly, we rejected events with a coincident triggered event in the veto detector.
The cut efficiency as a function of phonon energy was determined using the laser-calibration data after applying the live-time cuts to pulses coincident with the laser trigger signal. Pile-up events that occurred within a laser-event trigger window were included as part of the efficiency calculation. The binned efficiency (E i ) is the fraction of events in the i-th bin that pass the quality cuts. We expect the efficiency to be smooth; however, our measurement of (E i ) shown in Fig. 1 (bottom panel) has both statistical fluctuations and systematic uncertainty. In order to avoid folding these effects into the final results, we fit a smooth function to (E i ) and assigned a conservative envelope that accounts for the statistical and systematic uncertainties. The systematic uncertainty is due to pile-up events that were not rejected by the live-time cuts, resulting in a misreconstruction of the energy. Although this envelope was propagated as part of the total experimental uncertainty in the final results, we verified that it is not the dominant source of uncertainty.
The energy region of interest (ROI) for this analysis is 50-650 eV. The lower bound guarantees inclusion of the full single-e − h + -pair peak at 100 V bias and a trigger efficiency consistent with unity. We set the upper bound at 650 eV to focus on the corresponding low-mass ranges where this analysis has competitive sensitivity.

A. DM Signal Models
The blinded DM-search data were analyzed to set exclusion limits on light DM χ scattering as well as dark photon V and axion-like particle (ALP) absorption. The DM models for χ, V , and ALPs are identical to those used in Ref. [12] and [34]. We set limits on the kinetic mixing parameter ε for V , the axio-electric coupling g ae for ALPs, and the effective DM-electron cross sectionσ e for χ. In all cases we assume that the respective DM candidate constitutes all of the galactic DM with a local mass density of ρ DM = 0.3 GeVc −2 cm −3 .
The upper and lower σ pe curves are derived from tracing upper and lower bounds of the published data after applying temperature corrections, along with the nominal curve data that did not have temperature corrections applied. The corrections account for the temperature dependence of indirect, phonon-assisted absorption that occurs at energies below the direct band gap (∼3 eV). We followed the methodology and analytical model for photon absorption found in Ref. [52] to extrapolate the data below 4 eV to a temperature of 50 mK.
This analysis adopted the same ionization production model as used in Ref. [12] to compute the mean number of e − h + pairs n eh produced for an interaction with a given recoil/absorption energy. For recoil/absorption energies above the Si band gap E gap = 1.2 eV but below the average energy per e − h + pair eh = 3.8 eV, n eh = 1; for energies above eh , we determined e − h + pair probabilities from binomial distributions using selected Fano factor values, F .
The total phonon energy measured for an event, E ph , is the recoil/absorption energy E r plus the energy produced by the NTL effect: where the ionization production model and Fano statistics determine the distribution for n eh . We combined the signal models with a charge trapping (CT) and impact ionization (II) likelihood model, which mainly contributes to the distribution of events between quantized e − h + -pair peaks [53]. Charge trapping occurs when an electron or hole falls into a charge vacancy in the crystal, reducing the total number of electrons or holes that traverse the entire detector and lowering the measured energy for an event. Impact ionization occurs when a moving charge in the crystal liberates an additional loosely bound charge, thereby increasing the measured energy for an event.
We determined the CT and II probabilities by fitting the model used in Ref [53] to the laser-calibration data. The results from the fit are 11 ± 3 % and 2 +3 −2 % for the CT and II probabilities, respectively, and were subsequently used to generate the signal models. Because we were unable to determine an energy dependence of the energy resolution within the ROI for this analysis (50-650 eV), the signal models were convolved with a weighted average energy resolution: σ E = 3.6 eV. We determined σ E by averaging over the resolutions of the first six, Gaussian-fitted e − h + -pair peaks from the com-bined laser-calibration data weighted by the corresponding uncertainty in each peak. Lastly, we multiplied each signal model by the efficiency curve (bottom panel of Fig. 1) as well as the exposure (1.2 gram-days). An example of a 1 GeV/c 2 light DM signal model is shown in the top panel of Fig. 1. and FDM ∝ 1/q 2 (bottom) and with Fano factor of 0.155 (solid-blue curve). The light blue band represents our estimate of the systematic uncertainty, which is dominated by varying the Fano factor assumption in the ionization model from F = 10 −4 to 0.3. Other direct detection constraints shown include SuperCDMS HVeV R1 [12] (red), DAMIC [54] (green), SENSEI [23] (orange), XENON10 [55,56] (teal), and XENON1T [57] (pink).

B. Limit Setting
The Poisson exclusion limit for each DM model was calculated independently for the first six e − h + -pair peaks using a limit setting window of ± 3 σ E centered on each peak. While taking into account the look-elsewhere effect, we selected the lowest limit amongst the individual e − h + -pair peaks at each DM mass to obtain a final limit with a 90 % confidence level (C.L.). The light blue band represents our estimate of the systematic uncertainty, which for masses 4×10 −3 keV/c 2 is dominated by varying the Fano factor assumption in the ionization model from F = 10 −4 to 0.3; for masses 4 × 10 −3 keV/c 2 , the uncertainty is dominated by the discrepancy in the photoelectric absorption cross section. Other direct detection constraints shown for V and ALPs include SuperCDMS Soudan [34] (maroon), XENON10 (teal), and XENON100 (purple) [58]; additional constraints on V include SuperCDMS HVeV R1 [12] (red), DAMIC [54] (green), SENSEI [23] (orange), and anomalous energy loss mechanisms in the Sun [24]. For the axioelectric coupling, the entire region shown is disfavored by the observed cooling of red giant [25,26] and white dwarf stars [26,27].
This limit calculation differs from Ref. [12], which determined the limits using the Optimum Interval (OI) method [59,60]. Due to the improved energy resolution of this analysis compared to Ref. [12], the OI method was found to be overly sensitive to the shape of the expected DM signals measured in the detector and thus to the effects of CT and II, leading to systematic uncertainties that are difficult to estimate. In contrast, the Poisson method applied to this analysis is insensitive to these systematic effects. A comparison of the two methods finds up to a factor of 2 stronger limits with the OI method due to the sensitivity to the model shape (the same comparison performed on the Ref. [12] analysis results in no such difference due to the poorer energy resolution). In this analysis we used the more conservative Poisson limit setting method, as it is more effective at constraining the systematic uncertainties.
To quantify the impact of systematic uncertainties, the limits were recalculated with Gaussian distributed random variates for the energy calibration, energy resolution, CT and II fractions, and efficiency, according to their corresponding means and uncertainties. For the photoelectric absorption cross section, we made a random choice between the lower, upper, and nominal curves. At each mass, we took the average from all trials and used the ± 1 σ equivalent values from the resulting limit distribution as the limit uncertainty. The limits and their propagated uncertainty are calculated separately using three different values for the Fano factor: the one measured at high energy, F = 0.155 [61], and the values of F = 10 −4 and F = 0.3 assumed to cover the systematic uncertainty of the Fano factor at these energies. Figures 2 and 3 show the limits on χ scattering and V /ALP absorption, respectively, compared to existing limits. The limits on χ assume a DM form factor of either F DM = 1 or F DM ∝ 1/q 2 [62]. The light blue bands representing our estimates of the systematic uncertainty envelops the ±1 σ values of all three limits obtained using the different Fano factor assumptions in the ionization model. At most masses, the uncertainty bands are dominated by the varying Fano factor assumption; the exception is for 4 eV/c 2 in the V and ALP absorption models, where the uncertainty is dominated by the discrepancy in the photoelectric absorption cross section.

VI. DISCUSSION AND OUTLOOK
The limits in Figs. 2 and 3 are remarkably close to those from our previous run [12] despite the ∼ 2.5 times larger exposure. They are in fact weaker at some masses due to the higher observed event rate in the 3 rd e − h +pair peak coupled with the higher CT and II probabilities in this measurement, as well as the use of a more conservative limit setting method (see Ref. [53] for recent CT and II measurements of the detector used in Ref. [12]). Table I compares the efficiency-corrected event rates for each e − h + -pair peak within a ± 3 σ E window. The event rate observed in each peak is similar in this run compared to Ref. [12] despite a different detector design, cryostat, location, overburden, and shielding.
Another 0.39 gram-days of data were taken with a bias of 60 V across the detector in order to determine if the results are voltage-dependent. Table I shows that the resulting event rate for each number of e − h + pairs is similar to the corresponding rate from the 100 V data, suggesting a voltage-independent result. Furthermore, the event rate above the first e − h + -pair peak is comparable to that seen in other charge-readout experi- in each e − h + -pair peak between this work and Ref. [12]. The event rates displayed from this analysis are calculated from the DM-search data measured with a bias voltage of 100 V, as well as from the additional dataset measured with a bias voltage of 60 V. For each number of e − h + pairs, the event rate is determined by counting the number of observed events within a ± 3 σE window centered on the peak. The uncertainty shown is the 3 σ uncertainty in the number of observed events assuming Poisson statistics.

This Work
Ref. [ ments [12,23,54,63], and adds to the growing narrative of unexplained, O(Hz/kg) low-energy excesses measured in many sub-GeV DM searches (Refs. [64][65][66] and references therein). This result from our detector with unparalleled energy resolution provides a new dataset that can contribute to understanding the origin of these unknown background events. A third run with an identically designed detector is planned in a dilution refrigerator in a shallow underground site with 255 m.w.e. overburden (NEXUS Facility [11]) to probe the correlation between the unknown events and known environmental background sources. Finally, due to the significant impact that charge trapping and impact ionization have on the signal reconstruction, there is an ongoing effort toward understanding these charge propagation effects and investigating factors that influence them. A DM model spectrum with CT and II included is shown in the top panel of Fig 1. The black curve shows the DM signal model for a 1 GeV/c 2 light DM particle with form factor F DM ∝ 1/q 2 , scaled to the limit of excludedσ e , using an II probability of 2 %. The CT probability shown by the hatched region is varied from 0-15 %. For this model and limit-setting scheme, these processes do not determine the ultimate sensitivity. However, lower CT and II rates combined with a more robust understanding [67,68] will allow us to use the region between the peaks in the limit-setting procedure to improve the sensitivity of future analyses, as well as to fully utilize the improved resolution of this detector for additional background discrimination.

VII. ACKNOWLEDGEMENTS
We would like to thank Rouven Essig and Tien-Tien Yu for helpful discussions and assistance with using QEdark [62] to generate the dark matter model used in this analysis. We thank Noemie Bastidon for her work in the preliminary design of our optical fiber setup and wire bonding. We gratefully acknowledge support from the U.S. Department of Energy (DOE) Office of High Energy Physics and from the National Science Foundation (NSF