surface backgrounds on the CCDs. The measured bulk ionization spectrum confirms the presence of an excess population of low-energy events in the CCD target with characteristic rate of ∼7 events per kg-day and electron-equivalent energies of ∼80 eV, whose origin remains unknown. DOI: 10.1103/PhysRevD.109.062007

The DArk Matter in CCDs (DAMIC) experiment searches for the interaction of dark matter particles in the galactic halo [1] with silicon atoms in the fully depleted active region of scientific charge-coupled devices (CCDs). Between 2017 and 2019, DAMIC acquired 11 kg-days of data with an array of seven CCDs with conventional readout (pixel noise σ pix ∼ 1.6 e -) installed in a low-background setup deep underground at SNOLAB [2]. In December 2020, DAMIC reported a statistically significant (3.7σ) excess of events above its background model between the threshold energy of 50 eV ee and 200 eV ee [2]. The population of the excess events was best described by ionization events uniformly distributed in the bulk of the CCDs, with an exponentially decaying energy spectrum with decay energy ε ¼ 67 AE 37 eV ee and a rate of 5.1 AE 2.3 events per kg-day. The spectral fit was performed to all ionization events above the energy threshold in depth vs energy space, with a full background model including different templates for bulk and surface events. This approach was necessary because the bulk and surface populations could only be distinguished statistically, since the determination of the depth of the interaction was limited by readout noise. Consequently, the dominant systematic uncertainty was identified as the modeling of surface backgrounds, and the statistically significant rise was reported only as an unmodeled excess of events over background.

In November 2021, the DAMIC setup was upgraded with two skipper CCDs, which can achieve single-electron resolution. This article presents results from a 3.25 kg-day target exposure acquired throughout 2022 with the upgraded detector [3]. The order-of-magnitude improvement in readout noise, from σ pix ∼ 1.6 e -to 0.16 e -, results in a significant improvement in the depth localization of ionization events. This allows for a new analysis based on a fiducial selection to obtain a clean sample of bulk events down to a threshold of 23 eV ee . In addition, the lower threshold is expected to increase the measured rate of reconstructed excess events from 1.7 to 3.0 events per kg-day.

The DAMIC skipper upgrade features a few notable differences from its predecessor [2]. The previous seven DAMIC CCDs were decommissioned to install two 24-megapixel DAMIC-M skipper devices [4] in a new oxygen-free high conductivity (OFHC) copper box. The stock copper was stored in SNOLAB for five years to suppress cosmogenic activation [5], and brought to the surface at the time of machining. The total sea-level exposure of the machined components is less than nine days. The box with the two CCDs is connected to the existing cold finger to bring the sensors to a temperature ≲140 K. The two ancient lead bricks employed in the previous installation are positioned above and below the CCD box to provide additional shielding from environmental backgrounds. CCD 1 is located at the bottom of the box, and CCD 2 at the top, with their front surfaces (on which the pixel array is patterned) oriented upward. There is a 2.5 mm vertical gap between CCDs, with no material in between. The previous vacuum interface board (VIB) was replaced with a new one designed for skipper devices, which remains shielded by 18 cm of lead located above the CCDs [2]. CCD control and readout is performed with the Low Threshold Acquisition (LTA) electronics designed for SENSEI [6]. In the setup, two LTA boards are synchronized on the same clock signal allowing low-noise readout of both CCDs. Three-meter-long coaxial cables carry the signals from the LTAs through the 42-cm polyethylene neutron shield to the air side of the VIB.

