# Polarization induced by charged particles in real solids

## I Introduction

A charged particle penetrating a solid causes a distortion in the electronic density around the particle and behind its position: this is what Bohr called the induced density wake. Related to this induced density there is an induced potential. For a sufficiently high velocity the induced wake shows an oscillatory behaviour. When the velocity of the projectile is larger than the average velocity of target electrons (typically $`v_F`$, the Fermi velocity), one may consider a linear response of the medium. However, in the case of projectiles moving with smaller velocities, nonlinearities may play an important role for metallic densities ($`2<r_s<6`$, $`r_s`$ being the average electron distance).

Much research has been oriented to the study of these quantities. Pioneering work on dynamic screening was performed by Neufeld and Ritchie. They evaluated the induced potential and density wake by using a local dielectric function as the linear response function of the medium. Based upon a linear response of the target, different aspects related to the induced polarization such as wake riding states or the spatial distribution of the induced potential and density have been investigated. The induced potential and density have been studied beyond linear-response theory, in the static electron gas approximation and within hydrodynamical formulations. The quadratic induced polarization by an external charge in the full random-phase approximation (RPA) has been recently reported.

In previous works the induced potential and density have been evaluated on the basis of a jellium model of the target, which consists of an isotropic homogeneous electron gas embedded in a uniformly distributed positive background. However, in a more realistic approach valence electrons move in a periodic potential and one-electron excitations split into the so-called energy bands. The impact of band-structure effects on plasmon dispersion curves, dynamical structure factors, stopping power, and hot-electron lifetimes has been investigated only very recently. It has been shown that these effects can be important even in the case of free-electron metals such as aluminum. The wake potential induced by swift protons through different solids has been studied using a model dielectric function.

In this paper we report, within linear-response theory, an ab initio evaluation of the induced potential and density of ions moving through real solids. Within our description we can obtain the induced potential under channeling conditions, i. e., when the projectile penetrates through different symmetry directions of the solid. In section II we derive explicit expressions for the position-dependent induced potential and density. In section III numerical calculations of the induced potential in Al and Si are presented, for both random and crystal-symmetry incident directions. In section IV the most relevant conclusions are summarized. Atomic units are used throughout, i. e., $`m_e=e=\mathrm{}=1`$.

## II Theory

We consider a point particle of charge $`Z_1`$ moving in an inhomogeneous system with velocity $`𝐯`$ and impact vector $`𝐛`$, such that

$$\rho ^{ext}(𝐫,t)=Z_1\delta (𝐫𝐛𝐯t).$$
(1)   
For a periodic crystal we can write:

$$\rho _𝐆^{ext}(𝐪,\omega )=2\pi Z_1e^{\mathrm{i}(𝐪+𝐆)𝐛}\delta \left[\omega (𝐪+𝐆)𝐯\right],$$
(2)   
where $`𝐆`$ is a reciprocal-lattice vector and $`𝐪`$ lies in the Brillouin zone (BZ). The external charge induces a density, so that the total density variation of the medium is given by the sum

$$\rho _𝐆(𝐪,\omega )=\rho _𝐆^{ext}(𝐪,\omega )+\rho _𝐆^{ind}(𝐪,\omega ).$$
(3)   
Poisson’s equation allows us to write the potential related to a given density. In our case,

$$\varphi _𝐆(𝐪,\omega )=v_𝐆(𝐪)\rho _𝐆(𝐪,\omega ),$$
(4)   
where $`v_𝐆(𝐪)=4\pi /|𝐪+𝐆|^2`$ and $`\varphi _𝐆(𝐪,\omega )`$ represent the Fourier components of the bare Coulomb potential and the total potential, respectively. Within linear response theory, we can write the total potential in terms of the external charge:

$$\varphi _𝐆(𝐪,\omega )=\underset{𝐆^{}}{}ϵ_{𝐆,𝐆^{}}^1(𝐪,\omega )v_𝐆^{}(𝐪)\rho _𝐆^{}^{ext}(𝐪,\omega ).$$
(5)   
$`ϵ_{𝐆,𝐆^{}}^1(𝐪,\omega )`$ are the Fourier coefficients of the dielectric function, which within the RPA read

