Electronic stopping power in a narrow band gap semiconductor from first principles

The direction and impact parameter dependence of electronic stopping power, along with its velocity threshold behavior, is investigated in a prototypical small band gap semiconductor. We calculate the electronic stopping power of H in Ge, a semiconductor with relatively low packing density, using time-evolving time-dependent density-functional theory. The calculations are carried out in channeling conditions with different impact parameters and in different crystal directions, for projectile velocities ranging from 0.05 to 0.6 atomic units. The satisfactory comparison with available experiments supports the results and conclusions beyond experimental reach. The calculated electronic stopping power is found to be different in different crystal directions; however, strong impact parameter dependence is observed only in one of these directions. The distinct velocity threshold observed in experiments is well reproduced, and its non-trivial relation with the band gap follows a perturbation theory argument surprisingly well. This simple model is also successful in explaining why different density functionals give the same threshold even with substantially different band gaps.


I. INTRODUCTION
The study of fast-moving charged particles shooting through solid materials started with Rutherford's famous experiment of showering a gold foil with alpha particles to substantiate the nuclear model of the atom 1 . Such fast-moving particles strongly perturb the target material. The perturbed state of the medium either relaxes to its original state or to a new state with structural defects, depending on the nature of the interaction. The study of such defects, generally referred to as 'radiation damage', is of great interest from the point of view of applications ranging from nuclear engineering 2 to biological soft matter for medical applications 3 and materials engineering for space electronics 4,5 .
Stopping power is a quantitative measure of the interaction between the projectile and the target medium, defined as the energy transferred from the former to the latter per unit distance traveled through the material. The fast-moving charged particle dissipates its kinetic energy by collisions with the nuclei and the electrons of the medium. Therefore, it is traditional to differentiate between these two distinct dissipation channels; the loss of energy to electronic excitations is known as the electronic stopping power S e , and the loss of energy to the nuclear motion is known as the nuclear stopping power S n .
There is a growing interest in modeling the stopping power of ions with velocities between 0.1-1 atomic unit (a.u. hereafter) 6 . In this regime the electronic stopping power (ESP) is generally dominant; however, at lower velocities the contribution from nuclear collisions also becomes sizable 7 . Fermi and Teller 8 , using an electron gas model, found the ESP to be proportional to the projectile velocity for v < 1 a.u. Lindhard 9 and Ritchie 10 , applying a linear response formalism to an electron gas model of simple metals, predicted a linear velocity dependence within the low projectile velocity limit. Almbladh et al. 11 showed, by calculating the static screening of a proton in an electron gas using density functional theory (DFT), the significant limitations of the linear response treatment. Using DFT Echnique et al. 12,13 proposed a full non-linear treatment to account for nonlinear effects such as the presence of bound states and the complex electronic structure of the heavy projectiles in the low velocity limit. Recently, the modeling of proton and antiproton stoppings in metals, using jellium clusters as a model of the target, has been extended to intermediate and high projectile velocities using realtime time-dependent density functional (TD-DFT) 14,15 simulations 16,17 .
Fermi and Teller 8 pointed out that, in case of insulators, the linear velocity dependence of the ESP is only valid in the limit in which the kinetic energy of the projectile is greater than the band gap. An extensive amount of interesting work has been carried out on the problem of ESP within the linear response theory 18-21 and nonlinear formalism 22 . A detailed background on the subject can be found in Ref. 23 and references therein. A vast majority of these approaches is limited to an electron gas model of metals and do not take into account important features such as the local inhomogeneity of the electron density, core state excitations and band gaps in case of insulators and semiconductors. These features become increasingly important at low velocities. Radiation damage in metals has also been studied, obtaining very interesting qualitative results describing the processes in model systems, using explicit electron dynamics within a tight binding model 24,25 .
Relatively recently, TD-DFT based first principles calculations of ESP [26][27][28][29] have been performed for insulators and noble metals to explain some interesting effects observed experimentally 30-33 which do not fit the known theoretical models 12,23 . These TD-DFT based calculations have successfully reproduced the expected threshold behavior in wide band gap insulators, and the role of d electrons in the non-linear behavior found in gold. In contrast, there has not been much work done on semiconductors, except for a study 34 which investigated oscillations in the ESP by varying the atomic number Z. However, no systematic velocity-dependent investigation has been attempted at this level of theory. Recent experiments show a possible small velocity threshold for protons in bulk Ge, a system with very small band gap 35 . The band gap of Ge is almost 20 times smaller than that of LiF while the observed threshold velocity in Ge is only 2 to 3 times smaller. Very little is known about the velocity threshold in small band gap materials.
Experimentally it is almost impossible to measure directly the ESP at velocities 0.2 a.u., as usually the total stopping power S = S n + S e of the medium is measured. The ESP can then be extracted from the measured spectrum using different models 36,37 . However, a quantitative knowledge of all possible mechanisms contributing to the total stopping power is necessary to extract the electronic component properly. At velocities not much higher than 0.1 a.u. it becomes rather difficult to disentangle the two contributions 38 . However, in simulations it is possible to directly access the ESP using TD-DFT based non-adiabatic electron dynamics simulations. In such simulations the projectile is directed along a crystal direction, where it does not get too close to any of the target nuclei. The nuclear contribution to the stopping power, therefore, is negligibly small and can even be completely suppresed by constraining the host atoms to be immobile and the projectile velocity to be constant.
In this study we have investigated the ESP of H in Ge. A small band gap and relatively low packing density makes Ge particularly interesting for the investigation of the threshold behavor which has been observed in wide band gap insulators 29,33 . The simulations have been carried out using an equivalent method to Refs. 28 and 29. Furthermore, we have systematically studied the direction and impact parameter dependence of the ESP, for which very little is known. The accuracy offered by this method, as verified in the satisfactory comparison to experiments below, allows us to explore these aspects explicitly.