The 24-megapixel skipper CCDs were designed by Lawrence Berkeley National Laboratory [7] and fabricated by Teledyne/DALSA. The CCD substrate consists of 670 μm of fully depleted, high-resistivity (>10 kΩ • cm) n-type silicon. Ionizing particles produce electron-hole pairs in the silicon as they lose energy. The charge carriers are drifted across the substrate, in the vertical direction (ẑ), by means of a bias voltage applied to a backside contact, V sub ¼ 60 V. Carriers experience thermal motion as they drift in the substrate electric field, which leads to a spatial variance in the transverse plane (xŷ) that is proportional to the transit time. Holes are collected at a potential minimum created by the buried p-channel below the pixel array. The pixel array of each CCD consists of 6144 × 4128 (columns × rows) pixels, each of size 15 μm × 15 μm, for a total sensitive mass of 8.9 g per CCD. Clock signals are sent to the three-phase polysilicon gates on each pixel to transfer the collected charge across CCD rows (ŷ direction) and into the horizontal register (bottom row), where it is transferred along x to two skipper amplifiers located at opposite ends (U and L) of the horizontal register. Skipper amplifiers can perform multiple nondestructive charge measurements (NDCMs) of the charge in a single pixel. Taking the average of N sample pixel samples improves the readout noise by a factor 1= ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi N sample p , which enables single charge resolution with enough samples. The energy depositions from low-energy ionizing particles-whose tracks are much shorter than the pixel size-result in pixelated 2D Gaussian distributions of charge in the images. The lateral spread of the charge distribution, σ xy , is positively correlated with the depth of interaction, with a maximum value σ max ∼ 1 pix.

Sources of single-electron events relevant in skipper CCDs have been characterized in Ref. [8]. DAMIC-M data at low charge multiplicities are consistent with a Poisson background from leakage current in the CCDs [4]. Uncorrelated single-electron events arise from a combination of instrumental and radiation backgrounds, and are hereafter referred to as dark counts. To minimize leakage current, CCD potentials were inverted to fill surface traps in several steps during cooldown, starting at 160 K until the final cold finger temperature of 110 K. CCD biases were optimized to suppress light emission by the amplifiers, and the clock values were selected to minimize spurious charge [8] while maintaining high charge transfer efficiency. Following commissioning and cooldown, the pixel readout noise for a single measurement (N sample ¼ 1) was ∼4 e -in all four skipper amplifiers. The background dark count rate in the images was between 2.4 × 10 -3 and 3.0 × 10 -3 e -=pix=day. The rate of accidental noise clusters, which sets the analysis threshold, depends on a combination of pixel readout noise, shot noise from dark counts, and the exposure per pixel. While increasing the signal integration time and the number of NDCMs decreases the readout noise, it also increases the overall readout time and, thus, the exposure per pixel. An increase in the exposure per pixel increases the dark counts per pixel, and hence the shot noise. We performed studies on simulated images with the measured noise profile of the CCDs to establish the readout parameters. For science data, CCDs were read out continuously with a signal integration time t int ¼ 8 μs and N sample ¼ 460. 1 This resulted in a pixel exposure time of ∼50 hours and a pixel noise σ pix ∼ 0.16 e -. The output data were stored as arrays of size 3300 × 210 × 460 (N columns × N rows × N sample ), with one file per amplifier. The CCD physical array is divided by two in the x direction since two amplifiers are used for readout. Each amplifier reads 3072 physical pixels per row and an additional 228 pixels past the physical extent of the array, which constitute the overscan. Overscan pixels have significantly shorter exposure than physical ones, and allow a check on noise and charge transfer efficiency. Unlike previous DAMIC data runs, the charge was read out for every physical pixel, rather than after summing the charge from a group of pixels, since this provides better spatial resolution.

A total exposure of 4.81 kg-day was acquired between February 2022 and January 2023 in seven data runs. For the first run, in which a 1 kg-day exposure was collected, images with a smaller overscan (N columns ¼ 3100) were acquired. Images were processed as follows. Images with pixel values in analog-to-digital units (ADUs) were constructed by averaging the value for each pixel across NDCMs. The first NDCM was excluded to avoid the noise transient at the start of pixel readout. This procedure produces a 3300 × 210 (3100 × 210) average image. The median pixel value is a good estimator of the baseline since >96% of pixels have an occupancy of 0 e -. Thus, for every row in the average images, the median value was subtracted from every pixel value to correct for any baseline shift between rows. The distribution of averaged pixel values (PVD) exhibits discrete peaks corresponding to the pixel charge content (0 e -, 1 e -, 2 e -, etc.) that is dominated by dark counts. The PVD of every image was fit to a Poisson distribution convolved with a Gaussian function to evaluate the gain (from the mean of the Gaussian, in ADU=e -), noise (Gaussian standard deviation) and background dark counts (the Poisson mean). A total of 10124 images were inspected during daily shifts and 142 were excluded based on several quality criteria, i.e., if the gain was outside the 3σ range defined for each amplifier and run, the pixel noise did not decrease with the expected 1= ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi N sample p , or there was a prominent gradient in the baseline across the image. Most of the discarded images were found to be correlated with electrical power instabilities at SNOLAB. Finally, 40 average images were concatenated into a joint 3300 × 8400 (3100 × 8400) image to reduce the number of ionization events split between images.

