Reassignment of magic numbers for icosahedral Au clusters: 310, 564, 928 and 1426

Icosahedral Au clusters with three and four shells of atoms are found to deviate significantly from the commonly assumed Mackay structures. By introducing additional atoms in the surface shell and creating a vacancy in the center of the cluster, the calculated energy per atom can be lowered significantly, according to several different descriptions of the interatomic interaction. Analogous icosahedral structures with five and six shells of atoms are generated using the same structural motifs and are similarly found to be more stable than Mackay icosahedra. The lowest energy per atom is obtained with clusters containing 310, 564, 928 and 1426 atoms, as compared with the commonly assumed magic numbers of 309, 561, 923 and 1415. Some of the vertices in the optimized clusters have a hexagonal ring of atoms, rather than a pentagon, with the vertex atom missing. An inner shell atom in some cases moves outwards by more than an \AA{}ngstr\"om into the surface shell at the vertex site. This feature, as well as the wide distribution of nearest-neighbor distances in the surface layer, can strongly influence the catalytic properties of icosahedral clusters. The structural optimization is initially carried out using the GOUST method with atomic forces estimated with the EMT empirical potential function, but the atomic coordinates are then refined by minimization using electron density functional theory (DFT) or Gaussian approximation potential (GAP). A single energy barrier is found to separate the Mackay icosahedron from a lower energy structure where a string of atoms moves outwards in a concerted manner from the center so as to create a central vacancy while placing an additional atom in the surface shell.


