Surface collective oscillations of metal clusters and spheres: Random-phase-approximation sum-rules approach

Using the once and thrice energy-weighted moments of the random-phase-approximation strength function, we have derived compact expressions for the average energy of surface collective oscillations of clusters and spheres of metal atoms. The L =0 volume mode has also been studied. We have carried out quantal and semiclassical calculations for Na and Ag systems in the sphericaljellium approximation. We present a rather thorough discussion of surface diffuseness and quantal size effects on the resonance energies.


I. INTRODUCTION
The study of surface modes of metal spheres has drawn some attention in the past' and is the subject of a renewed interest.
Very recently, some experimental information on surface modes ' and electric polarizabili- ties ' of similar systems, namely metal clusters, has been obtained.It seems that the gross features of the experi- mental findings are well reproduced by a self-consistent time-dependent local-density approximation (TDLDA), in which the small metal particles are modeled by a self- consistent spherical-jellium model." Unfortunately, these TDLDA calculations require a substantial numerical effort that is very likely quite formidable.Moreover, some of the physics might be ob- scured by the unavoidable numerical apparatus.These facts prompted Bertsch and Ekardt'" to use in this con- text the random-phase-approximation (RPA) sum-rule approach, which has proven to be extremely useful in describing giant resonances in atomic nuclei.' ' It is worth emphasizing that RPA is the small-amplitude limit of TDLDA and, consequently, both methods are com- pletely equivalent for studying the energy spectrum of these systems.
The purpose of the present work is to extend that of Ref. 14 to all multipole surface modes, and also to present some results for the L =0 volume mode.We will rely heavily on similar work recently done on isoscalar nuclear giant resonances.'   The basic hypotheses of our model are the following.
(i) The mean field created by the positive ions is obtained in the spherical-jellium approximation.
(ii) The valence electrons are described in a local- density approximation (LDA), either by a Kohn-Sham  (KS) or by a Thomas-Fermi (TF) method.Of course, the RPA sum rules will be deduced in the KS framework.
(iii) The external field acts only on the valence elec- trons.
These hypotheses are currently used in microscopic calculations of metal spheres and metal clusters." Thus, the method we use in this work applies to both systems.This paper is organized as follows.We shall describe the RPA sum rules in Sec.II.The KS and TF models we have used will be briefly discussed in Sec.III.The nu- merical results for the surface modes and the L =0 volume mode will be presented in Sec.IV.Finally, we shall make some concluding remarks in Sec.V.

