Electronic stopping power in insulators from first principles

Using time-dependent density-functional theory we calculate from first principles the rate of energy transfer from a moving proton or antiproton to the electrons of an insulating material, LiF. The behavior of the electronic stopping power versus projectile velocity displays an effective threshold velocity of ~0.2 a.u. for the proton, consistent with recent experimental observations, and also for the antiproton. The calculated proton/antiproton stopping-power ratio is ~2.4 at velocities slightly above the threshold (v~0.4 a.u.), as compared to the experimental value of 2.1. The projectile energy loss mechanism is observed to be stationary and extremely local.

The interaction of charged particles with matter has been a subject of extensive research since the discovery of subatomic particles [1]. Ions moving through solids gradually transmit their kinetic energy to electronic excitations of the host and deposit it along their path. The maximum of this deposited energy is the so-called Bragg peak and it occurs shortly before the particle stops. Hence the importance of studying the electronic energy loss of slow ions (with velocities below the Bohr velocity) travelling through solids. For metals, the understanding of this problem has been steadily progressing over the years [1]. For insulators, however, experimental results remain unexplained even for simple systems. This is particularly true at low velocities (the threshold effect) [2,3], where the contribution from nuclear collisions conceals the electronic stopping [4].
A detailed quantitative knowledge of these processes is required to understand the damage produced in materials when exposed to radiation. For ceramic materials devised for the encapsulation of nuclear waste [5] the prediction of durability over extremely long times is crucial. Radiation-damage simulations performed to date [6,7] rely on empirical force fields obtained from fits to lowenergy properties. The actual interatomic forces could be enormously altered, however, by the local electron heating produced by the electronic stopping.
In the semiclassical formalism, the electronic energy loss rate is given by the response of the system to the external potential, dE dt = −Zv · E ind , where Z and v are the charge and velocity of the projectile, and E ind is the induced electric field in the target material. Hence, to first order, linear response theory can be used to give the stopping power (SP) -the energy lost per unit path length-in terms of the dynamical dielectric function of the material. For metals this formalism shows that the SP is linear with small velocities, dE dx ∼ (Ze) 2 v [8,9], reflecting that no minimum energy is required to excite electron-hole pairs. It must be remembered, however, that the dielectric approach is not valid at low ion velocities, and non-linear effects cannot be neglected [1].
A different behavior is expected in wide band-gap insulators, given the finite energy required to excite electrons. For protons moving through noble gases [10], which also display a large minimum excitation energy, a velocity threshold has been observed experimentally [11], and explained in terms of the quantization of energy levels and charge exchanges. However, for solid insulators the situation is unclear. No threshold effects were originally observed in Al 2 O 3 , SiO 2 , or LiF, and the linear dependence dE/dx ∝ v was observed down to velocities of about 0.3 a.u. [2,12,13].
For protons under grazing incidence in LiF(001), and below ∼ 0.2 a.u. a threshold behavior was reported [14]. Under these conditions the proton does not penetrate the solid, and charge exchange is identified as the dominant mechanism for electronic stopping, with local electron capture from F − ions, giving rise to H 0 and H − . More recent experiments on thin LiF films show an apparent velocity threshold near 0.1 a.u. [3]. The different experimental setup (transmission) suggests a different physical stopping mechanism, based on electron-hole pair excitations. A threshold behavior is expected, and has been qualitatively predicted from linear-response and from perturbation theory calculations. These approximations, however, grossly underestimate the SP at low velocities, thus exaggerating the threshold [15]. A theoretical description beyond these approximations is needed.
In this letter we present a first-principles approach to the non-perturbative study of realistic solid-ion interactions. We follow the explicit time evolution of the electronic states of the host crystal as an external particle propagates through the system, by means of time- dependent density-functional theory (TD-DFT) [16]. We use here this scheme to understand the SP threshold effect and stopping mechanisms in LiF, a well-studied characteristic insulator, finding reasonable agreements with measured data. Most importantly, however, this study sets the scene for a promising line of theoretical simulations, assesses its possibilities and offers new insights into the stopping process.
The energy transmitted to the electrons from a constant-velocity moving ion is monitored. Total energy is thus not conserved since the energy loss of the projectile will be neglected (its large mass ensures a negligibly small decline in its velocity on the time scales of the simulations). As the center of the charged particle, r c , moves with constant velocity, r c = r 0 +v·t, the time-dependent Kohn-Sham (KS) equation [16] defines the dynamics of effective single-particle states (and thereby the electronic density and energy) under the external potential generated by the projectile and the crystal of Li and F nuclei.
The calculations were done using the Siesta ab initio method [17,18], in its time-evolving TD-DFT implementation [19], using the instantaneous local density approximation (LDA) to exchange and correlation [20]. The 1s core electron pair of F was replaced by a norm-conserving pseudopotential [21] in the fully nonlocal form. The 1s electrons of Li were explicitly included in the calculation as pseudized valence electrons. A double-ζ polarized basis was considered for the valence electrons (single-ζ for the 1s of Li), with an energy shift of 100 meV [18]. The grid cutoff [18] for integration was 118 Ry. The lattice parameter obtained for bulk LiF is 3.98Å, slightly smaller than the experimental value 4.03Å, as expected for LDA. The projectile (proton, p, or antiproton,p) was described by the bare ± 1 r potential [32]. It was not pseudized in order to treat p andp on the same footing.
Periodic boundary conditions were used throughout. The supercell size was chosen so as to minimize the spurious effects of the repetition while keeping manage-able computational demands. After convergence tests, a 4x4x4 supercell of 128 atoms was selected, such that a particle moves ∼11Å down a [110] channel before reentering the box. The Γ-point approximation was used for integrations in k. The projectile was initially put in the centre of a crystal cage and the time-independent DFT solution was obtained to define the initial state for the subsequent evolution. It was then moved with constant velocity along a [110] channel. The time-dependent KS equation was then solved numerically by discretizing time and applying the Crank-Nicholson algorithm as described in ref. [19]. Using a time step of 1 attosecond, the wavefunctions were then propagated for several femtoseconds. The electronic total energy was recorded as a function of time and the SP, dE dx , was extracted. from the average slope of this stationary regime. Stopping power results for p andp are presented in Fig. 2. A threshold effect at low velocities is apparent in the figure, in agreement with recent experiments [3], unlike the linear behavior observed for metals [1], but much smaller than that predicted for insulators by linear response theory [15]. The values obtained in this work (around 0.2 a.u.) are consistent with experiments [3,13]. The difference between p andp is also apparent, in contrast with the invariance under charge-sign change expected from perturbative treatments. Experiments show that the energy loss for protons is approximately twice as high as for antiprotons [13], consequence of the different kind of screening (one attracts electrons while the other attracts holes). In the neighborhood of v = 0.4 a.u. (slightly above the threshold) we obtain a SP ratio of ∼ 2.4, that compares well with the experimental value of 2.1 [13].
Since a q 2 charge dependence is expected from linear response, the SP over q 2 is plotted in Fig. 3 for both projectiles. The smooth behavior at the origin indicates that there are no substantial biases in the way the stopping mechanism is described for p andp. The slope near the origin corresponds to the q 3 dependence, while the bending corresponds to q 4 and higher terms. The displayed behavior, and particularly the sign of both terms, is similar to that observed for metals [22].
There is a discrepancy between measured and computed values for the SP for velocities above the threshold. The experiments show four times higher SP [3]. A factor of two is accounted for by the known relation between channelling conditions and the average over random trajectories [23,24]. Technical reasons account for a part of the remaining discrepancy. The inclusion of basis orbitals along the projectile's path (a sp single-ζ set every 0.5Å) increases the SP from 1.6 eV/Å to 2.8 eV/Å for v = 0.46 a.u.
[33] as shown in Fig. 2. At a more fundamental level, discrepancies are expected from the errors in electronic spectra around the band gap given by instantaneous LDA in TD-DFT (even if it includes a RPA correction on KS eigenvalues) [25]. We now analyze the locality of the energy absorption in terms of Mulliken charges [26]. Although they are known to be unreliable in their explicit dependence on the basis set, their changes for fixed basis set are much more indicative. By comparing the atomic charges obtained dynamically and adiabatically (a static self-consistent calculation with the external potential at the same instantaneous position), insights are obtained into the stopping process. Fig. 4 shows this comparison for the proton case, where the following is observed: (i) The screening of the proton is enhanced in the dynamic case (our basis set describes this screening by a transfer of charge from F to Li), due to the increased polarizability of the host at the frequencies explored by the moving ion. (ii) A delay is apparent in the dynamical screening, which is the main responsible of the dissipation. (iii) The electronic excitation process is extremely local: the changes observed in atoms closest to the trajectory are much larger than any other, and the dynamical screening effect is only noticeable when the projectile is very close to the ion.
The locality of the energy transfer is confirmed when calculating the energy absorbed by a small cluster of LiF (Li 6 F 5+ ). The inset of Fig. 2 shows a striking similarity in the overall v-dependence of the energy absorbed by the cluster and the SP in the solid. If we took the effective path length in the cluster as 1.4Å (the number of valence electrons in the cluster corresponds to one formula unit), both SPs would be indeed very similar [34]. In the solid, the energy accumulated along the path would then diffuse away at longer time scales (corresponding to the effective band width), defining a wake. A tentative definition of local energy [35] shows well differentiated time scales for the excitation by the projectile, on one hand, and the ensuing out-diffusion, on the other (not shown). The short-ranged initial excitation can be rationalized in terms of the electronic localization length scale relevant to dielectric response [27], which is expected to be very short for LiF. Similar locality in the energy loss mechanism has been recently found in a confined twodimensional electron gas, however [28].
A characterization of the charge state of the projectile [29,30] has not been attempted. The present calculations allow for the establishment of charge states of any kind within the constraints that the electronic charge and momentum are conserved. It is clear that the mid-gap state travels with the proton, and that it becomes partly populated. It is not clear, however, that the charge associated to that state represents a meaningful definition of the charge state of the projectile, since whether the screening charge builds up as it passes or travels with it is not determined by the population of that state. This effect will be explored in further works.
In conclusion, we have presented a general approach for the non-perturbative first-principles study of the electronic stopping power in solids. New insights into the electronic SP in insulators have been provided for pro- tons and antiprotons in LiF. We thank P. M. Echenique, E. Salje and I. Nagy for interesting discussions. EA acknowledges the hospitality of the Donostia International Physics Centre. This work has been funded by the EC's Marie Curie program, UK's BNFL and NERC, the Spanish MEC, the UPV/EHU and the Basque Government through the Etortek program. r instead of 1 r and converging the SP with respect to ρ. This is equivalent to a compensating homogeneous background. The finite-cutoff filtering of the 1 r cusp customary in plane-wave calculations [31] has been replaced here by a Gaussian charge smoothening of a width dictated by the grid cutoff used in siesta [18].
[33] The basis set used throughout this work has been the one defined by atomic orbitals of target atoms, so that no explicit bias was introduced in the comparison of p andp. There is a finite-basis saturation effect at high velocities, but that is observed at higher velocities than the ones studied here. Although the absolute value of the SP is not completely converged by our basis set, the convergence for both the relation proton/antiproton and the respective thresholds is satisfactory.
[34] TD-DFT is known to provide better results for small systems [25], we thus expect the spectra of our cluster to be more accurate than that of bulk.
[35] The energy associated to a basis function |µ is defined as Eµ = ν ρµνHνµ, where ρµ,ν and Hµν are the density and Hamiltonian matrices, respectively. Note that it does not include double-counting terms.