A mask was defined per amplifier to exclude defects and regions with increased dark counts. Defects in the silicon lattice alter the local band gap, which leads to increased charge leakage at specific locations in the CCD [9]. Because of continuous readout, all defects appear as "hot" columns in the science run images. Any columns with λ i > 35 × 10 -3 e -=pix were excluded from the analysis since their PVD is typically not well described by a Poisson distribution, as estimated from a chi-square goodness-of-fit test. Few columns with smaller λ i but exhibiting a non-Poisson charge distribution (p-value < 3.2 × 10 -5 ) were also excluded. Columns with x > 1565 in the 1L amplifier (L amplifier of CCD 1) were discarded due to inefficient charge transfer past this point, possibly due to a charge trap in the horizontal register. To maximize the sensitivity to defects, a set of images with N sample ¼ 1 was acquired at a cold finger temperature of 160 K, where charge leakage from defects is orders-of-magnitude larger. These "warm" images were read out after a 30-minute exposure of the pixel array, with defects appearing as localized hot pixels due to the relatively short readout time. We then generated a "median image," where each pixel contains the median value of the pixel over all warm images, and identified defects as contiguous pixels > 5σ above the background noise. In addition, we evaluated the median pixel value for every column over all images in the science data, and estimated a moving average of the column medians in a 200-column window. Any column in the science data with median value 1.3 times (2σ) higher than the moving average and coincident with a defect in the warm images was excluded. Finally, to reject events that are split between adjacent images, we masked a border of 10 pixels in the x direction and five pixels in the ŷ direction around the joint image. The amplifier masks account for a 52%, 28%, 12% and 36% loss in exposure in the 1L, 1U, 2L and 2U amplifiers, respectively.

To remove backgrounds associated with high-energy ionization events, pixels were also masked on a per-image basis. Pixel clusters were identified as groups of contiguous pixels with charge ≥3 e -. The cluster energy was estimated as the sum of pixel values in the cluster assuming 3.8 eV ee =e - [10]. For every cluster with energy >10 keV ee , we masked the smallest rectangular region of pixels containing the cluster plus a 2-pixel wide border around it. Trails of charge from inefficient transfer of large charge packets from high-energy depositions were removed by excluding the 800 (100) trailing pixels in the horizontal direction in the L (U) images, and 20 trailing pixels in the vertical direction. These selections were obtained by studying the PVD in the trailing regions until there was no evidence of excess charge above background dark counts. We observed Oð10 -4 Þ cross talk between amplifiers, with high-energy events read out by one amplifier resulting in few-electron events in the other amplifiers. While this is most relevant for α decays, we conservatively masked all pixels that were read out at the same time as pixels within the rectangular box surrounding high-energy clusters. In addition, pixels for which the standard deviation of the NDCMs was larger than 10 e -were excluded. This selection accounts for occasional inaccurate estimates of the pixel charge due to noise fluctuations. The per-image masks result in only a 0.6% loss in exposure with small variations per amplifier. Overall, image selection and masking results in a final exposure of 3.25 kg-day.