Introduction
Gold nanoclusters have many intriguing properties and can be used in various applications such as catalysis [1][2][3][4] , plasmonics 5,6 and biosensors. 7 In order to gain an understanding of their functionality, the atomic structure of the clusters needs to be identified. A striking example of the variation of functionality with size is the remarkably facile oxidation of CO on Au clusters with diameter between 2 to 3 nm adsorbed on TiO 2 . 4 Clusters outside this size range, both larger and smaller, have been found to be less active per surface site. The ordering of the atoms in the clusters, in particular on the cluster surface, strongly influences the catalytic properties. Low coordinated Au atoms on the surface are believed to be responsible for the catalytic activity (see, for example, ref. 8 ), but the reduced efficiency for clusters with diameter smaller than 2 nm remains a puzzle.
Small Au clusters are planar but clusters with 13 or more atoms have three-dimensional structure. 9 Various pyramids and cages have been deduced from experimental and computational studies of clusters with up to a few tens of atoms. [10][11][12][13][14][15] The rich size dependence of Au clusters in this size range has been ascribed to relativistic effects. 13 When the interaction between atoms does not depend strongly on the directionality of the chemical bonds, atomic clusters tend to have icosahedral shape and five-fold symmetry up to a certain cluster size. This applies, for example, to rare gas atoms and simple metals, where the interaction energy mostly depends on the distance between atoms. Calculations of clusters where the atomic interaction is described with a simple pairwise Lennard-Jones potential function have shown that the icosahedral structure motif dominates for clusters with up to 5000 atoms. 16 Eventually, for large enough clusters, the atoms order in a way that is characteristic of the crystalline phase as it becomes lower in energy than icosahedral ordering. Mackay presented the atomic structure of icosahedral clusters of increasing size as shells of atoms are added to the basic 13 atom cluster. 17 The surface layer is complete and compact when the number of atoms in the cluster is 13, 55, 147, 309, 561, 923, 1415 . . . These are often referred to as 'magic numbers' since the energy per atom can then be particularly low and they can, therefore, be expected to be more abundant when clusters of various size are generated experimentally.
Clusters with Ino-decahedral 18 and cuboctahedral shapes turn out to have the same magic numbers as the Mackay icosahedra. Marks decahedra 40,41 involve different truncations of the twinned tetrahedra and give rise to different magic numbers. 36,43 For noble gases and simple metals, a correlation is indeed found between the abundance of clusters as a function of size and the magic numbers. Some theoretical and experimental studies have assumed that the same magic numbers apply to Au clusters.
The electronic structure of metal clusters also favors clusters of a certain size corresponding to complete valence electron shells leading to extra stability for clusters with 8, 18, 34, 58, and 92 atoms. [19][20][21] But, for larger clusters this effect becomes relatively less important than the packing of the atoms.
Many studies have been devoted to Au clusters with 55 atoms, one of the traditional magic numbers. The optimal structure was unexpectedly predicted to be mostly amorphous but nevertheless chiral. 10,22,23 Local structure analysis has identified part of it as a distorted icosahedron. Photoelectron spectroscopy measurements of ionic Au − 55 clusters combined with density functional theory calculations 14 as well as aberration corrected scanning transmission electron microscope (AC-STEM) measurements of Au 55 clusters deposited on a amorphous carbon surface 15 have lent support for this structural prediction. The origin of the of the low symmetry of the Au 55 cluster has been ascribed to a strong contractions of the surface of the cluster, thereby creating vacancy islands. As the cluster grows up to Au 64 , additional atoms get introduced into the surface layer and the energy per atom is reduced as a result. 14 This is analogous to gold crystal surfaces, such as Au(100), where the additional atoms get inserted into the surface layer in the so-called hex reconstruction. 24,25 As the cluster size increases, the icosahedral (Ih) structure becomes less stable with respect to the decahedral (Dh) and cuboctahedral structures, the latter having the same local atomic ordering as the face centered cubic crystal, exposing (111) and (100) facets at the surface. The (100) facets are higher in energy so a lower energy structure is generally obtained by truncating the cuboctahedron in such a way as to reduce the area of the (100) facets to form a truncated cuboctahedron (TCo). It has often been assumed that the crossover from one structural type to another is a sharp function of cluster size and that it follows the progression Ih → Dh → TCo. 26,31 The determination of these crossovers has been the topic of many computational studies and various predictions have been made. The crossover from Ih to Dh has generally been predicted to occur for relatively small clusters, with less than 150 atoms. A larger variety of predictions have been made for the Dh to TCo crossover, with estimates ranging from 200 to 700 atoms [31][32][33][34] or even larger than 5000 atoms. 35,38 Other calculations have, however, indicated that the crossover between structural motifs is not so simple. 36,37,39 While the energy per atom of Dh clusters varies smoothly with size, as the structure can be that of Ino or Marks decahedra or some variation on those motifs, there are larger fluctuations for the TCo motif where the regular crystalline facets provide the highest stability. This results in alternations of the most stable structural motif over a wide size range, up to around 4000 atoms. 37 As a result, different magic numbers are obtained for Dh and TCo clusters than one would expect from the shell closing of the ideal structures mentioned above.
Recent advances in experimental imaging techniques, in particular the development of the aberration corrected scanning transmission electron microscope (AC-STEM), have led to direct observations of three-dimensional, atomic-scale structure of Au clusters. [26][27][28][29][30] The AC-STEM method has been used to study larger Au clusters deposited on a support, even as large as 1000 atoms. While the preparation method can be tuned to promote the formation of certain structures, 29 the inherent metastability of clusters and energy barriers for transitions from one structure to another make it difficult to attribute experimentally observed isomeric populations to equilibrium populations. 30 A further complication is the possible temperature dependence of the optimal structure due to entropic effects. 41,42 The atomic-scale structure of Au clusters with up to around 1000 atoms has been investigated extensively in recent years by Palmer and coworkers using AC-STEM. 15,26,[28][29][30] All three structural motifs, Ih, Dh and TCo, have been observed over this size range. The icosahedral structure is deduced to be the least stable of the three motifs for clusters with 923 atoms, since it has lower abundance 26,29 and is found to transform to one of the other motifs under beam irradiation. 28,30 Determining the relative stability of the Dh and TCo motifs has turned out to be more difficult. For Au 923 , the Dh and TCo clusters generally retain their structural motif under irradiation suggesting that both are stable structures 28 but for Au 561 the Dh isomer was found to transform to TCo. 30 Interestingly, the populations of Dh and TCo for 561, 742 and 923 atom Au clusters were found to be nearly the same, within experimental uncertainty. Based on classical dynamics simulations of cluster growth where transformations between structures were not observed during the growth process, it was concluded that the structure of a cluster becomes 'locked in' at a small size after which templated growth takes place without transformation to a more stable structure. 30 An alternative explanation can be formulated based on the close but alternating stability of the Dh and TCo structures with cluster size, mentioned above. 37 Despite the large amount of research to date, the motif preference for medium size (300 to 1000 atoms) Au clusters remains an open question. Furthermore, the prediction of the most stable atomic structure for a given structural motif has turned out to be a challenge. We study here the atomic structure of Au nanoclusters with the icosahedral motif in this size range and search for those that have particularly low energy per atom. The global minimization of the energy with respect to atom coordinates for clusters of this size is challenging and it is hard to be sure that the optimal solution has been found, but it is clear from the results obtained that the usual magic numbers do not apply to Au clusters.

