Chapter 4: A Deformed Sphere Model of the Electronic Structure of Small Metallic Particles

Chapter 4

A Deformed Sphere Model
of the Electronic Structure of Small Metallic Particles

Introduction
A Spherical Approximation: The Shell Model
   Theoretical Description
    Results
A Perturbed Sphere and the Loss of Symmetry
    Ideas from Previous Models
    Difficulties with Traditional Models
    Theoretical Description
    Results
Deformations and Local Minima
Computational Method

Introduction

The availability of ionization potential (IP) measurements of silver clusters offers a unique opportunity for theoretical investigations that might explain the complex structure observed over the size range of 1-100. Certainly the patterns in the IP's must be manifestations of an underlying physical phenomenon. However it is not obvious whether the phenomenon is chaotic in nature – implying features that are difficult to predict and sensitively dependent on the characteristics of individual clusters – or whether it is the result of a few basic forces whose interaction is predictable but nonetheless result in a non-smooth behavior.

Fortunately, the complexity of the patterns suggests that if a model where to predict this behavior it would certainly confirm our understanding of such systems. The question is what approximations, if any, can be made to model the phenomenon and yet still explain the detailed structure observed. For the smallest clusters, Ag₂-Ag₉, the ionization potential pattern represented in this series has been investigated at the ab initio pseudopotential level, but, as is well known, this model is not feasible for the simulation of much larger particles at present. Although rapid improvements are being made with pseudo-potentials,⁵⁷ the calculation of electronic structure is still difficult for transition metals. Possibly more problematic are the very large number of geometric isomers that any structural model must deal with in some reasonable fashion. Density functional methods have extended the reach to larger particles and the dynamic Car-Parinello method has to some extent addressed the problem of large numbers of geometric isomers,⁵⁸ but a consistent study of clusters up to 100 atoms is not presently feasible.

At the other extreme, there are simple "shell" models that treat metallic clusters as spherical containers of electron density within an effective potential. The delocalized nature of the valence electrons in metallic clusters has suggested approximations that treat the "valence" or atomic s1 electrons independent of the "core" electrons.⁵⁹ This model arose historically from a similar treatment of particles in the nucleus, and has been applied to metal clusters since 1981.^{59 60} It uses a background potential to simulate the ion core and electron interaction potentials. The model owes its success to the surprising postulate that both the electron-ion core interaction and the rather strong electron-electron interaction can be approximated by an effective potential that is independent of the explicit electron-electron or electron-ion distances. Some forms of this model, the 'jellium' models, treat the electron interaction potential within the local density framework.⁶⁰ Since the core electrons are ignored in both the shell and jellium models, they are expected to be applicable to metals such as the alkali metals that have a large energy gap between valence and core electrons. However for the s1 transition metals, Cu, Ag and Au, the gap is smaller – small enough in the bulk that the underlying d-electrons contribute to conduction. In this situation the validity of the jellium and shell models is not so clear, but has still been successful at least for ground state properties, as will be shown.

In an attempt to explain, using the simplest possible model, the patterns exhibited in the experimentally observed ionization potentials, we begin with a comparison to the shell model. Then, after showing that this is not sufficiently accurate and suspecting that the IP patterns indicate geometric structure that is ignored in the sphere, we investigate distortions to the spherical model. These suspicions are based on previous work that has suggested the importance of distortions to the sphere.^{62 62} We have not worked within the jellium model to treat the electron interaction potential more explicitly as has been done elsewhere.⁶³ Such a treatment is not expected to significantly change the features observed. In this sense, we follow in a recent tradition of similar ideas, ⁶⁴ but with a refined treatment of the model. A successful comparison of this new model with experimental results will show that "the most vexing question concerning metal clusters" –⁶⁵ their structure – is partially answered in the process.

A Spherical Approximation: The Shell Model
Theoretical Description

In the shell or jellium model an electron is assumed to move in a spherically symmetric effective potential caused by two components: the ion-core and the electron-electron interaction. The potential can be assumed to have some particular form or, alternatively, be separated into its components and then evaluated within the local density approximation. The restriction to self consistency can also be required of the electron density . This is commonly referred to as the self consistent field (SCF) jellium model. In either case, the single particle energy levels are found as solutions to the Schrödinger equation within spherical polar coordinates,

(4.1)

As in the case of the hydrogen atom, the equations are separable into radial and angular parts. The wave function simplifies as a product of two functions,

(4.2)

Angular functions take the form of the well-known spherical harmonics with eigenvalues l, m. Expression of the Laplacian operator Ñ 2 in spherical polar coordinates leads to another eigenvalue equation for the radial functions,⁶⁶

(4.3)

Analytical solutions to this equation within arbitrary radial potentials are not feasible, so that solutions, Rnl(r) and e nl, are usually formed numerically. Solution of this eigenvalue ordinary-differential equation forms the core routine in programs used to compute the energy levels of spherical jellium clusters.

Orbitals and their energies are then distinguished by quantum numbers n and l. (The convention for n,l is based on the nuclear convention with n=1,2,etc. rather than n=l+1,l+2,etc. as in the Coulomb potential.) Large differences in the energies of neighboring clusters are assigned as 'shell-closings'. In the independent electron model these correspond directly to gaps in the eigenvalue spectrum of the Hamiltonian. Within most potentials, closed-shell stable species occur for 1—100 electron clusters as follows:

