Low-energy collective electronic excitations in Pd metal

V. M. Silkin,1,2,3,4 I. P. Chernov,4 Yu. M. Koroteev,2,4,5 and E. V. Chulkov1,2 1Depto. de Física de Materiales and Centro Mixto CSIC–UPV/EHU, Facultad de Ciencias Químicas, Universidad del País Vasco, Apdo. 1072, 20080 San Sebastián, Basque Country, Spain 2Donostia International Physics Center (DIPC), P. de Manuel Lardizabal, 4, 20018 San Sebastián, Basque Country, Spain 3IKERBASQUE, Basque Foundation for Science, 48011 Bilbao, Spain 4Tomsk Polytechnical University, pr. Lenina 30, 634050 Tomsk, Russia 5Institute of Strength Physics and Materials Science, Russian Academy of Sciences, pr. Academicheskii 2/1, 634021 Tomsk, Russia Received 17 August 2009; revised manuscript received 17 November 2009; published 21 December 2009


I. INTRODUCTION
Single-particle electron-hole and collective electronic excitations in solids have been investigated for many years on the basis of the momentum-and frequency-dependent dynamical electron-density response. 1Basic knowledge on it was gained from a study of a free-electron gas ͑FEG͒ model in which an electron valence density of real solids is represented by an average electron density, characterized by an interelectron distance-r s parameter. 2For years this model was used in many studies, e.g., in order to obtain information on the dynamical exchange-correlation effects.However, direct comparison of predictions of FEG theories with measurements of real solids was rather difficult.In particular, very important and frequently dominant effects in dynamical dielectric properties of solids are introduced by electron band structure which in real solids can substantially deviate from the FEG one.For instance, the presence of more than one energy band can be a source of numerous interband transitions presented in real solids.Sometimes these transitions produce only little effect on the excitation spectra, introducing, e.g., a finite width in measured plasmon peaks ͑which must be zero up to certain momenta according to the FEG models, i.e., plasmon would live once excited for ever͒.However, frequently these transitions produce dramatic impact on the dielectric and optical properties of solids.A wellknown examples of this are, e.g., appearance of an additional low-energy plasmon peaks in layered materials [3][4][5][6][7][8][9] and carbon nanotubes, [10][11][12] dramatic redshift of plasmon frequency in Ag, [13][14][15] or a negative momentum dependence of plasmon peak in Cs. 16,17 For three-dimensional ͑3D͒ solids the FEG predicts that for collective excitation of an electron system or a plasmon mode, there is some threshold which is determined by average valence density.As in real metals r s is in the 1-6 a.u.range, these frequencies normally arise in the ϳ3-20 eV interval.Hence these plasmons can not participate directly in low-energy dynamical processes occurring in an electron system in the vicinity of the Fermi surface.On the other hand, in the 1950s it was pointed out that once in an electron system more than one energy bands with different Fermi velocities cross the Fermi level there exists possibility of appearance of an additional low-energy mode with peculiar sound-like dispersion-acoustic plasmon ͑AP͒.At small momenta q the energy AP of this mode disperses linearly with q as AP = u AP • q, 18 i.e., AP tends to zero as q → 0.Here the coefficient u AP is a group velocity of the acoustic plasmon and it is determined by the Fermi velocity of the carriers in the energy band with lower Fermi velocity, u F slow , being u AP Ϸ u F slow .A similar idea of existence of an acoustic-like mode was proposed in the 1960s for explanation of elevated temperatures for transition to the superconducting state, T c , in transition metals, [19][20][21] where two kinds of carriers in both slow ͑d bands͒ and fast ͑s-p bands͒ are usually presented. 22he idea of acoustic plasmon existence was regularly evoked for explanation of elevated T c in superconductors with unusually T c ͑Refs.23-25͒ which can be explained by the fact that the presence of an acoustic plasmon might dramatically affect the low-energy dynamical screening properties of an electron system.Nonetheless, to the best of our knowledge, this kind of plasmons has been shown to exist in any bulk material neither experimentally nor by ab initio calculations.Only very recently a detailed ab initio calculation of dynamical momentum-dependent response demonstrated its possible existence in the MgB 2 compound for the momentum trans-fers along the direction perpendicular to the basal hexagonal boron planes. 26,27Note that a different type of acoustic plasmon characteristic for layered conductors can be caused by the out-of-phase motion of carriers in the neighboring layers.][30][31] In the present paper we demonstrate that the first kind of a low-energy plasmon with acoustic-like dispersion can exist in a real 3D metal where more than one energy band crosses the Fermi level.Here, in a particular case of Pd it is shown that this acoustic plasmon is a truly 3D phenomenon and exists in all three high-symmetry directions. 32Note that the excitation spectra of palladium have been a subject of intense study in the past.2][43][44][45][46][47] Besides some spread in the reported data one can observe that dielectric properties in Pd are determined by a dominant bulk plasmon mode located in the 7-8 eV energy range and several other less-pronounced peaks in energy-loss spectra due to numerous interband transitions. 38,48At the same time we are not aware of any similar detailed study of excitation spectra of this metal which would be focused on the low-energy ͑less than ϳ1 eV͒ region.
The outline of the paper is as follows.In Sec.II we present the details of the ab initio calculation of dynamical momentum-dependent dielectric function in 3D solids.Description and discussion of the obtained results are given in Sec.III.Finally, the main conclusions are summarized in Sec.IV.Atomic units ͑ប = e 2 = m e =1͒ are used throughout unless otherwise is stated.

