Energy loss due to defect formation from ²⁰⁶Pb recoils in SuperCDMS germanium detectors

Crystal defects can occur when incident radiation recoils off of an atom transferring sufficient energy to displace the atom from its lattice site, thus creating a vacancy. If the displaced atom remains in the crystal, it is referred to as an interstitial atom (or an “interstitial”). The combination of the vacancy and the interstitial is referred to as a Frenkel pair or Frenkel defect.¹ Energy can also be lost through creation of defect clusters, dislocations, and amorphous zones. The creation of these defects permanently stores energy in the crystal, with the fraction of incident energy that goes into defect formation depending in part on the mass of the impinging particle, the deposited energy, and the crystal properties.² Collectively, the total energy lost to formation of defects is referred to as the Wigner energy.³ 

The energy required to displace an atom from its lattice site is the displacement threshold energy. For germanium, previously determined displacement threshold energy values from theory and various molecular dynamics simulations are inconsistent, ranging from 7 to 30 eV.^4–10 

The value of the displacement threshold energy has implications for physics experiments that employ solid-state detectors to search for nuclear recoil events with sub-keV energy depositions. In this letter, we focus on defect formation with data from the Super Cryogenic Dark Matter Search (SuperCDMS) experiment,^11–17 which aims to detect nuclear recoils from weakly interacting massive particles (WIMPs)¹⁸ by measuring the energy deposited when a WIMP scatters off of an atomic nucleus in a detector's crystal lattice. The SuperCDMS program is targeting low-mass WIMPs¹⁷—from a few hundred MeV/c² to several GeV/c²—using advanced detector designs having detection thresholds on the order of the Ge-atom displacement energy. Because the energy that goes into the formation of defects is not directly observable, an accurate determination of the energy loss to defect formation is important for understanding the low-energy detector response and thus for discerning the ultimate low-mass-WIMP sensitivity reach.

To measure the Wigner energy associated with ²⁰⁶Pb-on-Ge interactions, we consider data from the most recent phase of the SuperCDMS experiment,^11–16 when it was located in the Soudan Underground Laboratory. ²¹⁰Pb sources were deployed adjacent to two detectors to evaluate their in situ response to non-penetrating radiation from the decays of ²¹⁰Pb and its daughters ²¹⁰Bi and ²¹⁰Po (see Ref. 11 and Fig. 2 therein). These data include ²⁰⁶Pb-on-Ge recoils, for which a significant disagreement between the simulated and measured spectra is evident near the expected 103 keV endpoint energy (cf. Fig. 4 in Ref. 11). In this letter, we reconsider this discrepancy while allowing for the possibility that the measured recoil energy is effectively reduced due to formation of defects. The measured and simulated ²⁰⁶Pb spectra are compared using a χ2 statistic to find a best-fit energy-loss fraction that brings the two into agreement. The results from each detector are calibrated for events near the detector surface (vs. in the bulk) to obtain a value for energy loss due to defect formation in Ge.

The SuperCDMS Soudan experiment operated 15 cylindrical, interleaved Z-sensitive Ionization and Phonon (iZIP) Ge detectors at ∼mK from 2012 to 2015,^11,12,16 arranged in five stacks of three detectors each. Data from the top and bottom detectors of the third such stack—called T3Z1 and T3Z3, respectively—are used in this study.

Each iZIP detector had several independent phonon and ionization readout channels on both of its flat faces. The ionization electrodes on the top and bottom faces were biased at +2 V and −2 V, respectively, while the interleaved phonon sensors were held at ground. The resulting electric field caused positive and negative charge carriers from particle interactions in the detector bulk to drift to opposing faces, whereas within ∼1 mm of either face, most of the charge carriers were collected by the electrodes on that face of the detector. This asymmetry in charge collection between the two detector faces makes it possible to distinguish energy depositions near a detector face from those in its bulk (cf. Fig. 3 in Ref. 11).