$$ϵ_{𝐆,𝐆^{}}(𝐪,\omega )=\delta _{𝐆,𝐆^{}}v_𝐆^{}(𝐪)\chi _{𝐆,𝐆^{}}^0(𝐪,\omega ),$$
(6)   
where $`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )`$ represent the Fourier coefficients of the density-response function of non-interacting electrons,

$`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle \underset{BZ}{}}{\displaystyle \underset{n}{}}{\displaystyle \underset{n^{}}{}}{\displaystyle \frac{f_{𝐤,n}f_{𝐤+𝐪,n^{}}}{\epsilon _{𝐤,n}\epsilon _{𝐤+𝐪,n^{}}+(\omega +\mathrm{i}\eta )}}`$ (9)   
$`\times \varphi _{𝐤,n}|e^{\mathrm{i}(𝐪+𝐆)𝐫}|\varphi _{𝐤+𝐪,n^{}}\varphi _{𝐤+𝐪,n^{}}|e^{\mathrm{i}(𝐪+𝐆^{})𝐫}|\varphi _{𝐤,n}.`$   
Here, the sums run over the band structure for each wave vector $`𝐤`$ in the first BZ, $`|\varphi _{𝐤,n}`$ and $`\epsilon _{𝐤,n}`$ are the one-electron wave functions and energies, and $`f_{𝐤,n}`$ are Fermi factors, $`f_{𝐤,n}=\theta (E_F\epsilon _{𝐤,n})`$. $`\mathrm{\Omega }`$ is the normalization volume.

Combining Eqs. (3), (4) and (5) we obtain the following equation for $`\rho ^{ind}`$ in terms of $`\rho ^{ext}`$:

$$\rho _𝐆^{ind}(𝐪,\omega )=\underset{𝐆^{}}{}\left[|𝐪+𝐆|^2ϵ_{𝐆,𝐆^{}}^1(𝐪,\omega )|𝐪+𝐆^{}|^2\delta _{𝐆,𝐆^{}}\right]\rho _𝐆^{}^{ext}(𝐪,\omega ).$$
(10)   
Substituting Eq. (2) into Eq. (10) and Fourier transforming back to real space yields

$`\rho ^{ind}(𝐫,t)=`$ $`{\displaystyle \frac{Z_1}{\mathrm{\Omega }}}{\displaystyle \underset{𝐪}{\overset{BZ}{}}}{\displaystyle \underset{𝐆}{}}e^{\mathrm{i}(𝐪+𝐆)(𝐫𝐛𝐯t)}{\displaystyle \underset{𝐊}{}}e^{\mathrm{i}𝐊𝐛}e^{\mathrm{i}𝐊𝐯t}`$ (13)   
$`\times \left[|𝐪+𝐆|^2ϵ_{𝐆,𝐆+𝐊}^1(𝐪,(𝐪+𝐆+𝐊)𝐯)|𝐪+𝐆+𝐊|^2\delta _{𝐆,𝐆+𝐊}\right],`$   
where we have made use of

$$\frac{\mathrm{d}^3𝐪}{(2\pi )^3}\frac{1}{\mathrm{\Omega }}\underset{𝐪}{\overset{BZ}{}}\underset{𝐆}{},$$
(14)   
and we have set $`𝐆^{}=𝐆+𝐊`$.

If we consider a definite trajectory of the projectile, only the $`𝐊`$ vectors such that $`𝐊𝐯=0`$ contribute to the sum, and