II. RPA SUM RULES
A. General description RPA sum rules are described in great detail in Refs.15 and 16.We refer the interested reader to these references where, besides a clear exposure of the fundamental re- sults, one may find references to earlier work in this field.
Sum rules mk are moments of the strength function s(E) = )~6(F-E")l(n l g lp) l '-, where the sum (integral in the case of continuum spec- trum) extends over all the excited states of the system.Q is the external field acting on the system; E", ln ), and l P ) are the excitation energies, the excited states and the ground state (g.s.) of the system, respectively.The mo- ment of order k is defined as m, = f E"s(E)«=g '.~""l(.Igloo&l'.From these moments, the average energy and variance of 39 8247 1989 The American Physical Society the strength can be obtained: Obviously, a direct computation of mk from Eq. ( 2) is hopeless in most practical cases because one should know the whole spectrum.However, odd moments of S(E) can be obtained with RPA precision as expectation values on the KS ground state lP & of suitable operators.
Defining Ek --(mi, /mk z)', it has been shown that' with similar expressions for ~, 7 ], and ~2.Developing the expectation value (7) in powers of ri, the above expressions allow for a straightforward (but usually cumbersome) evaluation of m 3. B.
we get The double commutator entering the definition of m & is easy to evaluate when Q only depends on the position.This is not quite so for m 3, which is more easy to obtain by scaling the Slater determinant lP&.Indeed, defining the scaled wave function u= and p is the g.s.equilibrium density.Equation ( 12) shows that the velocity field is irrotational.Since b, Q=0, we which is peaked at the surface.The general expression for p2 is which in the present case becomes 2 VQ Vpi. (15) Equations ( 13) and ( 15) explicitly show that Jp, dr=0, i=1,2.It is worth noting that the condition divu=0 is the origin of the surface character of the oscillations gen- erated by Q.That will be more clear when we discuss the where p(r) is the KS g.s.equilibrium density.For L = 1, it is proportional to the electron number X: p"(r)=&4"lplk"&=p+qp, +g p, +. . ., ~"(r)=&y"lr"ly"&=~+~-~, +~r, + where and for L=2, proportional to the square of the (rms) ra- dius mi, (L =2)= fi 10K &r '&,   2m 4m   p=g 5(r -r, ), &r'&-: -J drr'p(r) . [m 3 (j-e ) ]. Altogether, both Coulomb contributions will be referred to as m3(C).It has been proven' that pure volume terms, i.e. , terms that only depend on the electron density, like the Coulomb-exchange contribution in the Slater approximation or the correlation energy in the Lang-Kohn approximation, do not contribute to m3.That can be easily seen using the technique we will out- line below to obtain m3(e-e ) and m3( j-e ).
A, (r) represents a centrifugal kinetic-energy density, and in the A' -order TF approximation one has A, =2~/3.
It is worth mentioning that in the dipole case, only m3(j-e) contributes to the "restoring force" represented by m 3. Indeed, the electron-electron kinetic and Coulomb contributions vanish because the operator r Y, 0 corresponds to a translation of the electron cloud as a whole.Only the translational symmetry-breaking jellium field gives a nonzero contribution.
C. The sum rules m "m "aud m, for Q = r This field will generate monopole volume oscillations.
It has been thoroughly studied in the nuclear case.' It is easy to show that while the second integral reads m, =2 N(r ) .
The m, sum rule, which is closely related to the po- larizability of the system, can be obtained by solving the constrained Hamiltonian H+A, Q.It can be shown ' that It is easy to check that fp; dr=0, i=1,2.Notice that for l. =0 the transition density p, is not peaked at the surface, having a node at the value of r such that 3p(r)+ rp'(r) is equal to zero.
For an electron functional consisting of a kinetic term (T), a direct Coulomb term (E, ), a Coulomb exchange term of Slater type (Ec, ), and a correlation term of Lang-Kohn type, P dr, b +( 3/4 prr) '   (33) plus a jellium-jellium term (E/ ) and a jellium-electron term, E. , = f V, (r)p(r)dr, (34) a b+(1/a)(3/4m p)' (35) where V. is the jellium Coulomb potential, the scaling of the total energy yields E(a)=a T+aEc"+aEc+ f V/a p(ar)dr This is easy to evaluate only in the monopole case Q = r 2 because it is spherically symmetric.Since we are mainly interested in surface modes, we shall leave the calculation of the polarizabilities for a forthcoming study, except for this I. =0 mode.