II. CALCULATION METHOD
Within the first Born approximation the inelastic scattering cross section of x-rays and electrons for momentum transfer q and energy transferred to the system is proportional to the dynamical structure factor S͑q , ͒ of manyelectron system.S͑q , ͒ is related by the fluctuation-dissipation theorem 1 to the dielectric function ⑀͑r , rЈ , ͒ which, in general, is a nonlocal energy-dependent function.For a periodic solid S͑q , ͒ is where ⍀ represents the normalization volume, G's are the reciprocal lattice vectors, and Im͓⑀ G=0,G Ј =0 −1 ͑q , ͔͒ is the energy-loss function related to the Fourier coefficients of the density-response function

͑2͒
Here G ͑q͒ =4 / ͉q + G͉ 2 is Fourier transform of the Coulomb potential.Within linear-response theory relates the electron density, n ind ͑r , ͒ induced in the system, to an external perturbation potential, V ext ͑r , ͒, through the equation In the framework of time-dependent density-functional theory, 49,50 ͑r , rЈ , ͒ obeys the integral equation ͑r,rЈ͒ = 0 ͑r,rЈ͒ + ͵ dr 1͵ dr 2 0 ͑r,r 1 ,͓͒͑r 1 − r 2 ͒ where 0 ͑r , rЈ , ͒ is the density-response function for a noninteracting electron system, ͑r − rЈ͒ is the bare Coulomb potential, and K xc ͑r , rЈ , ͒ accounts for dynamical exchangecorrelation effects.For a periodic crystal the integral equation ͑4͒ transforms in the reciprocal space into a matrix equation where the factor of 2 accounts for spin, the sums over n and nЈ run over the band structure for wave vectors k in the first Brillouin zone ͑BZ͒, f nk is the Fermi distribution function, nk and nk ͑r͒ are Bloch eigenvalues and eigenfunctions, respectively, of the Kohn-Sham Hamiltonian, and a positive infinitesimal.Numerical evaluation of GG Ј 0 ͑q , ͒ matrix is the most time-demanding part of such calculations and one should put some finite value 51 for which introduces an additional broadening in the calculated spectra.In practice, instead of direct use of Eq. ͑6͒ the matrix G,G Ј 0 ͑q , ͒ can be evaluated by performing calculations at imaginary 52 or complex 53 frequencies with subsequent analytical continuation to the real axis.An efficient approach based on the tetrahedron method was realized as well. 54Here we use an alternative route consisting in calculating, at the first step, of the spectral function matrix S G,G Ј 0 ͑q , ͒ ͑Refs.55 and 56͒ according to expression Present calculations have been performed for fcc Pd with the use of experimental lattice parameters a c = 7.3512 a.u.The one-particle energies nk and wave functions nk ͑r͒ were obtained as a self-consistent solution of the Kohn-Sham equations with the use of exchange-correlation potential in the form of Ref. 57.The electron-ions interaction was described by a nonlocal norm-conserving ionic pseudopotential.58 Wave functions nk ͑r͒ were expanded in plane waves up to a kinetic-energy cutoff of 70 Ry.Up to 50 reciprocal vectors G were included in the Fourier expansions of 0 , , and ⑀ matrices.Hence, the local-field effects 59,60 were incorporated in the evaluation of the energy-loss function, Im͓⑀ G=0,G Ј =0 −1 ͑q , ͔͒, through the inclusion of nondiagonal matrix elements of G,G Ј 0 ͑q , ͒ in matrix Eq. ͑5͒.In Eq. ͑7͒ the sum over k includes a 144ϫ 144ϫ 144 sampling, that corresponds to 1 492 992 points in the BZ.The sums over n and nЈ included energy bands up to energy of 50 eV above the Fermi level.The S G,G Ј 0 ͑q , ͒ matrices were calculated at the discrete mesh of energies ranging from 0 to 50 eV with a step of 0.01 eV.In the numerical calculations the ␦ function in Eq. ͑7͒ was replaced by a Gaussian with a broadening parameter ␥ = 0.025 eV.We solve Eq. ͑5͒ using a random-phase approximation, i.e. taking K xc =0.It was demonstrated 61 that for transition metals this is a rather good approximation for the description of low-energy electronic excitations, i.e., in an energy range of primary interest here.