$`\rho ^{ind}(𝐫,t)=`$ $`{\displaystyle \frac{Z_1}{\mathrm{\Omega }}}{\displaystyle \underset{𝐪}{\overset{BZ}{}}}{\displaystyle \underset{𝐆}{}}e^{\mathrm{i}(𝐪+𝐆)(𝐫𝐛𝐯t)}|𝐪+𝐆|^2{\displaystyle \underset{𝐊}{}}{}_{}{}^{}e_{}^{\mathrm{i}𝐊𝐛}`$ (17)   
$`\times \left[{\displaystyle \frac{ϵ_{𝐆,𝐆+𝐊}^1(𝐪,(𝐪+𝐆+𝐊)𝐯)}{|𝐪+𝐆+𝐊|^2}}{\displaystyle \frac{\delta _{0,𝐊}}{|𝐪+𝐆|^2}}\right].`$   
The prime in the sum over $`𝐊`$ accounts for the $`𝐊𝐯=0`$ condition. Therefore, Eq. (17) gives the density induced by an external charged particle moving in a definite trajectory through a crystal.

In order to get the induced potential we just apply Poisson equation, which yields

$`\varphi ^{ind}(𝐫,t)=`$ $`{\displaystyle \frac{4\pi Z_1}{\mathrm{\Omega }}}{\displaystyle \underset{𝐪}{\overset{BZ}{}}}{\displaystyle \underset{𝐆}{}}e^{\mathrm{i}(𝐪+𝐆)(𝐫𝐛𝐯t)}{\displaystyle \underset{𝐊}{}}{}_{}{}^{}e_{}^{\mathrm{i}𝐊𝐛}`$ (20)   
$`\times \left[{\displaystyle \frac{ϵ_{𝐆,𝐆+𝐊}^1(𝐪,(𝐪+𝐆)𝐯)}{|𝐪+𝐆+𝐊|^2}}{\displaystyle \frac{\delta _{0,𝐊}}{|𝐪+𝐆|^2}}\right].`$   
We can define an average or induced random potential as the mean value over impact vectors $`𝐛`$, which results in the $`𝐊=0`$ term of Eq. (20):

$$\varphi _{random}^{ind}(𝐫,t)=\frac{4\pi Z_1}{\mathrm{\Omega }}\underset{𝐪}{\overset{BZ}{}}\underset{𝐆}{}\frac{e^{\mathrm{i}(𝐪+𝐆)(𝐫𝐯t)}}{|𝐪+𝐆|^2}\left[ϵ_{𝐆,𝐆}^1(𝐪,(𝐪+𝐆)𝐯)1\right].$$
(21)   
The most important contribution to the position-dependent induced potential of Eq. (20) is provided by the $`𝐊=0`$ term. For those directions for which there are not reciprocal vectors satisfying the $`𝐊𝐯=0`$ condition, we have the average potential of Eq. (21). For a few highly symmetric or channeling conditions, non-negligible corrections to the random result are found. The random induced potential exactly coincides with the well-known jellium result when one replaces the inverse dielectric matrix entering Eq. (21) by the inverse dielectric function of a homogeneous electron gas.

The position-dependent stopping power can be obtained directly from Eq. (20), by simply taking into account that

$$\frac{\mathrm{d}E}{\mathrm{d}x}=\frac{Z_1}{v}\varphi ^{ind}|_{𝐫=𝐛+𝐯t}𝐯,$$
(22)   
which gives an expression that exactly coincides with the position-dependent stopping power derived from the knowledge of the imaginary part of the projectile self-energy.

The main ingredient of our calculation of the induced potential and density is the inverse dielectric matrix, which we have evaluated in the RPA by inverting Eq. (6). Hence, at this point we have only considered the average electrostatic interaction between the electrons. We have not found differences in the induced potential by including many-body short-range correlations in the response of the system. The one-electron Bloch states entering Eq. (9) are the self-consistent LDA eigenfunctions of the Kohn-Sham equation of density-functional theory (DFT). We first expand the states in a plane-wave basis and then solve for the coefficients of the expansion self-consistently. The electron-ion interaction is described in terms of a non-local, norm-conserving ionic pseudopotential. The XC potential is computed with use of the energy first calculated by Ceperley and Adler and then parametrized by Perdew and Zunger. We subsequently evaluate the $`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )`$ polarizability and invert Eq. (6) as we sum over $`𝐪`$, $`𝐆`$, and $`𝐊`$ to obtain the induced potential and density.

## III Results

