Band Structure versus Dynamical Exchange-Correlation Effects in Surface Plasmon Energy and Damping: A First-Principles Calculation

176801-1 A first-principles parameter-free calculation that includes full three-dimensional band structure and dynamical exchange correlations is reported for the dynamical surface response and surface plasmon (SP) on a simple metal prototype surface Mg(0001). We demonstrate that band structure effects have a more profound impact on the SP characteristics than dynamical exchange correlations. A comparison with jellium and one-dimensional potential evaluations shows that the band structure is of paramount importance for the correct description of the SP linewidth and also leads to a better description of the SP energy dispersion. The inclusion of the exchange-correlation kernel results in a better agreement with experimental data. We show that lateral crystal local field effects have a negligible impact on the SP properties. Significant anisotropy is predicted for the SP linewidth.

Since the introduction of the surface plasmon (SP) by Ritchie [1] remarkable progress has been achieved in the understanding of collective electronic excitations at metal surfaces [2 -4].These excitations play an important role in such areas as surface dynamics [2 -5] including oneparticle (electron and hole) dynamics [6], SP microscopy [7], SP resonance technology [8], as well as in photonic applications [9].Conceptually the SP, a fundamental surface mode, can be traced to the peaks of the imaginary part of the surface response function [4] derived from the density response function of the interacting electron system [10].The density response includes both bulk and surface electrons and represents a very complex interplay between one-particle electron states (band structure) and dynamical many-body correlations [10] that affects both SP energy and linewidth.Thus one can expect that ab initio calculations of SP properties on some surface well studied experimentally can shed light on the relative importance of these two factors in the description of SP.
Up until now the most advanced theoretical calculations, based on a semiinfinite jellium model and timedependent local density approximation (TDLDA), accounted for the negative SP dispersion at small momenta for simple metal surfaces [11] which has been a subject of great debates for a long time [2,12].The calculations also led to the qualitatively good description of the SP dispersion at large momenta.However, these calculations showed that the evaluated SP energies are systematically too high compared to experimental data [11,13,14].Even for such simple metal as Mg, for which one can expect a minimal impact of band structure, the discrepancy of the calculated SP energy with measured results exceeds 3 times the experimental error bar [14].Also the measured dispersion is significantly flatter than the theoretical one [11,13,14].Another existing problem is the description of the SP damping, an important characteristic condition-ing, for example, the propagation of a SP polariton in surface plasmon subwavelength optics [15] or the functioning of proposed ''superlenses'' [16].Here the situation is even qualitatively wrong.The jellium models predict that for small two-dimensional (2D) momenta q the SP linewidth linearly decreases to zero [2,4], while experiments give a final linewidth value of the order of 1 eV at q 0 [11,13,17].Furthermore, the calculated linewidth increases with q much faster for small momenta and much slower for large q than the experimental one does [4].
The difficulty in the description of a SP stems from a gradual transition of main driving forces from the longwavelength limit to large momenta.In the longwavelength limit the dielectric response is controlled by bulk properties, while for large momenta the main effect is due to the surface dielectric response.Thus for the proper description of the dispersion of the SP energy and width it is necessary to include both the bulk and surface electronic structure on the same footing.
There were attempts to incorporate some aspects of band structure into the calculation of surface plasmons for real metals.These attempts include the stabilized jellium model [4], one-dimensional (1D) variation of the crystal potential in the direction perpendicular to the surface [14], and s ÿ d polarization model for Ag surfaces [4].The three-dimensional (3D) nature of adsorbate atomic layers was considered by Ishida and Liebsch [18] in the study of electronic excitations in thin alkali metal layers adsorbed on jellium substrate.
In parallel with the band structure the exchangecorrelation effects were shown to be important for the description of bulk collective excitations in simple metals [19].For surface collective excitations the incorporation of dynamical exchange correlations might be even more important than for bulk materials, since unlike the bulk In this Letter we present the first parameter-free ab initio calculation of the SP energy and linewidth for a real metal surface.The relative impact of band structure and dynamical exchange correlations on the SP is studied for Mg(0001).This surface can be viewed as a test one for the comparison of theoretical and experimental results because of the best quality of the Mg single crystal surface with respect to other simple metals [17].We also calculate the SP characteristics by using the jellium model and 1D potential that varies in the direction z perpendicular to the surface [20].This potential describes main features of the Mg(0001) surface electronic structure: an energy gap and a surface electronic state at the surface Brillouin zone (SBZ) center.The results obtained show that even for such a nearly free-electron metal as Mg an excellent agreement with the experimental SP energy and linewidth in a large range of 2D momenta is found if both the bulk and surface band structure on the same footing together with dynamical exchange correlations are taken into account.The inclusion of band structure into the theory is crucial for the description of the correct behavior of the SP linewidth as well as of its anisotropy.
Information on a surface plasmon can be obtained from the imaginary part of the surface response function gq; !defined as [4,21] gq; !Z drn ind r; !e qz ; (1) where r r; z, q is a 2D momentum parallel to the surface, and n ind r; ! is the density induced at the crystal surface by an external potential of the form e qz e iqr e ÿi!t : ( Here we take explicitly into account the 3D nature of the surface response without any assumption of a 2D translation invariance of the crystal, as it is a common practice in the jellium based models [4].In this case g depends on both the value and direction of q.
In time-dependent density functional theory [10] n ind and V ext are related by the equation The response function of interacting electron system, r; r 0 ; !, satisfies the integral equation [10] where 0 is the response function of noninteracting electron system, v is the Coulomb interaction potential, and K xc accounts for dynamical exchange-correlation effects.
At the surface the translational symmetry in the perpendicular direction is broken.To deal with this problem we employ the supercell geometry with slabs containing 16 Mg atomic layers and separated by vacuum intervals [22].In this case matrix elements of 0 are given by je ÿiqrGr j kq;n 0 i h kq;n 0 je iqrG 0 r j k;n i; ( where the sum over n; n 0 runs over bands up to energies of 30 eV above E F and the sum over k in the SBZ includes 7812 points; f k;n are the Fermi factors, E k;n ( k;n ) represent the one-particle energies (wave functions) of a Hamiltonian evaluated in the local density approximation with the use of ab initio norm-conserving pseudopotential [24].The exchange-correlation potential was taken in the form of Ref. [25].k;n were expanded in a plane wave basis set with a kinetic-energy cutoff of 12 Ry.S is the normalization area.Expansion of 0 includes 300 G vectors.Combining Eqs. ( 2) -( 5) we finally obtain Note, although in this formula the G;G 0 matrix elements with G G 0 0 enter only, the full 3D nature of is implicitly included via the evaluation of Eq. ( 4).
The key and most time consuming ingredient in the evaluation of g is the calculation of 0 .In direct methods of calculation of 0 G;G 0 one needs to use some finite value of the broadening parameter .As we are interested in the evaluation of the SP linewidth, we avoid the direct use of Eq. ( 5), doing faster calculations of the spectral function matrix with subsequent evaluation of the imaginary and real parts of 0 G;G 0 following Ref.[26].We have also checked whether lateral crystal local field effects are important for the description of the SP.For that we calculated 0 G;G 0 retaining reciprocal lattice vectors G f0; 0; G z g only.With this 0 G z ;G 0 z the surface response function was computed by using Eqs.( 4) and ( 6).The obtained energy and linewidth values are very similar to those found with the use of full 3D 0 G;G 0 : indeed these quantities are indistinguishable from those presented in Figs. 1 and 2. As these calculations are much faster than ones based on 3D G's, it allows us to drastically reduce the computation time retaining all important 3D effects that are included in Eq. ( 5) due to 3D Bloch functions k;n and one-electron energies E k;n .
The calculated SP dispersion for two symmetry directions in SBZ is shown in Fig. 1 together with the experi-mental data [17].We also show the results obtained from the jellium and 1D potential models.In the latter evaluations we have used much thicker slabs that allow us to describe the SP at very small q and study the dependence of the results obtained on the slab thickness.Thus, comparing these three models one can discriminate the relative role of such factors as modulation of one-electron potential in the direction perpendicular to the surface and full inclusion of 3D band structure effects.For these models we have done calculations with three kinds of exchange-correlation kernel K xc [see Eq. ( 4)]: random phase approximation (RPA), when K xc 0; TDLDA; and the momentum dependent kernel parametrized by Corradini, Del Sole, Onida, and Palumno (CDOP) [27].One can see in Fig. 1 that the best agreement with the experiment is obtained from an ab initio calculation with the CDOP kernel, which lowers the dispersion curve for high momenta compared to the TDLDA one.This influence of the CDOP kernel with momenta is similar to that obtained in the study of a bulk plasmon dispersion [28].The RPA significantly overestimates the surface plasmon energy at large q.Similar dependence on K xc (not shown in Fig. 1) is also obtained from evaluations within both the jellium and 1D potential models.One can see that the use of the 1D potential significantly reduces the SP energy at small q with respect to the jellium model.Nevertheless, this model does not reproduce well the experimental data at large q, approaching the jellium model results, that emphasizes the increasing importance of the 3D band structure effects with the increase of momentum.Note, that despite the explicit inclusion of lattice effects into ab initio theory the calculation leads to almost perfect isotropy of the SP dispersion that is also the case in the jellium and 1D potential evaluations.
The calculated SP linewidth dispersion is shown in Fig. 2. Again, the best agreement with the experiment [17] is obtained from ab initio calculation with the CDOP kernel while the RPA kernel significantly lowers the linewidth value.For small q the agreement with the experimental linewidth is worse than in the case of the plasmon energy.It is explained by difficulties in the extraction of the linewidth values from the calculated

Ab initio,(
) ΓΜ ΓΚ Jellium model 1D model potential Experiment, Ref. [17] Linewidth (eV) Linewidth (eV) surface response function for these momenta, where slab effects are manifested in a complicate manner [4,23].Note that ab initio theory gives a negative slope for the linewidth dispersion at small q in agreement with the measured data [17] and for larger q the theory results in a clear dependence of the linewidth on the momentum direction.In contrast to the energy dispersion the 1D potential evaluation does not improve the jellium model results, leading to qualitatively wrong behavior of the linewidth at small q.
We have presented the first ab initio calculation of dynamical surface response with full inclusion of band structure and dynamical exchange correlations using Mg(0001) as an example.We have shown that even for a nearly free-electron metal such as Mg band structure effects are crucial to understand SP characteristics, dispersion and linewidth.In particular, 3D band structure effects are crucial for the description of the SP linewidth, whereas jellium and 1D potential models fail to reproduce the experimental linewidth for small and large momenta.These effects account for the significant anisotropy predicted for the SP linewidth.They are also important for the energy dispersion of the SP.The inclusion of the exchange correlations improves the theoretical results leading to excellent agreement with the experimental data [17].We have found that lateral crystal local field effects have a negligible impact on the SP properties.Our conclusions about the relevance of band structure can be extended to systems with more complicated electronic structure than Mg, for instance, to Ag surfaces, which are of particular interest due to an unusual energy and dispersion of a SP, [3,4] and to Pd surfaces [3].Already phenomenological inclusion of d bands of Ag into the surface response via s ÿ d polarization model led to better description of the SP dispersion on silver surfaces [4].One therefore can expect that the exact inclusion of band structure into the calculation of dynamical surface response can account for unusual features of the SP characteristics for these metals.The presented results demonstrate a feasibility of ab initio calculations of collective excitations on surfaces of different materials such as simple, noble, and transition metals as well as collective excitations in free standing and adsorbed thin films as was illustrated recently for metallic monolayers [29].
density experiences fast variation at short distances from bulk values to zero in a vacuum.

FIG. 2 ( 4 FIG. 1 (
FIG. 2 (color online).(a) Surface plasmon linewidth dispersion for Mg(0001) calculated with CDOP kernel.Ab initio data for the ÿ M (ÿ K ) direction of SBZ are shown by the filled (open) diamonds.Solid (long dashed) line is the corresponding best data fit.Dashed (dotted) line presents the jellium model (1D model potential) dispersion.Stars indicate experimental data [17].(b) Ab initio surface plasmon linewidth dispersion calculated with the RPA (circles), TDLDA (squares), and CDOP (diamonds) kernels.Open (full) symbols correspond to the ÿ K (ÿ M ) direction.The long dashed (solid) lines are the best fit for the open (full) symbols.