III. CALCULATION RESULTS AND DISCUSSION
The calculated band structure as a basic ingredient of the dielectric matrix evaluation is shown for Pd in the vicinity of the Fermi level in Fig. 1.This band structure is in good agreement with previous theoretical studies [62][63][64][65][66] although some small quantitative differences in energy positions of energy bands resulted from the use of different calculation methods can be noted.In the figure one can see how three energy bands cross the Fermi level.As will be shown below this produces strong effect on the low-energy response properties in palladium.
In Fig. 2 we present the calculated dielectric function for Pd, ⑀͑q , ͒ϵ⑀ G=0,G Ј =0 ͑q , ͒, and imaginary part of inverse dielectric matrix, Im͓⑀ G=0,G Ј =0 −1 ͑q , ͔͒-a so-called energyloss function-at q = 0.0712 a.u.−1 along the ͗100͘ symmetry direction as a function of energy for small .In the figure one can clearly see that in this energy region these quantities are qualitatively different from that calculated with the use of the FEG model.In particular, one can observe that Re͓⑀͔ instead of one zero-crossing in this energy region ͑as it has Re͓⑀͔ calculated within the FEG model͒ indeed has three such zero crossings.This is explained by the fact that there are two dominant peaks in Im͓⑀͔: a peak A centered at = 0.09 eV and a peak C at = 0.54 eV. 67As Im͓⑀͔ and Re͓⑀͔ are interrelated via the Hilbert transformation, this makes Re͓⑀͔ to reach zero three times.Whereas the zero crossings at 0.25 and 0.55 eV do not produce any peak structure in the energy-loss function because of coincidence with the local peaks at fairly the same energies in the Im͓⑀͔ part, the zero crossing of Re͓⑀͔ at 0.32 eV ͑highlighted in Fig. 2 by a circle͒ combined with a local minimum in Im͓⑀͔ at the same energy does produce a sharp peak in Im͓⑀ −1 ͔.The presence of this peak is a signature of existence of a well-defined collective mode at this energy.Similar calculations performed for other small momenta in all three symmetry directions demonstrate that the energy positions of local maxima A, B, and C in Im͓⑀͔ and corresponding minima between them disperse quasilinearly with ͉q͉.This behavior is explained by the next observation: at small ͉q͉ the energy for intraband transitions in Eq. ͑7͒ can be expressed as where u nk = ‫⑀ץ‬ nk / ‫ץ‬k is the group velocity in band n at a given wave vector k.As the maximum number of such intraband transitions in a given band n occurs with the u nk components close to its maximal value, the corresponding structure in Im͓⑀͔, which is a measure of a number for electron-hole ͑e-h͒ transitions, is peaked at the high-energy side of this e-h region. 1 As a result, the slope of dispersions of maxima in Im͓⑀͔ is determined by the corresponding maximal Fermi velocity components in a given direction for each energy bands crossing the Fermi surface.
Figure 3 shows density of states of the energy bands at the Fermi level versus their energy and magnitude of the group velocity components.One can note significant difference between the group velocity components in the s-p and d bands.
As the positions of peaks in Im͓⑀͔ disperse fairly linearly with the magnitude of q, the same do the peaks, corresponding to the acoustic plasmon, in Im͓⑀ −1 ͔.The collection of the peak positions in Im͓⑀ −1 ͔ for transfer momenta in the three main symmetry directions in Pd, is presented in Fig. 4. In the figure the dispersion of the acoustic plasmons is shown by solid lines, whereas the peak positions corresponding to interband transitions are presented by dashed lines.One can see that whereas along the ͗100͘ and ͗111͘ symmetry directions there exist only one branch corresponding to the acoustic plasmon, along the ͗110͘ one there are two such branches.The reason for this is explained in Fig. 5 where dielectric ⑀͑q , ͒ and energy-loss Im͓⑀ −1 ͔ functions are presented for momentum transfer q along the ͗110͘ symmetry direction with q = 0.0504 a.u.−1 .Unlike Fig. 2, here one can observe three clear local maxima A, B, and C in Im͓⑀͔ which are explained by the fact that along this direction the strong anisotropy of the group velocity components of all three bands crossing the Fermi level occurs.The presence of the three peaks in Im͓⑀͔ in this energy interval leads to the appearance of five zero crossings in Re͓⑀͔.As a result, the two of them ͑denoted by circles in Fig. 5͒ in combination with the local minima in Im͓⑀͔ produce two clear AP 1 and AP 2 peaks in Im͓⑀ −1 ͔ corresponding to two distinct acoustic plasmon modes.The calculated group velocities of acoustic plasmons in Pd for the three main symmetry directions are presented in Table I.One can see that these values are of two orders of magnitude higher than the sound velocity in metals.Hence these acoustic plasmons in Pd can not interact directly with an acoustic phonon mode. 25ere we would like to stress the importance of achieving of convergency in the calculations in order to obtain the unambiguous results regarding the low-energy dielectric properties of real metals.In Fig. 6, as an example, we demon- strate how the evaluated dielectric function at a small momentum transfer q crucially depends on value of the broadening parameter ␥ employed in the evaluation of the spectral function via Eq.͑7͒.Thus, at the smallest ␥ value a well-defined AP peak corresponding to the acoustic plasmon appears at = 0.18 eV.When the ␥ value increases the energy of this peak is blueshifted and its width increases significantly.As ␥ increases further more this mode even disappears.This is accompanied by dramatic changes in the dynamical dielectric function.Note that at large ␥, instead of three zero crossings in this energy range, the ab initio Re͓⑀͔ part crosses zero only once and has behavior similar to that expected from the FEG model.In this respect quite instructive situation has occurred with the MgB 2 compound.After discovery of superconductivity in this material, 68 it was suggested on base of simple tight-binding calculations that an acoustic plasmon could exist there. 24Nevertheless, subsequent ab initio calculations 69,70 did not confirm that claim.However, very recently this plasmon mode has been found in very detailed ab initio calculations. 26,27he results presented here signal about paramount importance of achievement of convergency in the evaluation of the dynamical momentum-dependent dielectric function and consequently dynamical Coulomb interaction in real solids.In particular, it would be of great interest to investigate possible relation to the evaluated superconducting properties of materials within an "ab initio" approach, 71,72 where the dynamical Coulomb interaction is a key ingredient.Thus the difference of order of 10 K in the predicted T c in MgB 2 between two calculations 71,72 which used different approximations for the description of dynamical Coulomb interaction points out to the importance of inclusion of the proper dynamical screening.
This work presents for the first time an unambiguous ab initio evidence of existence of an acoustic plasmon mode in an elemental 3D bulk metal such as Pd.Recent successful experimental observation in the energy-loss experiment 73 of an analogous acoustic mode 74,75 at a metal surface in a nice agreement with ab initio predictions 73 points out to the principal possibility of experimental detection of this kind of excitations in bulk as well.Regardless of a possible relevance of acoustic plasmons for superconductivity, 76 the present findings amplify our basic knowledge of the low-   energy dynamical electronic screening properties of solids which can qualitatively differ from the expectations of simplified theories.