In this section we present the results of our calculations of the induced potential when a charged particle penetrates through aluminum and silicon. The average electron densities of Al and Si are similar: the corresponding free-electron gas (FEG) is characterized by $`r_s=2.07`$ for Al and $`r_s=2.01`$ for Si. However, aluminum is a metal and silicon a semiconductor and the induced potential will exhibit a different behaviour, which we would not obtain on the basis of a FEG calculation.

Although Al is usually regarded as a jellium-like material, inelastic X-ray scattering experiments and theoretical analyses of the dynamical structure factor, stopping power, and hot-electron decay rates have revealed that it is necessary to take into account the full band structure for a proper description of its electronic properties. In this work, Bloch states have been expanded in a plane-wave basis with an energy cutoff of 12 Ry, which corresponds to keeping in this expansion $`100`$ plane waves. The sums over the BZ for both the polarizability $`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )`$ and the induced potential have been performed on $`10\times 10\times 10`$ Monkhrost and Pack meshes. The sum over reciprocal-lattice vectors in the potential has been extended to the first 15 $`𝐆`$ vectors, which corresponds to a cutoff in the momentum transfer of $`2.9q_F`$ ($`q_F`$ is the Fermi momentum). We have included 30 bands in the sums over the band structure for each $`𝐤`$ vector in Eq. (9), which allows us to calculate the induced potential and the stopping power up to velocities of the order of 2 a.u..

Silicon is a covalent crystal which shows strong valence-electron density variations in certain directions. This allows the formation of channels in which the density is very low. For example, the integrated density in the $`110`$ channel varies from $`r_s=1.49`$ at the atomic row to $`r_s=3.37`$ at the center of the channel, i.e., an $`80\%`$. This has an obvious impact on the induced potential and the stopping power. The covalent character of Si imposes a higher cutoff for the Bloch-state expansion than in Al, due to the higher degree of localization of the electronic states in the former material. This also results in Si having bands that are flatter than in the case of Al, so that the number of bands included in the calculation of the polarizability $`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )`$ must be larger in the case of Si if the same energy transfers as in Al are to be included. We have used a cutoff of 16 Ry ($`300`$ plane waves/state), and have included 100 bands in the calculation of $`\chi _{𝐆,𝐆^{}}^0(𝐪,\omega )`$. The sums over the BZ have been performed on $`8\times 8\times 8`$ Monkhorst and Pack meshes, and the sum over reciprocal-lattice vectors has been extended to the first 15 $`𝐆`$ vectors, which corresponds to a momentum-transfer cutoff of $`2.1q_F`$.

We focus on the spatial distribution of the induced potential along the incoming particle trajectory. The $`z`$ coordinate appearing in the figures is always relative to the particle position. The stopping power derived from the slope of the induced potential at the projectile position coincides with the stopping power reported in Refs.  and . An ion of $`Z_1=1`$ is always considered.

In Fig. 1 we have plotted the random potential induced by an ion with $`v=0.6`$ a.u. through Al (a) and Si (b). The solid line represents the calculation for the real solid, whereas the dashed line represents the FEG calculation. For this low-velocity regime (below $`v_F`$), differences between full band-structure and FEG calculations come from the presence of interband transitions and, also, the gap in the case of Si. The splitting of the bands makes the polarization easier and, for this reason, the induced potential at the origin is higher for the real crystal than for the FEG in Al. In the case of Si, the slope of the induced potential at the origin is smaller than in the corresponding FEG, due to the presence of the gap, which, for low-energy transfers diminishes the polarization of the electron system. The gap of semiconductors like Si brings about a non-linear stopping power for low velocities and an induced potential which presents, at these velocities, a lower slope.

In Fig. 2 we have plotted the random potential induced by an ion with $`v=1.5`$ a.u. through Al (a) and with $`v=1.6`$ a.u. through Si (b). Solid and dashed lines represent full band-structure and FEG calculations, respectively. These velocities are well over the plasmon-excitation threshold ($`v_t1.3`$ a.u. in both cases). An oscillatory behaviour appears behind the ion due to plasmon excitation, the wavelength of the oscillations being of $`2\pi v/\omega _p`$ ($`\omega _p`$ is the plasma frequency). As can be appreciated in Figs. 2a and 2b, the wavelength of the random potential is the same in the real crystal and in the corresponding FEG for both materials. This is a consequence of the fact that plasmon contributions to the stopping power are the same for both the real solid and the FEG. However, the oscillations in real Si are more damped than in real Al, due to the shorter lifetime of plasmons in Si which stems from the higher density of bands, thus increasing the decay channels.