The cluster search was performed using a likelihood clustering algorithm on the images [11]. The likelihood Lðq ij jN; ⃗ μ; σ xy ; λ i ; σ pix Þ defines the probability that the pixel values q ij in a specified window are described by a 2D Gaussian distribution of charge with mean position ⃗ μ, lateral spread σ xy and amplitude N on top of shot noise with mean occupancy λ i , where every pixel has white noise σ pix . The images were scanned with a moving window of 5 × 5 pixels, where, at every window position, we computed the likelihood that pixels come only from noise (L n ∶ N ¼ 0) or noise plus a Gaussian charge distribution (L g ∶ N 0 ; ⃗ μ 0 ; σ xy 0 ). The initial guesses for N 0 and ⃗ μ 0 are the total charge in the window and the center of the window, respectively. The initial guess for the lateral spread of the cluster σ xy 0 was set to 1.0 pix. The values for λ i and σ pix were fixed to the measured values for the specific image. At window positions where there is a clear preference for a Gaussian cluster (ln L gln L n > 7), we vary the window position until ln L gln L n is maximized. We then fix the window and perform a log-likelihood optimization with N, ⃗ μ and σ xy as free parameters. The test statistic ΔLL ¼ -ðln Lgln L n Þ, with Lg the best-fit likelihood, is a metric for the preference of the Gaussian hypothesis, with more negative values corresponding to a higher probability of ionization event within the selected window. The best-fit values for N, ⃗ μ and σ xy are taken as the best estimates for the energy, mean ðx; yÞ position, and lateral spread of the cluster. A second fit was performed with σ x and σ y as independent free parameters.

A selection on ΔLL was used to reject accidental clusters from noise. We started by simulating "blank" images, which contain only the measured pixel readout noise and shot noise from dark counts as a function of column number in the data. We ran the likelihood clustering on 30 kg-day of simulated blanks to get a distribution of accidental noise clusters, which had a maximum charge of 10 e -. We compared the distributions between data and blanks for clusters with q ≤ 10 e -and -25 < ΔLL < -10, where the spectrum is dominated by noise clusters. Excellent agreement in both cluster charge and the x coordinate of the cluster center was found. To obtain an accidental rate of R noise < 0.01 ðkg-dayÞ -1 uniformly across the charge range ½5; 10 e -, we select events with ΔLL ≲ -29, with small (at most AE2) variations between charge bins because of changes in the shape and amplitude of the ΔLL distributions. Above 10 e -, a fixed selection ΔLL < -28 was used.

A series of criteria was applied to select "valid" clusters. We select clusters whose fit window does not contain any masked pixels and which are properly centered (⃗ μ < 1.2 from window center). To select events with the expected 2D Gaussian topology, we first excluded extended clusters with σ x or σ y greater than 1.5 σ max . Ionization events interacting in the field-free regions past the horizontal register can produce horizontal clusters with degraded energy [12]. Such asymmetric clusters were rejected by excluding events with σ y < 0.2 and σ x > 0.5 pix. We applied an analogous rejection for vertical events (σ x < 0.2 and σ y > 0.5), which can originate from charge released by traps in the CCD pixel array. Finally, to exclude clusters from trails of charge that were missed by the mask, we require that no more than 8% of pixels within AE50 pixels about the center of the cluster along the row had a value ≥1 e -.

The variance of the lateral diffusion of charges (the lateral spread of a low-energy event) can be modeled as [2,7]:

where z and E denote event depth and energy, respectively, and parameters A and b depend on the microscopic properties of the substrate at operating temperature and V sub . Best-fit values for the parameters A ¼ 3.07 pix 2 and b ¼ 5.35 × 10 -4 μm -1 were obtained in a surface setup using straight cosmic muon tracks [11]. The values of α ¼ 0.889 and β ¼ 7.4 × 10 -3 keV -1 ee were obtained from a fit to backside events (predominantly 210 Pb and 210 Bi decays) as a function of energy in the SNOLAB data. We performed a selection on σ xy as a function of energy to remove 95.4% of ionization events originating from the surfaces of the CCDs, which was determined by simulating events < 10 μm from the front and back surfaces and adding them onto blank images. This selection was validated over the full energy range in the surface laboratory by illuminating with a 14 C β source (Q β ¼ 156 keV) the front and the back sides of a CCD from the same batch as those installed at SNOLAB. The CCD was operated with the same parameters, readout settings and dark counts as in SNOLAB. The leakage fraction of events following σ xy fiducial selection was 4 AE 1% and 3 AE 1% for backside and frontside illumination, respectively.