N	Configuration
2	1s2
8	[2] 1p6
18	[8] 1d10
20	[18] 2s2
34	[20] 1f14
40	[34] 2p6
58	[40] 1g18
92	[58] 2d103s21h22

The shell model has historically arisen as a model of atomic nuclei dating back to 1955. A great deal is known about its properties for nuclear wave functions.⁶⁷ (The filling order follows the pattern: SPuDS iF PuG DiSH oF PIG.⁶⁸) However the correspondence between cluster and nuclear wave functions is complicated by the importance of spin-orbit coupling in the nucleus and its relative unimportance in metallic clusters.^{69 70}

Results

Here we compare our experimental ionization potential determinations to those of the spherical shell model. The effective potential used is the Woods-Saxon potential, which also originated in nuclear physics but has been used extensively for metal clusters,⁶⁰

(4.4)

This is known to be a good approximation to radially averaged potentials from self-consistent local density calculations.⁷⁰ The radius of the cluster is R0=rs N1/3 with N the number of atoms and rs (3.02a0 for silver) the radius of a sphere that contains one electron in the bulk.⁴⁶ The depth of the potential U0 is the sum of the Fermi energy and the bulk polycrystalline work function, 4.26eV⁴⁶ and 5.49eV⁴⁹ respectively for silver. The variable m determines the steepness of the potential edge and is the sole parameter of the model. We chose a reasonable value m=rs/4 and did not vary it to match the experimental results. A similar value, 0.38rs, has been employed in the past for sodium.⁶⁰ Future research might investigate the effects of this term.

The single particle energies are obtained by solving the Schrödinger equation in the effective potential. To compute total energies, the orbitals are filled sequentially and the total energy is approximated as a simple sum of single particle energies. This includes the effective electron-ion core potential correctly but counts the electron interaction potential twice.

In the following discussion the vertical IP is assumed to be relevant for experimental photoionization. In the independent particle model this is simply the single particle energy.

Figure 4.1. A comparison of experimentally determined IP's() and those predicted by the spherical Woods-Saxon model() for silver clusters. Numbers indicate the predicted shell closings. Labels in italics indicate the shell filling order. Most of the shell closings are observed experimentally although the resemblance is weak otherwise.<full image>

Figure 4.1 shows a comparison of the experimentally determined IP with the results of the spherical Woods-Saxon shell model. The relatively high IP of some clusters indicates closed shell species. For the most part similar steps are observed in the experimental results: six shell closings are observed while three are weak or not observed. Clearly this suggests some validity to the model. The gradual increase in IP within a shell is a result of the increasing radius of the cluster. For an individual cluster these states are degenerate. Note however that the even-odd oscillation and a number of large steps, for example at 48, are not predicted. The failure to describe the details over this range suggests that either the presupposition of the effective potential in the shell model is incorrect or that these degeneracies are broken by the geometry of the cluster.

A Perturbed Sphere and the Loss of Symmetry
Ideas from Previous Models

Clearly the spherical shell model is a poor predictor of the ionization potential for those species that are not closed-shell. A number of refinements to the model have been proposed that consider non-spherical potentials. Presented here is a thorough investigation of the properties of a refined ellipsoidal model that avoids many of the approximations of other models. As a result, much stronger confirmation of the accuracy and validity of the model is obtained than has been reported previously. The model is shown to be accurate, not on an absolute energy scale, but in predicting the local patterns observed. In contrast to other ellipsoidal models there are no parameters that must be adjusted for each cluster size. This facilitates a realistic comparison to experimental findings. Since the original goal was to describe the patterns observed in the ionization potentials of silver clusters, we emphasize comparison with those experiments.

The large number of orbital degeneracies in the shell model is mathematically and computationally a pleasant feature; however removal of this orbital degeneracy through atomic motion can reduce the energy of the system through a Jahn-Teller distortion.⁷¹ (In any non-linear system there exists some vibrational mode that removes the degeneracy of an orbitally degenerate state.) Since these clusters are most likely close-packed it is probable that the breaking of symmetry would correspond to spheriodal or (even further) to ellipsoidal distortions. (A spheroid implies two equivalent axes and an ellipsoid no degenerate axes.) Such a distortion produces hybrid orbitals as mixtures of spherical orbitals, however in this model the extent of mixing is dependent upon the electronic configuration that forces the distortion.

The energy of such a system as a function of the distortion along each axis defines a potential energy surface(PES). Since volume conservation reduces the degrees of freedom to two, the surface is easily visualized. A computationally simpler model with a similar PES features independent electrons in a three dimensional box, again under volume conservation.

Figure 4.2. Two examples of axis distortions under volume conservation. Ellipses represent Jahn-Teller type distortions that occur when one p-orbital is occupied. The wave functions that result are hybrids of s and p orbitals. Elongation along the x-axis is limited by increased energy of the occupied s-like orbital. An analogous phenomenon occurs if the volume is a box containing electrons. Shown are the quantum numbers of the highest occupied orbitals nx,ny,nz and the approximate distortions that occur to minimize the total energy. The first is 'closed shell' with 16 electrons and the second with 17 or 18 electrons would, if it existed, have a low IP.<full image>

Figure 4.3 <full image>.⁷²

Analogies can also be made between this and Crystal Field Theory, regarding the breaking of symmetry of an ion by the electric field of its crystalline environment. However a group theoretical analysis of the splittings of the orbitals, here s through h, in an ellipsoid is unfortunately not very informative. The complete ellipsoidal distortion with point group D2h involves three parameters and separates the lower orbitals into their constituent suborbitals, e.g., p® px+py+pz and d® 4 terms. Splitting in a pure prolate potential is shown in Figure 4.3. This is the difficulty of the problem – its low symmetry.