Methods
The structure optimization of the clusters with three and four atom shells is carried out with a method referred to as global optimization using saddle traversals (GOUST). 55 There, the system is brought from one local minimum on the potential energy surface to another by identifying low lying first order saddle points. From each local minimum, several saddle point searches are conducted and the local energy minimum on the opposite side of the energy ridge is found by minimization after a slight displacement along the eigenvector corresponding to the negative eigenvalue of the Hessian matrix. The system is then advanced to the lowest local minimum found in the various saddle point searches carried out for this structure. This process is repeated for several different random initial configurations of the atoms.
The saddle points are found using the minimum mode following (MMF) method. [50][51][52] The EON software 53 and the EMT interaction potential function 54 (with Rasmussen parameter values) are used for the initial phase of the structural optimization. Subsequently, a local minimization of the energy is carried out using density functional theory (DFT) with either the PBE 56 or the PBEsol functional approximations. 57 The PBEsol functional is known to give better estimates of the surface energy of metals, so it is expected to be more accurate for small clusters where the surface energy is of great importance. The difference in energy obtained with the two functionals can be used as an indication of the accuracy of the DFT calculations. A plane-wave basis with a 230 eV kinetic energy cutoff was used to represent the valence electrons and the projector augmented wave formalism 58 for the inner electrons, as implemented in the VASP software. 59,60 Further analysis was carried out using a Gaussian approximation potential (GAP) 61,62 , fitted to an extensive data set of PBE calculations, including various types of crystal structures, surface slabs and clusters with random disorder (which mimics thermal excitation). GAP is a datadriven approach to interpolating the interatomic energies and forces based on Gaussian process regression. 63 The GAP potential used here for Au atoms 64 incorporates two-body and atom-centered many-body descriptors 65,66 to encode the local atomic environments. This leads to an atomwise decomposition of the total energy which can be used to analyze further the driving force for the structural changes in the gold clusters during optimization. Starting with the cluster structures obtained by GOUST and EMT, the GAP energy is minimized with respect to the atomic coordinates and the energy per atom evaluated. This also allows for a comparison of the DFT and GAP total energy changes to assess the accuracy of the machine learned potential function.
Guided by the results obtained for clusters with three and four shells, atomic structures were constructed for icosahedral clusters with five and six shells of atoms. This involved inserting rows of additional atoms in various ways and then minimizing the energy. The minimization was carried out using the GAP potential as well as with DFT.