The selection efficiency for bulk events was obtained by simulating ionization events uniformly distributed in z on blank images, applying the masks from data and running the likelihood clustering. Our diffusion model was previously validated for bulk events by irradiating a CCD with low-energy (< 24 keV) neutrons and γ rays [13,14]. Figure 1 shows the efficiency of our selections as a function of energy. After ΔLL and valid-cluster selections, we obtain ≥ 10% efficiency for bulk events down to 6 e - (23 eV ee ), which we set as the analysis threshold. The efficiency plateaus at 95% from the selection of valid clusters, with 4% of clusters lost because the fit window contains at least one masked pixel. The AE1σ uncertainty band is the statistical uncertainty in the number of selected events in the simulated blanks from the ΔLL selection. For the σ xy selection, we estimated the 68% confidence interval from simulations where the diffusion model parameters were varied within their uncertainties.

Although the copper box and CCD packages of the upgraded detector have a different geometry than in the previous installation, the apparatus components (i.e., CCDs, Kapton flex cables, glue, OFHC copper, brass screws and ancient lead), their composition and radioactivities remain nominally the same. To confirm that the ionization spectrum features the same dominant background components as before, we first analyzed the spectrum above 0.5 keV ee , where no excess above the background model was previously observed. The overall (bulk) background rate in the 1-6 keV ee energy range of 9.7 AE 0.8 ð4.4 AE 0.6Þ keV -1 ee kg -1 day -1 is comparable to the previous installation [2]. A spectral analysis up to 20 keV ee showed the characteristic spectrum of cosmogenic tritium (Q β ¼ 18.6 keV) in the bulk of the CCDs over an approximately constant background from Compton scattering of external γ rays. Surface events and x-ray lines from 210 Pb decays were also identified with comparable rates as before.

After establishing the final cluster selections and performing background studies, we unblinded the region of interest (ROI) below 0.5 keV ee . The clusters retained after the ΔLL and valid-cluster selections are shown in Fig. 2, with energy and σ xy projections beside the corresponding axes. The valid criteria remove 11 clusters that are contained in the fit window (with no masked pixels) and fall in the ROI. The σ xy projection shows the distribution of all ionization events from threshold to 1.0 keV ee , while the energy projection shows the spectrum of fiducial events after σ xy selections, demarcated by the black dashed lines in the scatter plot. Clusters in the ROI were individually inspected and were found to have the expected topology, located away from masked regions, and not spatially correlated with other ionization events. The event ðx; yÞ and time distributions were confirmed to be statistically consistent with uniformity, with Kolmogorov-Smirnov p-values between 0.36 and 0.78. The σ xy distribution was found to be consistent from a fit (p-value ¼ 0.73) with three distinct populations of bulk, front-and back-surface events. A total of 15 events pass the fiducial σ xy selection in the ROI, in contrast with the 4.8 AE 0.7 expected from the background rate measured in the 1-6 keV ee range. The increasing rate of events toward low energies at large σ xy is caused by partial charge collection of ionization events in the backside of the CCD, as described in Ref. [2].