In Fig. 3 we have plotted the position-dependent potential, Eq. (20), induced by an ion with $`v=1.5`$ a.u. along the $`100`$ direction in Al (a), and with $`v=2.0`$ a.u. along the $`110`$ direction in Si (b). The impact vector is $`b(0,1,0)`$ for Al and $`b(1,1,0)`$ for Si. $`b`$ is measured in units of the lattice constant $`a_c`$ of the target. Calculations for an impact parameter $`b=0`$ (atomic row) and $`b=0.25`$ (center of the channel) are represented by solid and short-dashed lines, respectively. For comparison, a local-density approximation (LDA) of the induced potential is also displayed in Fig. 3 by a dashed line for $`b=0`$ and a dotted line for $`b=0.25`$. In this approach the position-dependent induced potential is obtained as the potential induced in a FEG with an electron density equal to the average electron density along the projectile path. The slope of the induced potential gives the stopping power in each case. We can appreciate a substantial variation in the slope of Si as we move across the channel, due to the great variation of the density. A lower density at the center of the channel results in a smaller slope and a lower stopping power than in the atomic row and in the random case. Variations in the spatial distribution of the induced potential in Al are fairly small as $`b`$ changes, because the electron density is almost flat in this material. If we calculated the potential at the ion position as a function of $`b`$, we would obtain a qualitatively similar result in both the ab initio and FEG evaluation, i.e., following the local density variation. However, it is clear from Figs 3a and 3b that the spatial distribution of the induced potential is not well described by the LDA calculation. Above all in Si, Fig. 3b, the differences are more dramatic. The oscillations behind the ion have the same wavelength independently of $`𝐛`$ in the crystal calculation, whereas in the LDA calculation the wavelength depends on the average density for each path. Plasmon contributions appear in the $`𝐊=0`$ term of the position-dependent induced potential of Eq. (20). The remaining terms ($`𝐊0`$) are mainly corrections to the single particle e-h contribution, which are due to the presence of density variations away from the projectile path, i.e., the so-called crystalline local-field corrections. For this reason, the wavelength of the induced potential remains unchanged. Furthermore, the stopping-power peak is located at the same position for any velocity direction and for any impact parameter, since plasmons are collective excitations which involve all the electrons of the system. However, in an LDA calculation the stopping-power peak is located at the velocity that corresponds to the average density of the path and not to the average target density. These results lead us to the conclusion that LDA calculations are not suitable for the calculation, at these velocities, of the position-dependent induced potential and stopping power.

## IV Conclusions

We have presented full band-structure calculations of both random and position-dependent induced potentials in Al and Si. The linear-response formalism has been used to obtain the potential and density induced by an external charge penetrating through an inhomogeneous periodic system. The random potential has been evaluated in the RPA for velocities below and above the plasmon-excitation threshold. For low velocities, differences between the FEG and ab initio calculations come from the sensitiveness to the band structure of the target. For higher velocities, we have shown that oscillations behind the ion have the same wavelength in both the FEG and the real crystal, due to the fact that plasmon excitation remains unchanged. Finally, we have investigated the position-dependent potential induced by projectiles incident along the $`100`$ direction in Al and along the widest channel in Si, the $`110`$ direction. Variations in the spatial distribution of the induced potential are more pronounced in the case of Si. Besides, we have shown that the LDA calculation does not properly account, at the velocities under study, for the spatial distribution of the induced potential along the channel.

## V Acknowledgments

We thank P. M. Echenique and A. G. Eguiluz for stimulating discussions. We acknowledge partial support by the University of the Basque Country, the Basque Unibertsitate eta Ikerketa Saila, and the Spanish Ministerio de Educación y Cultura.