III. THE KGHN-SHAM AND THQMAS-FERMI MQDELS
We have seen in the preceding section that to evaluate the RPA sum rules one needs to solve the g.s.KS equations.In the context of metal clusters, this has been done some time ago (see, for example, Refs. 22and 23).The corresponding TF (or density-functional formalism) equa- tions have also been solved (see, for example, Ref. 24).
Consequently, we shall present here the main formulas and results just for the sake of completeness and further reference.
In this section, we use atomic units throughout, i.e. , energies are in units of 2 Ry=27. 2 eV, and lengths in units of the Bohr radius ao =0.529 A.
%"e shall apply the method described in the preceding section to Na and Ag systems.The radius of the jellium sphere is R =r, X', where X is the number of atoms and r, the ionic radius.We have taken r, =4.0 a.u.for Na and 3.0 a.u.for Ag and have used the following LDA functional: Again, the calculation of m3 only involves g.s.quantities, like the kinetic energy T or the electron density p.This is one of the merits of this approach.
Total energies E and rms radii (r )'/2 of some magic Na clusters in the KS and TF ap- proximations, in the latter case for two dilt'erent values of the von Weizsacker coelftcient 13 [Eq.{40}].r) and single-electron potential V(r) (in a.u.) corresponding to the %=20 Na cluster.The thick line is the KS density, whereas the thin one is the variational TF densi- ty.Also displayed are the adjusted TF density (dashed line) and the constant jellium density.The horizontal lines represent some of the single-electron levels.Starting from the bottom, they correspond to the 1s, 1p, 1d, 2s, 1f, and 2p states.
has been adjusted to reproduce the KS results on average.Indeed, P has a big inhuence on the electron surface diffuseness and on the exponential falloff of the electron density as well.Usally, a value of P=1 is taken.' Here we have employed P=0.65 for these two metals, al- though a slightly smaller value would have yielded better results in the case of Ag.To give an idea of how TF(P=1) and TF(f3=0.65)compare with KS, we have collected in Table I the total ene'rgies and rms radii (in atomic units) of some magic Na clusters.
Figure 1 shows the densities corresponding to the %=20 Na cluster, the KS single-particle potential V(r), and the electron single-particle levels, the upper two be- ing empty.The thick solid line is the KS density, whereas the thin solid line is the variational TF density.Also shown is the jellium constant density.
We have adjusted the variational TF densities to a two-parameter Fermi function: taking v = 1.For Na clusters, the average values R =roN' with ro=3.9 a.u.(slightly smaller than the ionic radius r, ) and a= 1.02 a.u.reproduce fairly well the variational TF densities for %~20 clusters.We have also let free the parameter v and found that v=1.The dashed line in Fig. 1 is the adjusted density, Eq. ( 41).The agreement between variational and adjusted TF densities is much better for bigger clusters.
Figure 3 shows the magic %=90 Ag cluster using the same convention as in Fig. 1.Since the A. g ionic radius is smaller than the Na radius, the clusters are more corn- pact for the former than for the latter.For silver, the two-parameter Fermi functions that better adjust the variational TF densities have ro=2.9and a=0.92 a.u.
Starting from the bottom, the single-particle electron lev- els for N=90 Ag are Is, lp, ld, 2s, 1 f, 2p, lg, 2d, lh, 2f, and 3p, the last two being empty.The 3s state, which is next to 1h in N=90 Na, lies at e(3s) = -0.046a.u. in the case of Ag.
To conclude this section we would like to say a word of caution about using sum rules in conjunction with the TF method.The deduction of m & and m3 explicitly uses the KS electron wave functions; only in the case in which we use them to build the g.s.Slater determinant ~P), Eqs.
(4) will give us the RPA m, and I 3 moments.As an ex- ample of the danger of first using a TF LDA approxima- tion and then scaling the electron density according to the first of Eqs. ( 8), notice that if we had proceeded this way, the leading term of the kinetic-energy density, which is proportional to p, would not have contributed to m&(T) for L )1 because it is a pure volume term.
However, if one uses a semiclassical model for p and ãfter arriving at the RPA expressions for m, and m 3, one may expect to obtain a fair estimate of the true RPA mo- ments.That was the point of view of Refs.17 and 18, which we shall also adopt here when we present the TF results.

IV. NUMERICAL RESULTS
Table II coHects the L = 0 results corresponding to three magic Na clusters, N=8, 20, and 92 (the numerical results presented in this section are in atomic units).Notice first that E] and E3 are rather close, indicating that some monopole strength is concentrated in this energy region.A similar concentration was found by Bertsch and Ekardt for the L=1 surface mode.' This table also shows that the main contribution to m3(L =0) comes from m3(j-e).m3(T) and m3(corr) [last term in Eq. ( 36)] contribute very little, altogether =12% for N=8 and = 1% for a cluster of 198 Na atoms.'   Tables III -V show the surface-mode results corre- sponding to the same clusters for several values of L. E3 decreases at very high values of L (see especially Table III).This is due to the finite potential well at which elec- trons are submitted, as is discussed in Ref. 25.It is quite interesting to observe that m3(j-e) and m3(e-e ) have op- posite signs and that for a given ¹ they tend to cancel each other when L increases.Indeed, for L )&1 the field r YIO is essentially acting on the external part of the electron density.For this part, electron screening almost cancels out the jellium Coulomb potential, and, conse- quently, the Coulomb contribution to m 3.
The following physical picture thus emerges.For a given cluster, the excitation energy corresponding to small values of L is mainly determined by the Coulomb energy, the electron kinetic energy playing a minor role.
When L increases, the Coulomb contribution becomes less important due to electron screening.Above a certain value of L, the kinetic energy is the only appreciable con- p=p, =ne(Rr) .
Then, for L =0, Eqs. ( 29 For L . ) 1, .usingEq. ( 42) and = -, '(3ir ) n:an ~, we get from Eqs. (1-6), ( 19), (26), and (28) tribution to the surface-mode restoring force.At larger L, one expects the mode to lose its collective character: the r YIO field will only probe the outermost single- electron contribution.Putting it in a different way, the main contribution to m3 will come from the most exter- nal electron shell and the excitation will be of electron- hole type.Consequently, there is a critica1 angular momentum L" for the existence of collective surface modes in a given cluster ¹.
Figure 4 shows the energies E3 versus ¹ for Na magic clusters of ¹=8,20, 40, 58, 92, 138, and 198  respectively.We want to stress again that E3 is an upper bound to the resonance peak energy.
The solid lines are the TF results obtained using the kinetic-energy density (40).TF calculations are straight- forward and actually we have performed them up to ¹=500.Notice that the TF results nicely average the KS ones.The main discrepancy is around ¹=50and it is due to a shell effect.Indeed, the ¹=58cluster is ob- tained by adding 18 1g electrons to the ¹ 40 cluster.
These 1g electrons in the more external shell contribute appreciably to the high multipolarity excitations.E3(N) has a quite different behavior depending on the value of the angular momentum.For high L the energies decrease with ¹ and tend to some asymptotic value.For small L the energies are rather independent of ¹, slightly increasing with it for small clusters.From the results shown in Tables III -V, we infer that these trends could be related to the dominance of m3(T) at high L and of m 3( C) at small L (actually, for L = 1 this is the only con- tribution to m 3).