Results
After carrying out a few tens of saddle point searches starting with the 309 atom Mackay icosahedron and small initial random displacements of the atoms, a concerted displacement of 6 atoms was found where the central atom moves into the first (12 atom) shell, taking the place of an atom that moves from the first shell to the second (42 atom) shell and the atom initially in the second shell moves into the third shell, thus adding a new atom to the surface of the cluster. This transition corresponds to a single saddle point on the energy surface. A similar process involving a concerted displacement of 7 atoms was found for the larger 561 atom Mackay icosahedron. The energy along the minimum energy path and an illustration of the displacement of the atoms of the 561 atom cluster is shown in Fig.  1. The final state of this transition is 0.8 eV lower in energy as both the compressive strain on the central atom and the tensile strain in the surface shell is relieved. As can be seen from the energy curve in Fig. 1, there is just a single maximum along the minimum energy path for this multi-atom transition mechanism. Subsequent steps in the global optimization result in the formation of vacancy patches in the surface shell as surface atoms group together, making room for additional atoms to be incorporated into the surface shell. Sequential addition of atoms in various places in or near the surface shell followed by structural optimization shows that the lowest energy per atom is obtained by adding one atom to the 309 atom cluster  . 1 Energy along the minimum energy path for a transition in a 561 atom Au cluster where the initial structure is the energy minimized Mackay icosahedron and the final state (shown in inset) has a vacancy in the center and an extra atom in the surface shell. The transition involves a concerted displacement of seven atoms (shown in red) and corresponds to a single saddle point on the energy surface. The saddle point was found using the minimum mode following method and EMT starting from small random displacement of the atoms from the Mackay structure, and the steepest descent path then found in both directions down from the saddle point. Blue spheres represent atoms in the first shell and red spheres represent the seven displaced atoms. and adding three atoms to the 561 atom cluster. The optimal structures were further refined by local minimization using DFT or GAP, followed by analysis of the structure and energy per atom. Fig. 2 shows the change in the energy of each atom as the central atom is removed from the 561 atom cluster. The lowering of the total energy corresponds to 1.6 meV per atom. The removal of the central atom results in lowering of the local strain. The effect is largest in the first shell of 12 atoms and, to a lesser extent in the surface shell. In the two subsurface atomic layers a handful of atoms increase slightly in energy. Table 1 summarizes the main characteristics of the 3 and 4 shell icosahedral clusters obtained in this way. The larger cluster is shown in Fig. 3.
The two surface facets containing extra rows of atoms are highlighted.
The DFT calculated energy per atom in these clusters compared with that of the FCC crystal, i.e., the excess energy per atom, is given in Table 2 and compared with Mackay icosahedral structures where the central atom has been removed but no extra atoms added to the surface shell. The addition of atoms to the surface shell is found to lower the energy of the cluster with four shells (containing 560 vs. 564 atoms) by ca. 10 meV per atom. Similar energy lowering due to the surface reconstruction and addition of surface atoms is found with the two density functionals, PBEsol and PBE, and with the GAP potential.
The vertex sites adjacent to facets containing an extra row of atoms are quite different from the vertex sites of a Mackay icosahedron. In the 564 atom cluster, such sites have the vertex atom sitting on a distorted hexagon instead of a pentagon. Vertex sites adjacent to two facets with -12 +12 0 +6 -6 Local energy difference (meV)   Table 1. Facets where additional atoms have been inserted are shown in blue. Vertex atoms at anomalous sites and belonging to the surface shell are shown in white, while subsurface atoms are shown in gray (seen when vertex atoms are missing). Atoms that have moved out from a subsurface shell to a location close to the center of a hexagonal ring are shown in red. The three types of anomalous vertex sites are shown in Fig. 4. Table 1 Characterization of the optimal structure of 3, 4, 5 and 6 layer Au icosahedra. In all cases, an atom has been removed from the central site to create a vacancy and atoms have been added to the surface shell so as to find the lowest energy per atom in DFT/PBEsol calculations. The anomalous vertex sites are defined and shown in Fig. 4    an extra row are vacant as the vertex atom has become a surface facet atom but a subsurface atom has been displaced outwards into the surface shell to a location near the center of the hexagon, see Figs. 4 and 5. The distance between this displaced atom and the underlying atom is large, 4.4Å, compared with the near neighbor distance of 2.93Å in the Au crystal. Anomalous vertex sites in the 310 atom cluster resemble those of the 564 atom cluster, except that the subsurface atoms are not displaced outwards at the vacant vertex sites. Fig. 3. shows the atomic structure of the three types of anomalous vertex sites. The intermediate size clusters, 564 and 928, are special in that they have vertex Au atoms that are only sixfold coordinated. Global optimization simulations of the larger, 5 and 6 shell, icosahedral clusters starting from the Mackay icosahedra is more challenging because of the large number of atoms. But, we have generated structures for these clusters by assuming similar motifs as for the smaller clusters, i.e. the removal of the central atom and introduction of extra surface atoms to form additional rows in some of the facets. Subsequent GOUST optimizations  Fig. 6 Distribution of nearest neighbor distances between surface atoms in the optimized icosahedral clusters (blue) and in the Mackay icosahedra with a central vacancy (red). The energy has been minimized using DFT/PBEsol. The dashed black line marks the calculated near-neighbor distance in the Au crystal at the same level of theory. The optimized icosahedral structures have a broader distribution and a smaller average distance between neighbors than the Mackay icosahedra.
were carried out using EMT, followed by local minimization using GAP or DFT. The optimal cluster structures found in this way are compared with Mackay icosahedra in Table 1 and illustrated in Fig. 3. The lowering of the energy per atom as compared with the corresponding Mackay icosahedron with a central vacancy is 8 meV per atom using DFT/PBEsol and 10 meV per atom using DFT/PBE for both 5 and 6 shell clusters (see Table 2). The GAP potential gives somewhat smaller change in the energy but the same trend. The largest cluster, containing 1426 atoms, has regained certain symmetry, as illustrated in Fig. 3. All 12 vertex sites are vacant there.
The distance between nearest neighbors in the surface shells of the optimized structures varies greatly, much more than in Mackay icosahedra, as illustrated in Fig. 5. Nearest neighbor distances ranging from less than 2.8Å to larger than 3Å can be found in the optimized structures. On average, the nearest-neighbor distances are smaller in the optimized structures as would be expected from the increased number of atoms in the surface shell.