Spheroidal distortions have cylindrical symmetry and are either prolate or oblate. A shell model for metallic clusters with such distortions was proposed by Clemenger in 1985 ^{61 62} based on a similar model for nuclei investigated by Nilsson. Experimental research on nuclear states had suggested the importance of quadrupolar distortions. In both cases the radial potentials are harmonic oscillators to which an additional empirical term is added to lower the energy of high angular momentum states and effectively reduce the steepness of the parabolic potential. Nilsson had computed the energy levels in this potential, however the nuclear model includes spin-orbit coupling that is not significant here. The cluster model is described in a recent review by de Heer.⁶⁴ The spherical Hamiltonian for an electron in this model is (in atomic units):

(4.5)

The first two terms establish a three dimensional harmonic oscillator of frequency w 0 and the final term is an empirical adjustment dependent on the parameter U and the quantum numbers l and n. Distortions are permitted by allowing two or three independent oscillator frequencies, w x, w y, w z. Under this distortion l and n are no longer good quantum numbers but approximate quantities.

Most implementations of the Nilsson-Clemenger model are restricted to either the ellipsoidal model without the anharmonic term⁷⁴ or the spherical model with the third anharmonic term. For small clusters it is common to ignore the third term completely. Inclusion of the anharmonic term in a spheriodal potential significantly complicates the analysis, coupling once-independent variables. Even so, the model has been extended by Saunders to include both the anharmonic term and ellipsoidal, i.e., non-axially symmetric, distortions.⁷⁵

Difficulties with Traditional Models

Although it should not be understated that the Nilsson-Clemenger model plays an important role in explaining the electronic structure patterns for small (N<25) metallic clusters, its usefulness for larger clusters is restricted. The anharmonic term U in both the Nilsson-Clemenger and in the Saunders Hamiltonian is shell dependent. Its parameters are adjusted to obtain the best match to experimental results or self-consistent jellium potentials. In nuclear physics its parameters have also been fit to spheriodal potentials.⁷⁶ The model is typically not extended beyond 40 atoms due to the necessity of varying U or the form of the anharmonic term to obtain a reasonably shaped well. Without extensive parameterization the harmonic oscillator potential fails to model the flat-bottomed potential of larger clusters.

An additional complication of the Nilsson-Clemenger model regards the conservation of cluster volume. This is accomplished by requiring: w 0=w xw xw x. This becomes a poor approximation for larger clusters since it ignores the l dependent portion of the potential that is important for high angular momentum orbitals and large clusters.⁷⁷

Theoretical Description

We present here an original model that employs an effective potential more realistic than the harmonic oscillator based models and yet allows for complete ellipsoidal distortions. Such a model is required by the larger clusters observed in the previous experiments. It should be emphasized that the model contains no parameters that are fit to experimental data. In addition, the model employs a numerical integration that could easily handle other types of radial potentials and, with some modification, potentials of lower symmetry simulating the ion cores. The disadvantage is that it relies more heavily on computational effort than the other models and thus lacks the analytical nature that encourages insights not so easily discovered in a numerical scheme. However in the future it may be of interest to introduce atomic cores as perturbations on the effective potentials. This would not significantly complicate the numerical model.

As cluster size increases, the potential begins to represent a flat-bottomed sphere or box more than a harmonic oscillator. Thus for this model the harmonic oscillator is not used as a starting point. Single particle energies are obtained by solving the Schrödinger equation and the total energy is approximated as the sum of single particle energies, as is done in the shell model. The Hamiltonian is expressed in Cartesian coordinates. For each electron this is simply,

(4.6)

and is defined on a basis to give a matrix eigenvalue equation,

(4.7)

Eigenvalues of the Hamiltonian are determined by projecting onto a basis set and diagonalizing. The Woods-Saxon potential is used for the effective radial potential. To effect the distortion to an ellipse of dimensions Rx,Ry,Rz, the potential is evaluated in scaled coordinates r' where x' = (R0/Rx)x and similarly for y' and z'. The ellipsoid dimensions, which approximate the cluster dimensions, are subject to the constraint of volume conservation, R03=RxRyRz. Although the model requires a numerical integration over the cluster, this is performed rapidly with a suitable basis set . A more detailed description of the computational method is given in the next section.

Note that the two parameters that are present in the Nilsson-Clemenger model, the anharmonic correction and the oscillator frequency, are absent in this model.

Calculations have been done on self consistent field local density (SCF-LD) jellium under the restrictions of spheroidal distortions, referred to as a 'Sommerfeld droplet' model.^{78 79} Here we relax the requirement of self consistency, use an effective potential, and allow complete ellipsoidal distortions. This results in a comparatively simple model.

Results

Thus the interaction between the ion-cores is ignored – except to conserve volume – and the combined effects of the electron-ion core potential and the many-body electron-electron potential have been approximated with a 'tapered' ellipsoid that distorts to reduce the total energy. Figure 4.4 shows that ellipsoidal distortions add what appear to be random noise to the IP pattern. The high degeneracy (2l+1) of the symmetric model is destroyed.

Figure 4.4. The effect of allowing ellipsoidal distortions to the spherical model. The sawtooth pattern () represents the ionization potentials of the spherical model and the curve (—) those of the distorted model. Closed shell properties at N are preserved while N+1 change significantly. <full image>