IV. CONCLUSIONS
We have presented a theoretical study of electronic excitation spectra in Pd with the focus on the low-energy domain.The calculations were performed with the full inclusion of the electron band structure obtained within selfconsistent pseudopotential approach.We demonstrate that the presence in Pd of two kinds of carriers ͑in s-p and d bands͒ at the Fermi level produces dramatic modification of the excitation spectra at these energies in comparison with FEG model predictions.Thus at low energies, additionally to several well-defined features in the loss spectra due to numerous interband transitions in this metal, at small momenta a peculiar plasmon mode, resulting from the presence of fast ͑in the s-p band͒ and slow ͑in the d bands͒ carriers, with characteristic sound-like dispersion-acoustic plasmon-is predicted.This plasmon posses a strong directional dependence on momenta: whereas for momentum transfer along the ͗100͘ and ͗111͘ symmetry directions it arises as an single mode, along the ͗110͘ symmetry direction two acoustic modes with different slopes appear.
Finally, our analysis points out that one should be very careful in performing such kind of evaluation.In particular one should pay special care regarding the convergence in the k sampling over the Brillouin zone in the evaluation of the density-response function o of noninteracting electrons.

FIG. 1 .
FIG. 1. ͑Color online͒ Calculated band structure of Pd in ity of the Fermi level, E F , along some high-symmetry directions of the Brillouin zone.The energy bands crossing the Fermi level are highlighted by colors according to their main symmetry-red ͑gray͒ and blue ͑dark gray͒ colors stand for the s-p and d bands, respectively.The energies are according to the Fermi level shown by dotted horizontal line.