Discussion
The optimal structures presented here for icosahedral Au clusters with 3, 4, 5 and 6 shells of atoms are significantly lower in energy than the Mackay icosahedra. The 'magic numbers' turn out to be 310, 564 , 928 and 1426 instead of the Mackay magic numbers of 309, 561, 923 and 1415. Since knowledge of the atomic structure is the basis for understanding the various properties of the clusters, such as catalysis, this is an important finding.
The fact that the surface shell includes additional atoms as compared with the Mackay structure is consistent with the experimental findings for smaller gold clusters with 55 to 64 atoms 14 and the known surface reconstruction of Au crystal surfaces. 24 The tendency of Au clusters and the gold crystal surfaces to accept additional atoms in the outermost layers can be related to relativistic effects. 14 The unusual atomic ordering at some of the vertex sites and the large variation in nearest-neighbor distances is likely to have significant effect on the catalytic activity of these clusters. The anomalous vertex sites formed where one or two of the adjacent facets include an extra row of atoms are particularly interesting. Previously, sites with a missing vertex atom and an open hexagonal ring had been found from DFT calculations of 55 atom Pt clusters. 67 There, these sites were referred to as 'rosettes'. Here, we also observe the possibility that an underlying atom moves outward to a location that is near the center of the hexagonal ring. This is the case for clusters of intermediate size, 564 and 928. The atom near the center of a hexagon is only sixfold coordinated, since the distance to the nearest inner shell neighbor is 4.4Å. Such low-coordinated atoms are likely to provide higher binding energy for adatoms and admolecules, and also have significant structural flexibility since a displacement toward the center does not involve large increase in energy. Such structural flexibility could be of importance in chemical reactions.
The saddle point searches used in the optimization method applied here readily revealed the formation of a vacancy at the center of the clusters. This lowers the energy per atom because the center is under compressive strain in large icosahedral clusters. Such vacancy formation had previously been predicted for icosahedral Pt clusters based on tight binding calculations 48 and is likely a common feature in icosahedral metal clusters.

Acknowledgements
This work was supported by the Icelandic Research Fund and the Academy of Finland under grants number 329483 (J.K. and M.A.C.) and 330488 (M.A.C.). Computational resources from CSC -IT Center for Science (NO-BLE Grand Challenge project) and Aalto University's Science IT Project are gratefully acknowledged as well as the Icelandic Research High Performance Computing facility. We thank Anna L. Garden for helpful discussions.