Figure 4.5 compares results for the smallest clusters studied. Although there is some discrepancy for the smallest clusters, N<9, the model imitates experimental results remarkably well. The prominent even-odd oscillations are also present. This suggests that the effect is due to a greater stabilization in the highest occupied orbital (HOMO) of the even electron cluster in comparison to the odd electron cluster. The doubly occupied HOMO has a stronger drive to distort the cluster and stabilize itself at the expense of lower energy orbitals. To a large extent electron spin is unimportant in this phenomenon. These are some of the first results on clusters of this size showing that even–odd effects are predicted reasonably well by the model .

Figure 4.5. Comparison of the distorted model () with experimental values () for the smallest silver clusters studied. Poor correspondence with the smallest clusters is a consequence of the independent electron model and a disregard for image charge effects as discussed in the text. The details of the pattern are, however, reproduced surprisingly well, particularly the odd-even oscillations and local features. <full image>

For clusters larger than 12, there exist local minima on the potential energy surface (PES) as a result of the flatness of the potential inside the cluster. It is common to find a number of stable structures for each size cluster. The IP difference between local and global minimum energy isomers is not dramatic but, as discussed later, is noticeable. In this analysis, tails on the photoionization efficiency curves (PIE) are attributed to the lower IP isomers and highest IP isomers are compared with experimental results. Note that this comparison is not always to lowest energy isomers. This is done because in all cases except Ag₂₁ the agreement with high IP isomers is better than with low IP or lowest total energy isomers. The following section discusses this further.

Over the size range of 6-40 atoms, Saunders' similar ellipsoidal model⁷⁵ has been applied to copper cluster electron affinity measurements with similar success up to the size of Cu40-.²¹

Figure 4.6. An extension of Figure 4.5 comparing model () and experimental () IP's for somewhat larger silver clusters. The absolute value is now more accurate and, with the exception of 49 and 50, the prediction is very accurate. A local transition from ellipsoid to prolate occurs at 49 that may explain this discrepancy if the real cluster does not make the transition at the same size that the model does. <full image>

Figure 4.7. A continuation of Figures 4.5 and 4.6 for the largest clusters studied showing model() and experimental() ionization potentials. The correspondence is still significant although there is less structure in the experimental values. This may be due to the difficulty of resolving large non-isotopically pure silver clusters through mass spectroscopy. The step at 74 is not predicted although 74 is predicted to have a high IP. This is a result of a structural transition that occurs in this region. <full image>

Figures 4.6 and 4.7 complete the set of comparisons. The correspondence between experiment and theory is exceptionally good for the larger clusters. In some cases even detailed features are reproduced. Beyond approximately Ag₉ the ionization potentials of silver clusters are described very well by the model. The clusters behave much like electrons bound to positively charged droplets that distort to reduce the electronic energy. We conclude, as other researchers have also concluded,⁷⁸ that the spherical jellium or shell model does not fail because it ignores the detailed ionic structure, but that it is very successful if the spherical assumption is relaxed. Certainly these results are strong confirmation that the electronic structure governs the cluster geometry even for clusters larger than 40. It appears that for ground state properties and clusters of this size the ellipsoidal distortions play a more significant role than the details of the ionic lattice.

There are a number of possible reasons for the relatively poor predictions for small clusters. As has been suggested by M. Kappes, this model does not account for the work required to remove the electron against the electrostatic force of the image charge in the cluster. The work required to remove an electron from a neutral metal sphere is a function of the cluster radius R,^{44 80 81}

(4.8)

Figure 4.8 shows the effect of adding this term to the calculated ionization potentials. Based on a comparison of the global behavior of the experimental and theoretical results, this appears to be the necessary correction. However, the correction leaves an approximately constant error of 1.5eV. This probably indicates an incorrect depth to the Woods-Saxon effective potential well, Equation 4.4. Naturally, the potential could be adjusted to match the results, but in avoidance of arbitrary modifications this has not been done.

Figure 4.8. Effect of adding an additional energy of e2/(2R) to the calculated IP's from the model. Both experimental () and calculated () photoionization potentials are shown. This term represents the work required to remove an electron from a positive charge enclosed in a metal sphere of radius R. It corrects the global trend in the theoretical results so that they closely imitate the experimental results although with a constant error of about 1.5eV. For all but the smallest n<7 clusters the model is successful in explaining the observed behavior. <full image>

There are additional reasons that could explain discrepancies between the model and experiment. The first is related to

the effective potential assumption in which IP's are equated with the single particle energy levels. Since ionization significantly changes the electron density of a small cluster, the effective potential within the small ion differs from that of the neutral. Hence the effective potential approximation breaks down and the IP no longer corresponds to single particle energies. The approximation is better suited to larger metal clusters where the relative change in electron density is less. Second, there may be an increase in the involvement of d-electrons that occurs over this size range. Third, the spherical or ellipsoidal shape approximation is rather crude for such small clusters and the underlying atomic positions could perturb the orbitals significantly.

Finally, another possibility is an incorrect estimate for the effective cluster size. From local density calculations it is known that the electronic density extends beyond the cluster radius delimited by rsN1/3. The absolute ionization potential is a sensitive function of the cluster radius. For example if, as is suggested by Perdew,⁴⁴ an additional electron spill-out distance of a=1.34a0 for silver is added to the cluster radius, so that R = a + rsN1/3, then the absolute IP's are closer to experiment — see Figure 4.9. The effect of this term diminishes as the radius increases. Since such a term is somewhat ad hoc and not clearly independent of cluster size, it was not included in other calculations.