II. METHOD
The calculations are carried out using an extension of the Siesta program and method 39,40 which incorporates time-evolving TD-DFT based electron dynamics 41 . The ground state of the system is calculated with the projectile placed at its initial position. The ground state Kohn-Sham (KS) orbitals 42 serve as initial states. Once the ground state of the system is known, the projectile is given an initial velocity and the KS orbitals are propagated according to the time-dependent KS equation 14 using the Crank-Nicholson method with a time step of 1 asec. The forces on the nuclei are muted and the projectile velocity is kept constant so that energy is transferred only through inelastic scattering to the electrons. In any event the projectile velocities are fast enough to leave little or no time for the nuclei to respond. The total energy of the electronic subsystem is recorded as a function of the projectile displacement. The gradient of the total energy against projectile displacement, after filtering out periodic oscillations due to the crystal periodicity, gives the ESP S e . The Kohn-Sham orbitals were expanded in a basis of numerical atomic orbitals of finite extent 43,44 . A doubleζ polarized (DZP) basis set was used to represent the valence electrons of the projectile and the host material, while the core electrons were replaced by norm conserving Troullier-Martins pseudopotentials 45 , factorized in the separable Kleinman-Bylander (KB) form 46 . The effect of the Ge pseudopotential was checked by introducing the core (3d ) electrons into the valence electrons. Details of the basis set and the pseudopotentials are given in Appendix A. The sampling of the real-space grid, for representing the electronic density and basis functions for the calculation of some terms of the Hamiltonian matrix 40 , was chosen to correspond to an energy cutoff of 200 Ry. We used a 96-atom supercell constructed by 2×2×3 conventional cubic cells of Ge as shown in Figure  1. A k -point mesh of 4 × 4 × 3 points generated with the Monkhorst-Pack method 47 corresponding to an effective cutoff length of 22.36Å 48 was used after testing its convergence. The exchange and correlation functional was evaluated using the local density approximation (LDA) in the Ceperley-Alder form 49 .
We used the theoretical lattice constant, which was found to be 5.59Å, compared to an experimental value of 5.66Å. This underestimation of ∼ 1% is typical for the LDA. An indirect band gap of 0.70 eV was found for bulk Ge, compared with an experimental value of 0.74 eV (at 0 K). However, it is important to note that this good agreement is fortuitous, as DFT with LDA generally either underestimates the band gap or does not produce one at all. Pseudopotential can be one the sources of cancellation of errors 50 along with a smaller lattice parameter which tend to open the band gap. Lee et al. 51 , using a plane-wave method, have reported an indirect band gap of 0.41 eV. Much larger band gap, up to 0.81 eV 52 , have been reported depending upon the details of the calculation. The dependence on the density functional was checked by repeating the calculations for the Perdew-Burke-Ernzerhof (PBE) functional 53 , for which the theoretical lattice constant was found to be 5.78Å with a direct band gap of 0.33 eV.
In order to check the convergence of our basis in Siesta, we have also computed the band structure with the plane-wave DFT code Abinit 54 , making use of exactly the same pseudopotential including the same choice of local potential and KB projectors, and a high kinetic energy cutoff of 95 Ry for the basis. The agreement for the valence and low-lying conduction bands is excellent, although we find a slightly smaller band gap of 0.58 eV with the plane-wave calculation (see Appendix B).

