Molecular Dynamics Simulations on Nuclear Recoils in Silicon Crystals towards Single Electron-Hole Pair Ionization Yields

Abstract

We have developed a novel methodology utilizing molecular dynamics simulations to evaluate the ionization yields of nuclear recoils in crystalline silicon. This approach enables analytical exploration of atomic-scale transport within the lattice without necessitating parameterization. The quenching factors across the nuclear recoil energy range from 20 eV to 10 keV have been thoroughly investigated. A remarkable agreement with experimental data is achieved, particularly for the minimal energy regime conducted to date, reaching the level of a single electron-hole pair. This work presents evidence of a crucial and fundamental distribution of the quenching factor, which can be associated to the collisional interactions underlying the transport phenomena. The region below 4 keV of the quenching factor, where discrepancies have been observed with the Lindhard's model, is found to be significantly attributed to the lattice binding effect and the specific crystal structure. In contrast, a gradual functional relationship is identified below approximately 100 eV, indicating that the quenching factor is influenced by the crystallographic orientation of the target material. From a distributional perspective, our analysis allows for the determination of the minimum exclusion mass for the dark matter nucleon elastic scattering channel at 0.29 $\mathrm{GeV}/c^2$, thereby significantly enhancing sensitivity for the sub-$\mathrm{GeV}/c^2$ mass region.

Figures (5)

• Figure 1: [left] The ionization production scenario involves the primary recoiled atom (red edge), followed by displaced atoms (orange circles) that leave their lattice sites during the cascade, and lattice-bound atoms (blue circles) that remain in position. The grey fog and shadows represent electronic gas and the ionization process, respectively. The dashed circles indicate a perfect lattice without thermal relaxation. [right] Traces of energetic atoms in a recoil event from MD simulation.
• Figure 2: Inherent distributions of QFs. Figures \ref{['fig:qf-dis-lowene']} and \ref{['fig:qf-dis-highene']} exhibit the dispersion properties in the sub-$\text{keV}_\text{nr}$ and $\text{keV}_\text{nr}$ regions, respectively. The expansion of the nonstructured QF distribution (above 200 $\text{eV}_\text{nr}$) is parameterized as $3.64 \times \mathrm{E[keV_{nr}]}^{0.43}$.
• Figure 3: [left] Schematic of Si lattice with crystallographic directions ([100], [001], [111]) and recoil angles $\theta/\varphi$ labeled. [right] Recoil angular dependence of the QF at 20 $\text{eV}_\text{nr}$, showing mean QF over $(\theta,\varphi)$. Attributed to the symmetry of Si, the $\theta$ and $\varphi$ can be reduced to $(0, 90^\circ)$.
• Figure 4: [top] MD-simulated silicon quenching factors (red points; the systematic uncertainties are present with a red shaded area) are compared with experimental measurements (colored crosses) izraelevitchMeasurementIonizationEfficiency2017azecherEnergyDepositionEnergetic1990albakryFirstMeasurementNuclearRecoil2023gerbierMeasurementIonizationSlow1990chavarriaMeasurementIonizationProduced2016doughertyMeasurementsIonizationProduced1992agneseNuclearrecoilEnergyScale2018 with Lindhard-like models superimposed (gray dashed lines, including the original Lindhard model from Eq. \ref{['eq:lindhard-qf']}) lindhard1963integralsarkisIonizationEfficiencyNuclear2023sarkisStudyIonizationEfficiency2020a. [bottom] The relative difference between our calculations and measurements.
• Figure 5: 90% C.L. upper limits for the $\chi$-N interaction derived from SENSEI spectrum adariFirstDirectDetectionResults2025 using our distribution perspective QFs and the Sarkis QF model sarkisIonizationEfficiencyNuclear2023, along with results from silicon-based experiment albakryInvestigatingSourcesLowenergy2022alkhatibLightDarkMatter2021angloherResultsSubGeVDark2023aguilar-arevaloResultsLowMassWeakly2020changFirstLimitsLight2025.