V. DISCUSSION AND CONCLUDING REMARKS
In this paper we have studied the surface modes of clusters and spheres of metal atoms within the RPA sum-rule approach.We have obtained with RI'A pre- cision the m& and m3 moments of the strength function for operators of the kind r Yl o.This allows us to define an average energy which, neglecting electron kinetic- energy contributions and surface-diffuseness corrections, computed the L =0 and L =1 modes for Na clusters.
Figure 5 shows the results we have obtained using different values of a.For each multipolarity, the upper curve corresponds to a = 0.72, the middle curve to a=1.02 (fit value), and the lower curve to a=1.32 a.u.
%'e see from this figure that the bending down increases when a increases.Besides the bending, increasing the surface diffuseness also causes a global shift of the E3 curves downwards.Thus, we conclude that the red shift- ing of the dipole mode from the classical Mie value is a surface-diffuseness effect.
We have also verified that the Na adjusted density (41) reproduces the TF results fairly well up to L =4.Above this value, it fails for clusters N 40, but still yields good results for bigger L, provided N & 40.
Figure 6 collects the E3 energies for silver.The trian- gles represent the magic Ag clusters N= 8, 20, 40, 58, and 90 results.As in Fig. 4, the dashed lines have been drawn only to guide the eye.The solid lines are the TF energies.
From bottom to top, the curves and points correspond to L = 1-7 and to L =0.We see that in this case the TF re- sults do not average the KS ones.Had we used for Ag a smaller value of P (=0.5), the agreement between TF and KS would have been much better.
Qualitatively, Figs. 4 and 6 look very much the same.
In the case of Ag the shell effect from N=40 to 58 is also quite distinct.The bending down of the L=1 mode is more pronounced for Ag than Na.This is so because the diffuseness corrections to the bulk E3, Eq. (46}, are of or- der (a lroN'~) and consequently bigger for Ag than for Na.The density (41) with the parameters corresponding to Ag can be used for a quick estimate of E3, with the same remarks as in the Na case.Finally, we can use Eq. ( 46) to estimate a lower bound of the critical angular momentum L"discussed at the be- ginning of this section.%'e have argued that the mode will not lose its collective character until the kinetic dom- inates over the Coulomb contribution to m3.Imposing that both contributions be the same, For high I values, it yields reduces to the classical Mie formula when N goes to infinity.
The expressions we have deduced are very simple to use and could be applied to other metals besides the ones we have studied here.The key ingredients for computing m, and m3 are the ground sta-te kinetic and particle elec- tron densities which can be easily obtained from a spheri- cal Kohn-Sham calculation.In this way, we are intro- ducing a realistic electron-density profile.The importance of surface inhornogeneities on the dispersion of sur- face modes, i.e. , the dependence of a surface-mode fre- quency upon the radius and surface diffuseness of the sphere, was stressed in Ref. 3 (see also Ref. 5}.This dispersion comes out from our model in a very natural way.For L & 1 we have an R (N) dependence from quantum-size effects associated with the electron ki- netic energy [Eq. ( 46)]; the surface-diffuseness effect is far less important, at least for small clusters and angular mo- menta.On the contrary, electron surface diffuseness is the only dispersion source for the dipole mode.
We have found that shell effects are important only in cases where high-angular-momentum electron orbits are involved.Thus, a semiclassical TF calculation gives re- sults that agree well with the KS ones, provided an effective Weizsacker correction to the TF kinetic-energy density is included.
The present method is of great simplicity and useful- ness.It can be readily used to obtain E3 collective ener- gies with RPA precision, or just the average trend by us- ing a semiclassica1 model like TF.Moreover, the expres- sion for E3(L) is simple enough to allow for a qualitative analysis of it as a function of the number of atoms ¹ For small values of L, we have stressed the important role played by the electron-density diffuseness to red shift the collective energies from the classical Mie value.The role played by the electron-density surface diffuseness on the dynamical polarizability of surface modes was already pointed out by Ekardt.'   For L ~3 the increasing contribution of the electron kinetic energy to E3 completely changes the picture and we may have a sizable blue shifting of the collective ener- gy for N ~100.It is worth mentioning that recent exper- iments performed on Ag spheres of 40 -70 nm have been reported, confirming the validity of the Mie model to study L = 1 and 2 optical-absorption spectra.Our analysis shows that no deviations from the Mie model should be expected for these Ag spheres and L modes.In order to observe a significant deviation from this model, one should appreciably excite L «3 modes on much smaller spheres.However, for such spheres retardation effects might be large.
The interplay between Coulomb and kinetic contribu- tions to E3 can be used to determine roughly at which angular momentum the modes are no longer collective.
We have argued that this could only happen after the ki- netic contribution takes over the Coulomb contribution.
Using this simple argument, we have found a value for L" that agrees well with previous estimates by Ekardt.Finally, we have also discussed the L =0 volume mode calculating with RPA precision the m &, m&, and m3 moments of the strength function.We have shown that there is some collectivity concentrated in this mode, which lies at much higher energy than the surface modes which we have thoroughly studied.
During the completion of this work, we became aware of a calculation by Brack, who arrived at the same ex- pressions for m3 using a similar technique.However, for the numerical applications he used a trial density TF method instead of the quantal KS and only discussed Na clusters.In this sense, we consider both calculations as complementary.
5(r -r;)V;, l The m3 sum rule will consist of three terms, one [m3(T)] coming from the kinetic-energy density in H and two from the electron-electron Coulomb energy m3(T)= [m 3(e-e )] and the jellium-electron Coulomb energy expansion of ~rl -r2~' , the first, in tegral in Eq. (21) reduces to f Pl I Pl 2 d d (r ) (r ) rp r FIG.1.Densities p(r) and single-electron potential V(r) (in a.u.) corresponding to the %=20 Na cluster.The thick line is the KS density, whereas the thin one is the variational TF densi- ty.Also displayed are the adjusted TF density (dashed line) and FIG.2.TF densities of some Na clusters.From inside out- wards, they correspond to X= 8, 58, 198, and 300.
FIG 4 The sum rules mi and m3 for Q=r YIO , one may estimate the centroid and vari- ance of the strength function by evaluating the three RPA moments m &, m &, and m 3.If most of the strength is in a narrow energy region, as it is in the case of some resonance states, then E, and E3 are estimates of the Consequentlymean resonance energy.Conversely, if E, and E3 are close, we deduce that some strength is concentrated around these values.For m & and m3 we have , atoms.Starting from the bottom, the curves and points correspond to L =1-7 and to L=O (upper curve).The L= 1 -7 KS results for these clusters are represented by triangles, whereas the L =0 KS results are represented by crosses.The dashed lines linking the KS points are only to guide the eye.The crosses lying near the L=1 curve at ¹=20,92, and 198 are the results obtained in Ref. 14, and the two crosses near the left bottom corner are the ¹=8and 20 experimental L=1 values of Refs. 9 and 8,

TABLE II .
RPA L =0 volume mode energies E& and E3 as well as m &, m &, and m3 sum rules {in a.u.) for several Na clusters.The different contributions to m 3 are also displayed.

TABLE IV .
Same as Table III for N=20.