Figure 4.9. The effect of adding an electron 'spill-out' term of 1.34a0 to the cluster radius. IP's are shown for the experimental results() and predictions with () and without () electron 'spill-out'. <full image>

We discuss first a number of detailed comparisons and then investigate the causes of the patterns in the ellipsoidal model.

There are three locations of partial failure in the model: Ag₂₁, Ag₃₁, Ag₄₉. Interestingly enough, these are the clusters with the lowest experimental IP's. They have uniquely long PIE tails that continue into the red end of the spectrum (Ag₂₁ PIE was shown in Figure 4.3). This supports the conjecture that there are a number of isomers present in the molecular beam and that the PIE linearization procedure does not differentiate between them. For example, calculations on Ag₂₁ suggest two stable ellipsoidal geometries, one of which has a low IP. The other isomer is a local minimum (meta-stable) and possesses a high IP. Theoretical clusters also make a distinct transition from triaxial to prolate at Ag₄₈–Ag₄₉. Experiments suggest that Ag₄₉ is still triaxial or that there are two isomers present, one of each configuration. (Alternatively, the ionization potentials of sodium clusters do exhibit a locally high IP at Na50 (See Figure 5 in the previous chapter and Ref. .) suggesting that they do follow the model and make triaxial-prolate transitions as the model suggests. The exact location of the transition is probably sensitive to the potential form.) As in the case of Ag₁₄ , these 'shape' transitions are most important in the vicinity of half way between closed-shell species, e.g. Ag₈–Ag ₁₈ and Ag₃₄ -Ag₅₈.

Of the largest clusters studied, two of them, Ag64 and Ag₇₄ are of particular interest. Interest was sparked in Ag₆₄ since it is the only even numbered cluster predicted to have an IP lower than its neighbors. This is not observed through linearization of the PIE curves; however, a subsequent examination of the mass spectrum suggests that some isomers of Ag₆₄ do indeed have an IP lower than their neighbors. See Figure 4.10. The cause is presently unknown.

Figure 4.10. Two mass spectra of silver clusters. Labels indicate the number of atoms in the cluster. One spectrum taken at high photon energies ionizes most of the clusters while the other at the low photon energies, and near the signal to noise limit, shows that Ag₆₄ may have isomers with a low IP as was found theoretically. <full image>

Ag₇₄ had been noted to be unique by its experimental IP in one of our previous publications.³⁸ The existence of a high IP for Ag₇₄ supports the validity of this model and the importance of triaxial distortions. In this size region halfway between the major shell closings at 58 and 92 there are a number of stable isomers with small energy differences. The high IP of Ag₇₄ is an indication of this transition. In this series Ag₆₉ is the last oblate structure, and the transition to prolate is completed at Ag₇₃. Ag₇₄ is then the first even electron species in the series to be prolate – stabilizing the HOMO and raising the IP. The experimental results suggest that this transition is more abrupt but nonetheless present. This may be the first confirmation of a simplification of this phenomenon recently proposed to exist in 'super-shell' sized clusters.⁸²

Deformations and Local Minima

There are at least two possible explanations for the surprising agreement between such a model and the experimental results. First, it is possible that the effective potential seen by the electrons is very much ellipsoidal. This would be due to the comparatively small size of the ion cores relative to the de Broglie wavelength of the valence electrons. A second possibility is that there are many isomers present in the molecular beam and that the 'average' ionization potential, taken over all geometric isomers, is equal to that of an ellipsoidal cluster. An additional possibility is that the surface of the cluster does not distort to resemble an ellipsoid, but that the ion cores within the cluster rearrange to reduce the potential energy of the partially filled shell with a similar effect. For small clusters, whose atoms all lie on the surface, this is the same phenomenon. For larger clusters these are distinctly different.

A look at the shapes of these model clusters gives some indication of the behavior that causes such erratic patterns. Figures 4.11 through 4.13 plot the ellipsoid dimensions as a function of the number of atoms in the cluster. These plots indicate the increase in axis length over that of the spherical cluster for each axis – d x,d y,d z – where Rx=d x+R0 with R0 the radius of the spherical cluster. This is drawn for one axis in the diagram to the left.

In most cases the spherical symmetry of standard jellium closed shells is still observed. Clusters with sizes N± 1, where N is the standard closed shell species, are generally prolate. Sizes N± 2 are generally oblate. For larger quantum numbers the connection to a spherical model begins to break down and this is no longer true.

Figure 4.11. Ellipsoidal distortions shown as increases in axis size for each of the three dimensions. Plotted is d versus N where axes have length R=d +R0, with R0 the radius of the spherical cluster and N the number of atoms. The three symbol types (,,) represent the three axes – x,y,z. The assignment of x,y, and z is arbitrary. In the range 2—8 the transition spherical-prolate-distorted-oblate-spherical occurs as expected for a p-orbital. Other patterns are not so clear but the same trend is observed for the filling of each shell. 'Shape' isomers often exist near major transitions, e.g., Ag₁₄ or Ag₁₅. Note that rs is 3.02a0 for silver. These are the shapes with a high ionization potential that have been compared to the experimental results and not necessarily minimum energy shapes. This is why Ag₁₉+ for example is spherical. <full image>