FIG. 3 .
FIG. 3. ͑Color online͒ DOS of the energy bands in vicinity of the Fermi level in Pd versus its energy and group velocity components along three main symmetry directions: ͑a͒ ͗100͘; ͑b͒ ͗110͘; ͑c͒ ͗111͘.Peaks denoted as A and C ͑B͒ correspond to DOS in d ͑s-p͒ bands.The energies are according to the Fermi level.

34 FIG. 4 .
FIG.4.͑Color online͒ Dispersion of peaks in calculated ab initio energy-loss function of Pd for the momentum transfer q along the ͑a͒ ͗100͘, ͑b͒ ͗110͘, and ͑c͒ ͗111͘ symmetry directions.Acoustic-plasmon dispersions are depicted by solid lines, whereas dashed lines show dispersion of less-pronounced energy-loss peaks having its origin in interband transitions.

FIG. 6 .
FIG.6.͑Color online͒ Dependence of imaginary ͑upper panel͒ and real ͑middle panel͒ parts of the dielectric function ⑀͑q , ͒, as well as the corresponding energy-loss function, Im͓⑀ −1 ͑q , ͔͒, on the broadening parameter ␥ employed in the evaluation of the S G,G Ј 0 ͑q , ͒ matrix.Momentum q is along the ͗100͘ symmetry direction with ͉q͉ = 0.036 a.u.−1 .Solid, dashed, and dotted lines stand for results obtained with ␥ = 0.025 eV, ␥ = 0.050 eV, and ␥ = 0.100 eV, respectively.The corresponding FEG results, obtained with the use of r s = 3.5 a.u., are presented by thin solid lines.

TABLE I .
Group velocities ͑in a.