III. RESULTS AND DISCUSSION
In an experiment with a polycrystalline sample the projectile gets channeled along different crystal directions. We have therefore taken into account the direction and impact parameter dependence. We have computed the ESP along three different channels. The calculated ESP is compared with experimentally measured values by Roth et al. 35 in Figure 2.  Figure 9).

A. The velocity threshold
The ESP varies linearly with projectile velocity, intercepting zero at a finite velocity. This indicates a definitive threshold, within an accuracy for the ESP determination of ∼ 0.01 eV/Å. Roth et al. 35  The threshold behavior has been observed in insulators both experimentally and theoretically. From perturbation theory a relationship between the projectile velocity and electronic transitions is given by (see, e.g., where v is the projectile velocity, ∆k is the change in momentum in electronic excitations, and ε g is the band gap and we are taking = 1 for simplification through out this article. This relation can be deduced by requiring the conservation of energy and momentum for a twoparticle collision event in the limit of mass of projectile M → ∞ (see Appendix C). Following equation 1, the velocity threshold for an indirect band gap modelled as in Figure 3 would correspond to the relation (see Appendix C) where m e and m h are the electron and hole masses, respectively, k 0 is the difference in crystal momentum between the valence band maximum and the conduction band minimum, ε g is the indirect band gap. It follows that for small k 0 the threshold returns to the direct band gap behaviour (see Ref. 55), and v th ∝ √ ε g . In the case when both parabolas are thin on the scale of k 0 , i.e., when k 0 ≫ (m e + m h )ε g , the threshold velocity rather goes as v th = εg k0 and is thus linear with ε g . This argument implies that, for parabolic bands, below a threshold velocity the ESP would drop to zero. For the case of periodic bands, however, this threshold would not be strict, but can still be defined within some accuracy depending on the smoothness of the projectile's potential convoluted with the relevant electronic wave functions 55 . From Equation 1, a threshold velocity in a given direction can be estimated from the band structure of the material by finding the gradient of the line which is a joint tangent to the valence and conduction bands, shown by the arrow in Figure 3. The threshold velocity estimated from the band structure in the [001] direction is found to be 0.053 a.u. as shown in Figure 4 Figure 4. However, the ESP calculated using LDA and PBE does not differ significantly at low velocities, and the two calculations produce almost the same threshold. This is a clear indication that the threshold phenomenon is not straightforwardly related to the band gap. The gradient of the joint tangent line of the valence and conduction bands in both cases is almost the same (shown by the solid arrows in  Figure 4). This suggests that the behavior of the ESP threshold at low velocities is rather related to the indirect band gap in the given direction regardless of its being the absolute gap. This further supports the above described model of the ESP threshold. The fact that the relation in Equation 2 is accurate using the unperturbed host band structure is somewhat surprising. Such agreement is due to the fact that the perturbing projectile potential does not significantly affect the band structure around the gap.

B. Direction and impact parameter dependence
We have found that the ESP strongly depends on direction in the crystal, particularly at high velocities. The difference in the ESP between the [111] and [001] channels is up to 3%, and between these two and the [011] channel it is up to 33%. The electron density along these channels is shown in Figure 6 in suitable planes. The electron density is then averaged over the z-axis, as shown in Figure 7. The direction with the lowest ESP for a channeled projectile ([011]) has a lower average density in the center of the channel compared with the two other channels. For [001] and [111] the averaged density is not significantly different, similarly what happens for the ESP. This suggests that the ESP in channeling conditions can be related to the average density along the trajectory, corroborating and supporting assumptions and approximations used in the literature 56,57 .
We have simulated five different trajectories in the [001] channel, as shown in the inset of Figure 8. The five trajectories are chosen to sample different impact parameters (different closest distance to any of the host atoms) within the channel. For each trajectory we show the total energy of the electronic subsystem versus distance for a given velocity of 0.5 a.u. in Figures 8 and 9. The plots in Figure 8 show the energy profile along the [001] channel; the periodic variation in the electronic energy reflects the periodicity of the crystal. A larger variation is seen for the trajectories with the lowest impact parameters, as should be expected; however, the base-lines of all the trajectories have the same gradient, which shows that, in this direction, the ESP is quite insensitive to impact parameters. A similar calculation in the [111] direction gives the same result (not shown). However, the ESP strongly depends on impact parameter in the [011] direction. The total electronic energy profile for five different trajectories in this direction is shown in Figure 9. The change in ESP from the highest impact parameter, i.e., the center of channel (empty circle data points in Figure 2) and the lowest impact parameter, i.e., close to the edge of channel (empty square data pionts in Figure 2) changes by a factor of two. Again looking at the average density in the [011] direction ( Figure 7) we can see that it changes by a factor of three from the center to the edge of the channel. This reflects the proposed strong correlation between the ESP and the averaged local density within a small radius of the impact parameter. It is to be expected that such radius (or cross section) would increase for slower projectiles. This is verified by the larger slope of the ESP for the center of the [011] channel trajectory for lower velocities. Indeed, the low velocity limit displays the same behavior for all trajectories, indicating that the larger cross section is seeing the same average electron density in all the cases. In experiment the ESP is naturally averaged over different directions and impact parameters, and precise knowledge of this averaging mechanism would be necessary to obtain a comparable average from our calculations. We have not attempted to do so, although it is clear from Figure 2 that any such averaging would result in a slight underestimation with respect to experiment, especially for high velocities.

IV. SUMMARY
We have systematically studied the different aspects of the ESP of H in bulk Ge, a representative narrow band gap semiconductor for which good experimental results are available. We have learned that the ESP is sensitive to the crystal direction and, in certain directions, to the choice of impact parameter. A detailed model is needed to average the calculated ESP over different directions.  The band structure and density of states of bulk Ge calculated using Siesta (LCAO) and Abinit (Plane Waves) is compared in Figure 10. The same pseudopotential (and its local and non-local components) is used in both codes. Appendix C

Threshold Velocity
This is a known relationship that can be obtained in several different ways; here, we present one such way of deriving it. If a particle of mass m and initial momentum k i collides with another particle of mass M and initial momentum K i , conservation of momentum requires that where k f and K f are the final momenta of the particles, respectively, and ∆k denotes the change in momentum. Conservation of energy requires that where ε i and ε f are initial and final energies of the particle of mass m, respectively. From equation C1, we can write On substituting equation C3 in equation C2, we obtain In the limit M → ∞, the second term in equation C4 vanishes, and the rest simplifies to where v = Ki M . The smallest excitation in the system would require ε f − ε i = ε g , where ε g is the band gap of the material, with an accompanying change in momentum ∆k of the electron undergoing the transition. The threshold velocity of the projectile at the onset of energy loss would therefore relate to the band gap as: (C6)

Indirect band gap
The argument for deducing the excitation condition in a direct band gap case can be extended to the case of parabolic bands with an indirect band gap. The condition for the direct band gap (ε g = 1 2 (m e + m h )v th ) can be found in Ref. 55. A geometrical way to proceed for the indirect band gap is to find the conditions for which a straight line (corresponding to the red arrow in Figure  3) would cross both of the parabolas, and from these derive the limiting velocity value below which there is no crossing. Considering first the parabola for electrons, we can write ε e = 1 2m e |k e − k 0 | 2 + ε g .
The transition line ε t = k e · v + ε 0 should cross the parabola ε e , where ε 0 is a constant defining the vertical positioning of the transition line of slope v (red arrow in Figure 3): 1 2m e |k e − k 0 | 2 + ε g = k e · v + ε 0 .
Here for simplicity we consider that k 0 and v are collinear. Furthermore, since we are interested in obtaining an equation for the threshold velocity, we can consider that k e is parallel to v without loss of generality. The equation C8 is quadratic in k e and can be solved to give k e = k 0 + m e v ± (k 0 + m e v) 2 − 2m e (ε g − ε 0 ) − k 2 0 .
Similarly, for holes we can write Again, the transition line ε t = k h ·v+ε 0 should cross this parabola. Equating the two gives a quadratic equation in k h which can be solved to give The two conditions C9 and C11 (for electrons and holes, respectively) can be combined as for that to be possible, leading to or, at v = v th ,