Figure 4.12. Similar to Figure 4.11 for the next range of clusters. Again the filling of the 2p-shell produces the expected distortions. A partially filled high angular momentum g-shell induces large distortions. Effects of the ellipsoidal to prolate transition at Ag₄₈-Ag₄₉ are observed experimentally. This is a transitional region where 'shape' isomers exist. <full image>

The absolute axis change remains relatively constant at <3a0 or <6Å, while the relative distortion ranges from 46% for Ag₄, down to 11% for Ag₁₀₀. This is still a rather large perturbation and explains the failure of the spherical model to account for features in this size range.

Figure 4.13. Similar to Figures 4.11 and 4.12 for the largest clusters studied. Significant distortions make the 2d, 3s, and 1h orbitals poor descriptions of the molecular orbitals but the 10-2-22 filling pattern is still discernible. <full image>

It should be noted that for the largest clusters studied there are some uncertainties in the precise value of the distortion due to approximations in the calculations. These are effects difficult to estimate and so no error bars are included in the figures.

Figure 4.14 contains a comparison with two other models over the range of available data. The comparison is only approximate since the other calculations are for Na clusters. For small clusters with N<12, similar patterns are found to the ellipsoidal harmonic oscillator model, but beyond Ag₃₄ this is no longer true. The results of the Nilsson-Clemenger Hamiltonian do not show a spherical shape to Ag₃₄. Note the transition between the ellipsoidal models that occurs at Ag₁₄. The shape of Ag₁₄ is known to be sensitive to the details of the model.⁸³ In our results two stable configurations are observed. The minimum energy configuration is prolate, although the two species have similar ionization potentials. These two stable isomers were also found in more sophisticated SCF local density jellium calculations by Brack et. al. ⁸³ The same workers also found Ag₁₂ to be triaxial, and this has been confirmed experimentally by three peaks in the surface plasmon.

Figure 4.14. Comparison of ellipsoidal distortions of this work for Ag clusters with ellipsoidal harmonic oscillator of K. Selby et. al. for Na clusters⁷⁴ and a more detailed tight-binding Hückel model of Poteau et. al. also for Na clusters.⁸⁴ Note the different y-scales since Na rs=3.93a0 and Ag rs=3.02a0. <full image>

For model clusters with ion cores shown in the third frame, the magnitude of the distortions depends on how one defines the edge of the cluster. Here the distortion has been defined in terms of the axes of an ellipsoid having the same moments of inertia as that of the cluster. Cluster geometries derived from calculations that explicitly treat the ion cores are in general not axially symmetric.

-----------------------------

While trying to perform geometry optimization, the author was surprised to find clusters trapped in local minima on the potential energy surface(PES). Naturally this complicated the problem. An investigation of the PES of the box-of-electrons model, previously mentioned, revealed similar features. For many clusters there are a number of stable configurations. A subsequent review of the literature revealed that spheroidal models occasionally yield stable prolate and oblate species.⁸⁵ Penzar and Ekardt have observed two isomers for Na₁₃ and Na₃₀ in an SCF spheroidal jellium model.⁷⁸ Some very recent work on large clusters (N<3000) noted local minima for the infinite square well potential.⁸⁶ In our experience all the configurations or local minima observed are located near the spherical shape in phase space, i.e., within ± 4a0. Figure 4.15 plots the simple PES for an eight electron cluster. This should be compared to the intricate surface discovered for Ag₁₅ in Figure 4.16.

Figure 4.15. A contour plot of the potential energy surface for the Ag₈ model. X and Y axes indicate the increase in axis length from the spherical length in atomic units, a0. Contour lines are 0.2eV apart and hash marks are directed downhill. For this closed shell species there is only one minimum. <full image>

Figure 4.16. A contour plot similar to Figure 4.15 but for Ag₁₅. Contour lines are 0.05eV apart and hash marks are directed downhill. The 7 additional electrons added to the 8 electrons of the closed shell reduce the symmetry of the system. There are three distorted symmetry axes that correspond to permutations of the 3 ellipsoidal axes and are related through volume conservation: R03=RxRyRz. The point (0,0) near the center of the figure is a local maximum. There are at least 3 unique local minima on this surface marked with labels 1-3. <full image>

In contrast to Ag₁₅, Ag₂₁ has an almost spherical stable geometry located in the center of Figure 4.17. There is also a meta-stable configuration corresponding to a more distorted cluster whose total energy is much higher (~0.8eV) than the global minimum.

Figure 4.17. Plot similar to Figures 4.15 and 4.16 for Ag₂₁. Contour lines are 0.1eV apart. The 3 additional electrons added to the 18 electrons of the closed shell introduce one unique local minimum with three unequal axes. Seen here are its six permutations. The global minimum is not spherical but is near (0,0) on this scale. For comparison, the IP's (not the total energies) are 4.21eV at the global minima near the center and 4.76eV at the six symmetrically equivalent local minima. This is a familiar pattern observed in the "shell+1" electron clusters. <full image>

Although there is considerable evidence for structural isomers of non-metal clusters such as Si ⁸⁷ and C, there is less evidence for isomers of metal clusters. Structural isomers of metal clusters have been suggested from reactivity studies⁸⁸ and some recent work on aluminum cluster ion mobilities has detected isomers.⁸⁹ A comparison of those results to this model has not yet been investigated although it may be of interest for future study. It is likely that 'shape' isomers could be distinguished spectroscopically by their different photoabsorption properties. The surface plasmon splits into components based on the symmetry of the three ellipsoidal axes.⁸³ An RPA (random phase approximation) calculation would be difficult in the low symmetry of the ellipsoid but may be worth the effort. Although less likely, the various isomers may also fragment differently in the process of photofragementation.