Our background model [2] predicts an approximately constant spectrum of bulk events at low energies. This is expected since electronic recoils can only be induced in the bulk by β decays or Compton scatterings, which both have FIG. 1. Selection efficiency for simulated bulk events after noise (blue line) and surface-event (red line) rejection with their respective uncertainty bands. The gray region below 5 e -was not considered in this analysis. approximately flat spectra in the ROI. 2 Nuclear recoils from fast neutrons and thermal neutron captures are expected to have orders-of-magnitude lower rates [2]. In addition, based on the number of surface (bulk) events rejected (retained) by the σ xy selection and the corresponding uncertainty from the diffusion model, we expect 1.8 þ1.9 -1.0 leakage events from the surface in the ROI. We performed an extended unbinned likelihood fit to the energy spectrum of fiducial events with a flat background component and an exponentially decaying spectrum, which adequately parametrized the observed excess in the previous analysis. The leakage is constrained to the expected value within uncertainty and assumed to contribute to the flat and exponential components as to have the same spectrum as the fiducial events in the ROI. The free parameters in the fit are the integrated counts in the background (b), leakage (l) and excess (s) spectra and the decay energy (ε) of the exponential. The spectra were corrected for the bulk event acceptance (red line in Fig. 1) and the fit was performed between 0 and 6 keV ee . The results are shown in Fig. 3 for E < 1 keV ee . The best fit finds s ¼ 11.0 þ4. 4 -3.8 excess events with ε ¼ 89 þ36 -24 eV ee , and b ¼ 3.57 þ0.52 -0.34 background events in the ROI. A likelihood-ratio test between the null hypothesis (s ¼ 0, with b and l free) and the best fit results in a p-value of 7.73 × 10 -4 , corresponding to a significance of 3.4σ. We translate the number of fiducial excess events to a total rate of 10.0 þ4.0 -3.4 per kg-day from the fiducial exposure. The parameter space of the excess compared to the result from the previous 11 kg-day FIG. 2. Distribution of clusters in energy vs σ xy space, for E ∈ ½0.023; 1 keV ee . Clusters with large (small) σ xy correspond to events occurring in the back (front) of the CCDs. Crosshatched gray regions are excluded by the 95% surface event rejection (dashed black), which results in the final sample of fiducial events (filled circles). The green-shaded region represents the energy ROI (E < 0.5 keV ee ) for the excess. The left (bottom) panel shows the σ xy (energy) projection of all (fiducial) events. The depth (z) scale from the nominal diffusion model is reported on the right-hand side axis. FIG. 3. Result of the spectral fit to fiducial events with all parameters free (red) and the null hypothesis (blue dashed). The fit prefers an excess low-energy exponential component. Counts in the legend are reported for the ROI. The top axis shows the nuclear-recoil energy scale from Ref. [16].

exposure [2] is shown in Fig. 4. The two results are statistically compatible and suggest a common origin of the excess population.

In summary, an excess population of bulk events above the expected background spectrum has been observed in two setups of DAMIC at SNOLAB, which employed CCDs with significantly different noise characteristics and readout modes. Both CCDs with conventional [14] and skipper [15] readout have demonstrated the capability to correctly resolve spectral features at the energies where the excess is observed. The topology of the events-which can be resolved with high resolution thanks to the subelectron noise of skipper CCDs-is consistent with ionization events. The modeling of surface backgrounds, which was the dominant systematic uncertainty in the previous analysis, has been addressed by analyzing a clean sample of bulk events selected using lateral spread of the clusters (σ xy ).

The observed excess ionization events likely arise from an unidentified constant source of radiation in the DAMIC detector or from the environment, which is common to the two experiments. As such, this excess is distinct from the excess of phonon signals reported in milli-Kelvin cryogenic calorimeters [17], which are likely caused by stress released by the crystal [18]. The only known interactions that could give rise to the observed excess spectrum are those from neutrons with silicon nuclei in the bulk of the CCDs. The observed spectrum could be reproduced by the scattering of neutrons with energies up to ∼17 keV and a flux of ∼0.2 cm -2 d -1 through the CCDs, but no such source of neutrons has been identified. Turning to more exotic interpretations, the bulk excess spectrum is well described by nuclear recoils from interactions of weakly interacting massive particles (WIMPs). For spin-independent WIMPnucleus coherent elastic scattering with standard galactic halo parameters [19], the excess corresponds to a WIMP with mass ∼2.5 GeV=c 2 and a WIMP-nucleon scattering cross section ∼3 × 10 -40 cm 2 . This interpretation is nominally excluded by results from CDMSlite [20] and DarkSide-50 [21]. Attempts to find a consistent interpretation between these experiments by systematically varying detector response (e.g., nuclear-recoil ionization efficiencies), WIMP speed distribution, as well as alternate particle interaction models (e.g., Ref. [22]) are beyond the scope of this article.

The SENSEI experiment [23], currently taking data at SNOLAB, is designed to acquire a larger exposure of ∼40 kg-day with skipper CCDs. The DAMIC-M detector [24], a 0.7 kg skipper-CCD array with an improved radioactive background rate of Oð0.1Þ keV -1 ee kg -1 day -1 , is scheduled to start operations at the Modane Underground Laboratory in 2025. If the bulk excess is detected in SENSEI and DAMIC-M, the significantly increased statistics will enable a high-resolution spectral measurement, studies of the time evolution of the excess, and investigations of its dependence on detector configuration and operating parameters to better understand its origin.

We are grateful to SNOLAB and its staff for support through underground space, logistical and technical services.