The SuperCDMS Soudan experiment had two in situ ²¹⁰Pb sources, with one installed facing the top side of T3Z1 and the other facing the bottom side of T3Z3. The sources were produced by exposing silicon wafers to a 5 kBq ²²⁶Ra source (which produces ²²²Rn gas) for 12 days inside a sealed aluminum box. After exposure, the wafers were surface etched to remove dust and radon daughters resting on the surface. This process resulted in a near-uniform implantation profile of ²¹⁰Pb to a depth of approximately 58 nm.^11,19 Based on the subsequent time of exposure to lab air, we estimate a 1.6 ± 0.1 nm oxide layer formed on the surface of each source wafer.^20,21

The data used in this analysis were collected from March 2012 to July 2014.^11,12 Ionization and phonon signals were measured for each event, and the ratio of these measurements (“ionization yield”) allowed for discrimination between event types. The detector responses were calibrated using ¹³³Ba gamma rays such that electron recoils in the detector bulk have ionization yield equal to one. ²⁰⁶Pb recoils have comparatively low ionization yield; in Fig. 1, they appear at a yield of ∼0.3 and they extend in energy to near the expected 103 keV endpoint.

In this study, ²⁰⁶Pb recoils are selected based on their ionization yield and the surface-event criteria developed in Ref. 11. Similar criteria are used to select near-surface electron recoil events (highlighted in Fig. 1) that correspond to gamma rays (top box) and betas (middle box) from decays of ²¹⁰Pb and ²¹⁰Bi. These event selections are used to estimate the detector resolution and energy scale for surface events, independent of the ²⁰⁶Pb recoils used to study energy loss from defect formation.

The SuperCDMS Soudan experiment was simulated with Geant4^22–24 version 10.1.p2 using the Screened Nuclear Recoil physics list.²⁵ A detailed simulation geometry was used including the detectors, all surrounding materials, and the ²¹⁰Pb source wafers. The source wafers were simulated with zero surface roughness. The full chain of ²¹⁰Pb decays was simulated according to the source wafers' implantation profile. One million primary ²¹⁰Pb atoms were simulated, with Geant4 allowed to handle the full decay chain for each event.

Selecting simulated ²⁰⁶Pb events in the 80–110 keV region of interest yields a total of ∼44 000 simulated events, approximately twice the corresponding number of measured ²⁰⁶Pb recoils. The simulation results show good agreement with the shape of the measured spectra up to 80 keV. However, as shown in Fig. 2, for larger recoil energies, there is a significant discrepancy between the measured and simulated spectra for each detector. The former are softer with substantially fewer events measuring the full 103 keV endpoint energy. This disagreement near the ²⁰⁶Pb-recoil endpoint is indicative of energy loss due to defect formation in the detector crystal, a process not taken into account in the simulation.

Two factors are considered to account for the discrepancy in Fig. 2: energy loss due to defect formation and energy smearing due to detector resolution. Each simulated event is first scaled to

where $E$ is the deposited energy, $f$ is the fraction of energy lost to defect formation, and $Ê$ is the remaining energy. $Ê$ is then treated with a Gaussian smearing function that has a standard deviation corresponding to the 1$σ$ detector resolution at that energy. The resulting smeared event energy $Ẽ$ is thus representative of the actual energy measured by a detector.

As demonstrated in Ref. 19, the resolution is an approximately linear function of energy in the range of 80–110 keV: $σE=0.63 keV+0.024E$⁠. The parameters are estimated by fitting to the 46 keV and 66.7 keV peaks in surface-event gamma-ray spectra. We assume that the energy resolution of ²⁰⁶Pb recoils has the same functional form, but the absolute value may differ slightly. To account for this difference, the resolution is scaled by a multiplicative factor $Ps$

where $σPb$ is the resolution function used to smear the simulated ²⁰⁶Pb recoil energies. Both $f$ and $Ps$ represent free parameters that are allowed to float in the fitting method outlined below.

After the simulated events are scaled and smeared, the resulting energies are compared directly to the measured energies as follows. Let $A$ and $B$ represent the set of measured and simulated event energies, respectively,