In reality the PES of a metal cluster includes the local ion-ion forces in addition to the delocalized electronic forces noted here. The local forces map out a high dimensional surface that is certainly covered with local minima. These results suggest the presence of electronic forces that emphasize certain meta-stable 'shapes' besides ion-ion forces that produce meta-stable atomic configurations. The combined forces, in a manner that is not obvious, determine the final atomic geometry. There are also higher multipole deformations and electronic configuration interaction near curve crossings that will perturb these states. Regardless of the precise details, the delocalized electronic forces should be important as the clusters form and cool, thus influencing final geometries.

Figure 4.18. A comparison of the ionization potentials of the local minimum energy geometries with the highest (Ä ) and lowest (—) IP's. Both curves represent stable structures. The IP's of the global minimum energy () configurations belong to neither curve consistently. Note the contrasts at Ag₉, Ag₂₁, and Ag₃₄. These clusters are observed experimentally to have long PIE tails. <full image>

Figure 4.18 depicts the range of IP's that were observed in the model. A key point is that experimental IP's have been compared to the maximum IP shapes (dark curve). There are indications that a number of isomers of silver clusters exist in the molecular beam as evidenced in the long tails of the PIE spectra. Further analysis is required before anything conclusive can be said regarding the existence of isomers. Interestingly enough, the IP's of sodium clusters⁵⁰ follow global minimum energy IP's. For example, Na₃₄ does not have a locally high IP and Na₄₂ has a low IP. This would certainly agree with suggestions that sodium clusters are rather liquid-like in comparison to silver clusters which would be frozen into local depressions on the PES. A better comparison would require additional calculations with parameters for sodium clusters.

Computational Method

The most natural basis on which to evaluate the eigenvalues would presumably be a product of spherical Bessel functions and numerically integrated radial functions. This may be possible; however the evaluation of matrix elements on this basis is not trivial and there would be coupling between many of the basis functions. The effect of quadrupolar distortions is to mix the radial and angular functions, reducing the symmetry of the problem and removing any simple separation of variables.⁹⁰ To avoid this problem a preliminary basis set is used employing sine functions that distort with the cluster. (A similar idea would be to use an FFT to evaluate the kinetic energy in momentum space,⁸³ although preliminary indications are that this is less efficient.) This particular implementation started from a simple model based on an electron in a box. It evolved as more terms were needed for accurate eigenvalues and as the need for speed arose when the potential energy surface was found to hold local minima.

For a Woods-Saxon potential with radius R0 distorted to dimensions Rx,Ry,Rz a set of basis functions distinguished by indices (nx,ny,nz = 1,2,...maxN) was used. Specifically,

(4.9)

These are effectively the solutions for a particle in a three dimensional box. They extend to Lx,Ly,Lz and are optimized to fit the potential. Near optimal values of Li=Ri+0.66rs were found by minimizing the total energy of a 25-electron cluster.

On this basis the kinetic energy portion of the Hamiltonian matrix is diagonal and easily computed. If i and j are used as indices referring to all maxN3 combinations of nx,ny,nz, then the kinetic energy matrix elements are,

(4.10)

The potential energy portion of the matrix,

(4.11)

is not diagonal and involves more computation. Nonetheless, it is independent of the ellipsoidal distortion. This is obvious if, for example, the integration coordinate is scaled to remove the L dependence,

(4.12)

These elements are often zero from the D2h symmetry of the triaxial ellipsoid that has symmetry planes s x, s y, s z. (For example nx and nx' must both be either even or odd for a non-zero element.) As a result the Hamiltonian matrix is separated into its eight irreducible representations. This reduces the computational effort by more than a factor of eight, from maxN2 to (maxN/8)2. When non-zero elements are encountered, a simple trapezoidal integration is used for evaluation.

A number of features in the model allow for rapid computation of the eigenvalues of the Woods-Saxon ellipsoid. First, the kinetic energy term in the Hamiltonian is already diagonal in this basis. Second, using a basis that distorts with the cluster implies that the potential energy matrix elements are independent of the cluster distortion. As a result the off-diagonal matrix elements need only be computed once and stored. Each energy determination involves only a modification of the diagonal kinetic energy elements and a matrix diagonalization.

The technique presented so far is direct and easily programmed. However a thorough investigation of the potential energy surface became necessary when it was discovered that local minima exist on the two dimensional energy surface. This required a more efficient evaluation of the total energy and was accomplished by minimizing the size of the matrix to be diagonalized. The initial sine-basis matrix, although easily computed, must be large to obtain accurate eigenvalues. By projecting onto a new basis that represents the final solutions more closely, a smaller matrix can be diagonalized.

To avoid repeated diagonalization of the large matrix for each set of distortion parameters, H is diagonalized once for the spherical cluster and then the eigenvalues and eigenvectors, E and Y , are stored and used as the new basis. On this basis the matrix elements of the distorted Hamiltonian H' are,

(4.13)

Expressing this as matrix elements,

(4.14)

one discovers that the first term has already been computed (in E) and that the final term is zero since the potential energy is independent of ellipsoidal distortions. The middle term is not difficult to compute since in the original basis the distortion effects only diagonal kinetic energy terms,

(4.15)

Finally the matrix H is assembled and the eigenvalues E' are computed on the new basis,

(4.16)

These are sorted according to energy and filled, two electrons per orbital, to compute the total energy. Since the Hamiltonian on the new basis more closely approximates the final solution, it can be truncated using the approximate energies Ei as a guide. Eigenvectors (or orbitals) with energies greater than +10eV were discarded.

For additional calculations on the same cluster, the process starting with Equation 4.13 was repeated. For these particular calculations, a random initial geometry was chosen and steps were taken in the two independent dimensions to follow the descent of the potential energy surface. To determine the global minimum this process was repeated until 10 additional iterations failed to find a lower total energy.

The computation runs rapidly on a PC. For each cluster the initial integration stage requires a few minutes. Subsequent total energy calculations on the same cluster require approximately 5-20 seconds on a 50MHz 486DX/2. Naturally, the time required increases rapidly as the basis set maxN is increased. The optimal value of maxN depends on the number of electrons, the desired accuracy and the patience of the user. Typically maxN was increased until the IP was converged to within 0.01eV. With this procedure maxN was set to 12 for a 100 electron cluster.

The program is written in C++ (approximately 1000 lines) and contains some inline assembly language for critical routines. The diagonalization routine is based on published code⁹¹ and algorithms.⁹² The code is available from the author on request.

57) K. Laasonen, A. Pasquarello, R. Car, C. Lee, D. Vanderbilt, Phys. Rev. B, 1993, 47, 10142.

