Examining real-time time-dependent density functional theory nonequilibrium simulations for the calculation of electronic stopping power

In ion irradiation processes, electronic stopping power describes the energy transfer rate from the irradiating ion to the target material's electrons. Due to the scarcity and significant uncertainties in experimental electronic stopping power data for materials beyond simple solids, there has been growing interest in the use of first-principles theory for calculating electronic stopping power. In recent years, advances in high-performance computing have opened the door to fully first-principles nonequilibrium simulations based on real-time time-dependent density functional theory (RT-TDDFT). While it has been demonstrated that the RT-TDDFT approach is capable of predicting electronic stopping power for a wide range of condensed matter systems, there has yet to be an exhaustive examination of the physical and numerical approximations involved and their effects on the calculated stopping power. We discuss the results of such a study for crystalline silicon with protons as irradiating ions. We examine the influences of key approximations in RT-TDDFT nonequilibrium simulations on the calculated electronic stopping power, including approximations related to basis sets, finite size effects, exchange-correlation approximation, pseudopotentials, and more. Finally, we propose a simple and efficient correction scheme to account for the contribution from core-electron excitations to the stopping power, as it was found to be significant for large proton velocities.

Figure 1

Electronic stopping power as a function of the velocity of a proton projectile in silicon. The solid black curve represents the empirical PSTAR model [45]. The sets of letters correspond to different experimental data sets that can be found in Ref. [46]. The green curves represent results from the dielectric response formalism using LR-TDDFT for calculating the dielectric matrix [25] with the core-electron correction (dark green) and without (light green). The solid blue curve corresponds to the values we obtained by calculating using RT-TDDFT with the “standard” parameters and approximations described in the Method section. The error bars represent standard deviations of the mean stopping power along the given proton trajectory, which is calculated as an average of the instantaneous force on the projectile proton.

Figure 2

Electronic stopping power as a function of the velocity of a proton projectile in silicon. The black curve represents the empirically fitted PSTAR model [45]. Convergence of the electronic stopping power curves with respect to various parameters is shown: (Left) Calculations with 80 Ry plane-wave kinetic cutoff energy and 50 Ry plane-wave kinetic cutoff energy. (Center) Calculations with four $k$ points in the Brillouin zone and only the gamma point in the Brillouin zone. (Right) Calculations with a 512-atom supercell and a 216-atom supercell.

Figure 3

Electronic stopping power as a function of the velocity of a proton projectile in Si. The solid black curve represents the empirical PSTAR model [45]. The colored curves correspond to the values we obtained by calculations using mixed Gaussian plane wave (GPW) RT-TDDFT using the PBE XC functional (dashed blue) and the PBE0 hybrid XC functional (red). The error bars represent standard deviations of the mean stopping power along the given proton trajectory, which is calculated as an average of the instantaneous force on the projectile proton.

Figure 4

Electronic stopping power as a function of the velocity of a proton projectile in Si. The solid black curve represents the empirical PSTAR model [45]. The colored curves correspond to the values we obtained by calculating using mixed Gaussian plane wave (GPW) RT-TDDFT in which core electrons are treated using pseudopotentials (blue) and treated explicitly using the all-electron approach (green). The error bars represent standard deviations of the mean stopping power along the given proton trajectory, which is calculated as an average of the instantaneous force on the projectile proton.

Figure 5

Population analysis based on projections of the time-dependent Kohn-Sham (TDKS) states onto the ground eigenstates of the silicon system. Percentages of population loss relative to the initial, fully occupied condition for the $3s$ and $3p$ atomic orbitals (red), the $2p$ orbitals (light blue), the $2s$ atomic orbitals (blue), and the $1s$ atomic orbitals (black) of the silicon atoms.

Figure 6

Total energy as a function of projectile proton ($v=6.0\mathrm{a}.\mathrm{u}.$) displacement in the simulation cell acquired from the RT-TDDFT simulations using the mixed GPW/GAPW approaches. In each plot, the results from three different treatments for the core ($1s, 2s, 2p$) electrons are shown: The black curve represents the simulation results in which the core electrons are represented by pseudopotentials (PP). The blue curve represents the all-electron (AE) results from using the GAPW method in which the core electrons are treated explicitly and allowed to be excited. The red curve represents the frozen-core (FC) results in which the AE wave functions are used but the core electrons are not propagated in time. (Left) The results from the simulations for the centroid path trajectory. (Right) The results from the simulations for a path parallel to the centroid path trajectory but 33% closer to the channel of target Si atoms.

Figure 7

Scheme for the single-atom collision simulations, where the red arrow indicates the trajectory of the proton (${\mathrm{H}}^{+}$). The simulation is carried out over a range of impact parameters h. The distance between the proton and the silicon target atom is given by r. The displacement of the proton is given by $l$.

Figure 8

Energy transferred in single-atom collision simulations as a function of impact parameter $h$ for the simulations using pseudopotentials (black), all electrons (blue), and frozen core (red) schemes. See the text for details.

Figure 9

Electronic stopping power as a function of the velocity of a proton projectile in Si. The solid black curve represents the empirical PSTAR model [45]. The sets of letters correspond to different experimental data sets that can be found in Ref. [46]. The dark blue curve corresponds to the values we obtained by calculating using mixed Gaussian plane wave (GPW) RT-TDDFT with the “standard” parameters and approximations described in the Method section. The light blue curve corresponds to the results after the core correction has been applied. The error bars represent standard deviations of the mean stopping power for the trajectory.

Figure 10

Electronic stopping power as a function of the velocity of a proton projectile in Si. The solid black curve represents the empirical PSTAR model [45]. The sets of letters correspond to different experimental data sets that can be found in Ref. [46]. The dark red curve corresponds to the values we obtained by calculating using plane-wave pseudopotential (PW-PP) RT-TDDFT with the “standard” parameters and approximations described in the Method section. The light red curve corresponds to the results after the core correction has been applied. The error bars represent standard deviations of the mean stopping power for the trajectory.