where the sets are of size $N$ and $M$⁠, respectively. Each set is binned by energy into $q$ bins

where $ai$ and $bi$ indicate the number of events in the ith bin. To gauge the level of agreement between these binned energy distributions, a χ2 statistic is calculated

We generate approximately one million sets $B$ for each detector, corresponding to different combinations of the scaling (⁠$f$⁠) and smearing (⁠$Ps$⁠) parameters. A χ2 value is determined for each set, creating a well-defined parameter space from which a minimum can be found yielding the best-fit values for $f$ and $Ps$⁠.

The measured event energies are based on the detectors' default energy calibrations, which are developed using gamma rays in the bulk of the crystal. The energy scale for surface events may be slightly different than for bulk events.¹¹ Consequently, the measured ²⁰⁶Pb recoil energies may differ from their simulated counterparts by an additional energy scaling factor that represents an intrinsic miscalibration and therefore is independent of defect formation. If present, a best-fit determination of the scaling factor $f$ in Eq. (1) would account for both this miscalibration and energy loss due to defect formation

where $fDF$ is the scale factor from energy loss due to defect formation, and $fsur$ is the surface-event scale factor. Because the total energy loss to defect formation depends on the mass of the incident particle, surface events from gamma rays and betas should have $fDF∼0$ to within the precision of this study. This allows for the determination of any intrinsic miscalibration via an independent examination of these alternate event classes. The energy loss to defect formation is thus

with $f$ determined from ²⁰⁶Pb events and $fsur$ determined from surface gamma-ray and beta events.

Application of the procedure outlined above to the measured and simulated ²⁰⁶Pb recoil energies gives a best-fit energy-scale parameter of $f=5.52±0.10%$ and $6.67±0.11%$ for detectors T3Z1 and T3Z3, respectively. Figure 3 shows the χ2 statistic as a 2-dimensional function of the smearing strength $Ps$ and the energy loss parameter $f$⁠. Statistical uncertainties (at 1$σ$ confidence) on the best-fit values of $f$ are determined by projecting the $Δχ2=1$ contours onto the “Energy Scale Parameter” axis. After application of the best-fit parameters, the simulated and measured ²⁰⁶Pb recoil energy distributions are in good agreement, as shown in Fig. 4.

The same analysis procedure is applied to simulated and measured distributions of surface-event gamma rays and betas highlighted in Fig. 1, with the results summarized in Table I. Because $fDF∼0$ for these event classes, the values in Table I are a direct measure of $fsur$⁠. A single value of $fsur$ is obtained for each detector by taking a weighted mean of the gamma-ray and beta results, which is then used to determine $fDF$ from Eq. (4); these results are summarized in Table II. The weighted mean of the two detectors is $6.08±0.08%$⁠, but because the individual measurements differ by 2.06 standard deviations (p-value 0.04), the uncertainty on the weighted mean is increased by a factor of 2.06 by increasing each uncertainty by that factor. This results in a more reasonable difference of one standard deviation and gives a weighted mean of $6.08±0.17%$ where the uncertainty is a combination of statistical and systematic uncertainty.

There is an additional systematic uncertainty of $%−0.03+0.04$ from the $1.6±0.1$ nm silicon dioxide layer on the surface of each source wafer. The best-fit energy loss is therefore $6.08±0.18%$ after adding the uncertainties in quadrature. Additional potential sources of uncertainty are discussed in the below paragraph.

Using our best estimate of the value for energy loss to the formation of defects, it is possible to determine the displacement threshold energy of a germanium atom. This is the average displacement threshold energy over all lattice angles²⁶ and is an important quantity for radiation detectors, WIMP-searches, and other applications.^10,27,28

For interactions involving the same species of incident and target atoms, the Kinchin-Pease equation estimates the number of defects formed²⁹ (with further refinement by Norgett et al.³⁰) In the case of an incident ²⁰⁶Pb recoil on Ge, displaced Ge atoms may be liberated with enough energy to form yet more defects; so there are two types of interactions to consider. TRIM-2013³¹ simulations were used to model the entire defect formation process, for a range of user-defined values of the Ge displacement threshold energy from 15 to 23 eV.