58) R. Car, M. Parinello, Phys. Rev. Lett., 1985, 55, 2471.

59) W. D. Knight, K. Clemenger, W. A. de Heer, W. A. Saunders, M. Y.Chou, M. L. Cohen, Phys. Rev. Lett., 1984, 52, 2141. See also ref. 60.

60) L. Martins, R. Car, J. Buttet, Surf. Sci., 1981, 106, 265.

61) K. Clemenger, Phys. Rev. B, 1985, 32, 1359.

62) K. Clemenger, Ph.D. thesis, University of California, Berkeley, 1985.

63) W. Ekardt, Phys. Rev. B, 1984, 29, 1558.

64) W. A. de Heer, Rev. Mod. Phys., 1993, 65, 611.

65) M.D. Morse, Chem.Rev., 1986, 86, 1049.

66) I. N. Levine, Quantum Chemistry, Allyn and Bacon, Massachusetts, 1983.

67) George F. Bertsch, The Practitioner's Shell Model, North-Holland, 1972.

68) Eat potatoes if the pork is bad.

69) W.-D. Schöne, Diploma thesis, Bielefeld University, 1991.

70) M. Brack, Rev. Mod. Phys., 1993, 65, 677.

71) H. A. Jahn, E. Teller, Pro. R. Soc. London Ser. A, 1937, 161, 220.

72) D.M.P. Mingos, T. Slee, L. Zhenyang, Chem. Rev., 1990, 90, 383.

73) S. G. Nilsson, (Danish) Det Kongelige Dananske Videnskabernes Selskab Matematisk-fysiske Meddelelser 1955, 29, 16.

74) K. Selby, M. Vollmer, J. Masui, V. Kresin, W. A. de Heer, W. D. Knight, Phys. Rev. B, 1989, 40, 5417.

75) W. A. Saunders, Ph.D. thesis, University of California, Berkeley, 1986.

76) A. Bohr, B. R. Mottelson, Nuclear Structure, 1985, W. A. Benjamin, Reading, Massachusetts; section 5-1.

77) M. A. Preston, R. K. Bhaduri, Structure of the Nucleus, Adison-Wesley, Reading, Massachusetts, 1975; section 9-9.

78) Z. Penzar, W. Ekardt, Z.Phys.D., 1991, 19, 109.

79) W. Ekardt, Z. Penzar, Phys. Rev. B, 1988, 38, 4273.

80) E. Schumacher, M. Kappes, K. Marti, P. Radi, M. Schär, B. Schmidhalter, Ber. Bunsenges. Physik Chem., 1984, 88, 220.

81) G. Makov, A. Nitzan, L.E. Brus, J. Chem. Phys., 1988, 88, 5076.

82) A. Bulgac, C. Lewenkopf, Phys. Rev. Lett., 1993, 71, 4130.

83) G. Lauritsch, P.-G. Reinhard, J. Meyer, M. Brack, Physics Letters A, 1991, 160, 179.

84) R. Poteau, F. Spiegelmann, J. Chem. Phys. 1993, 98, 6540.

85) M.L. Cohen, M.Y. Chou, W.D. Knight, W. A. de Heer, J. Phys. Chem., 1987, 91, 3141.

86) A. Bulgac, C. Lewenkopf, Phys. Rev. Lett., 1993, 71, 4130.

87) M.F. Jarrold, and J.E. Bower, J. Chem. Phys., 1992, 96, 9180.

88) M.B. Knickelbein, S.Yang, J. Chem. Phys., 1990, 93, 5760.don't have

89) M.F. Jarrold, J.E. Bower, J. Chem. Phys., 1993, 98, 2399.

90) Analyses using perturbation theory avoid this but require other approximations. See Ref. 72.

91) W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in C, Cambridge University Press, 1992, routines tred2() and tqli(); See also LAPACK and EISPACK Guide.

92) G. H. Golub, C. F. Van Loan, Matrix Computations, John Hopkins University Press, 1989.