The target material in the TRIM simulations was a solid mass of pure Ge with a thin layer of GeO₂ on top. As with Si, pure Ge reacts with oxygen in the atmosphere to create GeO₂ with a thickness that logarithmically depends on exposure time.^32,33 We estimated a GeO₂ layer thickness of $0.98±0.02$ nm.

TRIM predicts a monotonic, decreasing relationship between the percent energy lost to defects in the 80–110 keV energy range and the Ge displacement threshold energy. To match our best estimate of the energy loss value of $6.08±0.18%$⁠, TRIM simulations suggest using a displacement threshold energy of $19.7−0.5+0.6$ eV. The systematic error does not include modeling imperfections in Geant4 and TRIM.

This value is somewhat in tension with some molecular dynamics calculations.^10,34,35 However, TRIM uses simple potentials and includes tuned parameters to fit experimental implantation data. The more sophisticated potentials used in molecular dynamics simulations may yield different values. More experimental data are required to further investigate this.

Other sources of uncertainty were considered, resulting in no significant increase in the quoted uncertainty.

We investigated the effects of varying the thickness of the germanium oxide layer on top of the detectors. We estimated the thickness to be $0.98±0.02$ nm, and varying the thickness by 1$σ$ did not change our results at the precision given.

In this analysis, we make the assumption that all recoils are ²⁰⁶Pb events. However, there are some events where sputtered silicon atoms from the source wafers might contribute to the total event energy. Considering the incident energies and formation of defects for both Si and Pb ions shifts our best-fit result by less than one percent of the value obtained with the Pb-only assumption.

Nevertheless, the defect energy loss parameters for the two detectors are not quite statistically consistent (cf. f_DF in Table II). This inconsistency may represent a true physical difference due to differences in crystal properties between the two detectors. It may also be a result of an operational difference. The ²⁰⁶Pb recoils used in this analysis were incident on opposite faces of the two detectors (i.e., top vs. bottom), which were biased with opposite polarities and thus resulted in collection of predominantly positive or negative charge carriers by the ionization electrodes. A corresponding difference in charge collection efficiency (electrons vs. electron-holes) for ²⁰⁶Pb recoils relative to surface-event gamma rays and betas may explain the apparent inconsistency between the two detectors.

The ability of SuperCDMS iZIP detectors to differentiate event types was leveraged to find the Wigner energy following ²⁰⁶Pb implantation on Ge. We used this result with TRIM simulations to determine an average displacement threshold energy of $19.7−0.5+0.6$ eV for germanium. This value will play a critical role in understanding the sensitivity of future experiments designed to measure nuclear recoils (e.g., from dark matter interactions) in Ge detectors, especially as instruments move toward lower energy thresholds and better resolution. Our results also provide another important, empirically determined value from which the Stillinger-Weber potential³⁶ or others could be fit. Future detectors with thresholds on the order of the displacement threshold energy could confirm this result more directly with low-energy neutron calibrations.

The SuperCDMS collaboration gratefully acknowledges technical assistance from the staff of the Soudan Underground Laboratory and the Minnesota Department of Natural Resources. The iZIP detectors were fabricated in the Stanford Nanofabrication Facility, which is a member of the National Nanofabrication Infrastructure Network, sponsored and supported by the NSF. We would like to thank François Schiettekatte for his feedback in reviewing an early manuscript. Funding and support were received from the National Science Foundation, the U.S. Department of Energy, Fermilab URA Visiting Scholar Grant No. 15-S-33, NSERC Canada, the Canada First Research Excellence Fund, and MultiDark (Spanish MINECO). This document was prepared by the SuperCDMS collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359. Pacific Northwest National Laboratory is operated by Battelle Memorial Institute under Contract No. DE-AC05-76RL01830 for the U.S. Department of Energy. SLAC is operated under Contract No. DEAC02-76SF00515 with the U.S. Department of Energy.