Computational Characterisation of Structure and Metallicity in Small Neutral and Singly-Charged Cadmium Clusters

Computational Characterisation of Structure and Metallicity in Small Neutral and Singly-Charged Cadmium Clusters† Pablo ´ Alvarez-Zapatero and Andr´ es Aguado∗ Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX First published on the web Xth XXXXXXXXXX 200X DOI: 10.1039/b000000x Putative global minimum structures for neutral CdNand singly charged Cd+ Nand Cd− Nclusters in the small size regime up to N=21 atoms are reported. A global optimization approach based on the basin hopping method and a Gupta potential fitted to cluster properties is employed to generate a diverse databank of trial structures, which are then re-optimized at the densityfunctional level of theory. Novel, previously unreported, structures are found for many sizes. Our results successfully reproduce and interpret the size-dependent stabilities known from mass spectrometry, and strongly suggest that experiments aimed at determining the relative stabilities of neutral cadmium clusters are really measuring cation stabilities. We provide an in-depth analysis of electronic structure and use it to explain the gradual emergence of metallic-like behaviour as the cluster size increases. 1 Introduction Achieving an accurate structural characterization of atomic and/or molecular clusters is a crucial problem that needs to be solved in order to get a sound understanding of the physicochemical properties of these interesting nanoscale systems. This is true at an applied level, as the specific magnetic and optical responses, or the catalytic activities, of clusters (just to mention a few properties of interest in real applications using functional clusters) depend on the geometry of the atomic skeleton, and only with a detailed understanding of structure/property relationships will science be able to propose the ideal cluster system for a particular application. It is true at a more fundamental level as well, as clusters represent a sort of bridge connecting the atomistic and macroscopic worlds, and are thus ideal systems to analyse the emergence of typically bulk properties and statistical complexity as the number of atoms increases. In this connection, the cadmium clusters which are the object of study in this paper can be important in industrial terms as components of corrosion-protecting coatings1, and these protecting properties could be fine-tuned by nanoalloying it with zinc and locating an optimal composition for the nanoalloy. As Cd is a group XII metal with a closed-shell electronic configuration, cadmium clusters are also considered to be ideal model systems to study the nature of the insulator- to-metal transition that must occur along the path bridging the †Electronic Supplementary Information (ESI) available: Atomic coordinates (in xyz format and ˚ A units) and point group symmetries for the Global Minimum structures reported in this paper. See DOI: 10.1039/b000000x/ Departamento de F´ ısica Te´ orica, At´ omica y ´ Optica, University of Valladolid, Valladolid 47071, Spain; E-mail: [email protected] dimer with the bulk limit2. However, the problem of characterizing the structure of cadmium clusters has been addressed only in a very limited number of computational works3–10, which moreover do not generally agree in the proposed global minimum (GM) structures. Experimental results about relative stabilities of charged cadmium clusters are also available11–13, which have not yet been convincingly reproduced and interpreted. In this work, we report putative GM structures for Cd clusters in the small size regime up to 21 atoms, for both neutral (CdN) and singly-charged (Cd+ Nand Cd− N) clusters. The structures are obtained by a two-step method based on empirical potentials and density-functional-theory (EP-DFT). We locate for many sizes putative GM structures which are more stable than previously reported. The obtained cluster stabilities are found to reproduce the experimental mass spectra of charged clusters, wherefrom a physical interpretation of those spectra emerges. Finally, a detailed analysis of the electronic structure is offered to analyse the evolution of metallicity as a function of increasing size in the small size regime. 2 Computational Methods Structures of cadmium clusters are calculated in a two step procedure; the first one involves a global optimization using an empirical potential (EP) which provides a wide scan for local minima in the potential energy surface. Subsequently, the most promising structures obtained for each size are reoptimized at Density Functional Theory (DFT) level. The reliability of this combined EP/DFT approach relies on the quality of the structures provided by the global optimiza- 1–25 |1 tion during the first step. Since the DFT energies are obtained from a local optimization of trial EP structures, these should not be very dissimilar to the real ones. Thus, the EP should capture the physical behaviour of the real material at a reasonable level. We employ the Gupta potential14 defined as: E(CdN) = N ∑ i (Ea i+Er i) Ea i=−(N ∑ i6=j ξ2exph−2q(ri j r0 −1)i)1/2 Er i= N ∑ i6=j Aexph−p(ri j r0 −1)i (1) where ri j is the distance between atoms iand jand A,ξ,p,q,r0are constants traditionally fitted to experimental bulk quantities. In a recent work15, a simpler expression for this potential has been derived, by eliminating first the redundant parameter r016 and then adopting reduced energy and distance units. In brief, the reduced Gupta potential contains just two independent parameters, chosen as λ=q/pand ξ=Aep χeq. In terms of these new constants, the potential takes the form15: e E=1 2 N ∑ i(N ∑ i6=j χexp[eri j]−hN ∑ i6=j exp[−2λeri j]i1/2)(2) with the reduced distances and energy defined as eri j =pri j r0 and e Ei=Ei 2ξeq, respectively. We have presented in our previous works on zinc clusters15,17 a detailed protocol for finding the optimal λand χparameters for a given metal, which are expected to offer realistic trial structures to feed the DFT reoptimizations. In this work, however, we have decided to directly use the same reduced Gupta potential as for zinc clusters, rather than repeating the whole protocol which is quite complex. The rationale is that, Zinc and Cadmium being isovalent, we expect their bonding pattern to be similar, the main differences being related to different bond lengths and energy scales, to which the reduced Gupta potential is not sensitive anyways. Moreover, experimental observables such as abundance mass spectra or ionization potentials (discussed in the next sections) reveal a very similar size evolution for Zn and Cd clusters, providing further credentials for our approximation. Nevertheless, we admit that the accuracy of our strategy can only be confirmed a posteriori, and in this respect we will show that (a) our EP-DFT approach locates more stable minima than previous EP-DFT calculations for many sizes, which already confirms that the seed structures proposed by our EP model are better approximations to the DFT global minimum; (b) the DFT reoptimizations produce results in agreement with experimental observables. Taken together, these two results assess the accuracy of the whole protocol. Around 100 isomers were then selected from the Gupta database for each cluster size and reoptimized at the DFT level. In a final step, performed directly at the DFT level, we have manually generated structures for CdNby removing atoms with low coordination from the most stable CdN+1 structures or by adding one adatom to the CdN−1structures. Both this final step and the consideration of a large number of isomers for DFT re-optimization further reduce the possible errors due to approximate EP parameters. The DFT calculations had been performed with the SIESTA code18, under the PBE approximation19 to exchangecorrelation effects. An auxiliary Van der Waals forcefield was added to the DFT energy, of the form EVdW = −∑i<jf6(Ri j)C6R−6 i j , where f6(r) = 1−e−br ∑6 k=0 (br)k k!is a Tang-Toennies damping function20. The constant C6=279 eV ˚ A6is taken from a high-level ab initio study21 and the damping parameter is set to b=1.2˚ A−1. We have benchmarked this method by comparing Cadmium dimer and crystal phase features with those calculated with an explicitly non-local vdW-DFT exchange-correlation functional. In this case, the KBM functional22 was used (see Table I). Since both methods result in the same good accuracy, we decided to employ the “PBE plus Van der Waals auxiliary force-field” method, named PBE-D in Table I, due to its lower computational cost. We included 3dand 4selectrons into the active valence set, while the effect of core electrons were described with normconserving Troullier-Martins pseudopotentials27. The resolution of the real-space grid employed in Fast Fourier Transforms was defined by a cutoff energy of 300 Ry. The electronic wave function was expanded on a basis of localized atomic orbitals, including s,pand dangular symmetry with double-zeta plus polarization (DZP) quality. We explicitly checked that moving to DZP2 (i.e. doubling the number of polarization orbitals) or TZP (three different radial functions per angular momentum channel), modifies the near-neighbor distance by less than 0.01 ˚ A, and the cohesive energy by at most 0.02 eV, as compared to the DZP results shown in Table I. So the agreement with experiment would indeed be slightly better with a TZP2 basis set, but the analysis shows that the DZP basis is reasonably complete. We decided to show in Table I the DZP results as this is the basis employed in the cluster calculations. All equilibrium cluster geometries reported were obtained from unconstrained conjugate-gradients structural relaxation using DFT forces. The structures were relaxed until the force on each atom was smaller than 0.01eV /˚ A. 2|1–25 Table 1 Equilibrium bond distance dand dissociation energy Edis of the Cd2dimer, and lattice constants (a,c) and cohesive energy Ecoh of bulk Cd, as predicted by several exchange-correlation functionals, are compared to either experimental results or other theoretical calculations. Because experimental values about the spectroscopic constants of the dimer do not agree with each other, we choose to show high-level coupled-cluster benchmark results for comparison purposes, taken from Reference23, and the set of experimental results that shows the best agreement with the theoretical calculations, taken from24. Experimental results for bulk Cd are taken from Kittel25, and ab initio correlated calculation results for bulk Cd from26. Cd2d(˚ A) Edis (eV) PBE 3.840 0.051 PBE-D 3.765 0.043 KBM 3.761 0.037 CCSD(T) 3.873 0.040 Exp. 3.78 0.041 Cd(bulk) a,c(˚ A) Ecoh (eV) PBE 3.023, 5.865 0.80 PBE-D 2.984, 5.696 1.08 KBM 3.019, 5.875 1.09 Ref.26 3.04, 5.65 1.26 Exp. 2.98, 5,62 1.16 3 Results and Discussion 3.1 Putative Global Minimum Structures The putative GM structures that emerge from our study are shown in Figures 1 to 6. We shall start describing the neutral clusters, followed by cations and anions, focusing on the main differences that arise from the charge difference. Figures 1-2 show the neutral GM structures. The trimer is an equilateral triangle, and the tetramer is a perfect tetrahedron. The GM structures for N=5−8 are all based on different combinations of tetrahedral units: two tetrahedral units share one face for Cd5(forming a trigonal bi-pyramid), an edge for Cd6and a single atom for Cd7. Torsional excitations of the two tetrahedra in Cd7are very close in energy to the GM, so this cluster will be highly fluxional even at low temperatures. Cd8is more compact and is obtained by packing four tetrahedral units. A structural change occurs at N=9, which defines a new growth pattern for clusters with 10 to 16 atoms: the GM of Cd9is a tri-capped trigonal prism (TTP), and for N=10 − 14 the GM structures are obtained by adding atoms to this TTP unit. It is worth noticing, however, that for N=10 a perfect tetrahedron is competitive with the GM structure, and for Cd9a capped square antiprism is also a competitive isomer. Finally, Cd15 and Cd16 are formed by glueing together two TTP units. Cd17 is the first structure exhibiting a pseudo-internal atom, almost completely surrounded by a shell of more external atoms. Thus, Cd17 sets a borderline between structures with and without a core atom. The GM structure of Cd18 is obtained from a 13-atom decahedron with five equatorial adatoms, although severely distorted. Both Cd19 and Cd20 exhibit an inner atom surrounded by a chiral surface with C3symmetry, and removing one atom from Cd19 an isomer is generated for Cd18 which is essentially degenerate with the GM. The GM structure of Cd20 is identical to that previously obtained for Zn20 15,17, Mg20 28 and Na20 29 clusters, and it was the reference motif to explain the GM structures of zinc clusters with up to 24 atoms17. We see that Cd21 indeed displays this Cd20 structure with a low coordinated adatom. Different locations of the adatom result in nearly degenerate structures. The peculiar bonding pattern in Zn21, featuring a coexistence of insulating and metallic phases, was analyzed in a recent report17, and we will show later that the electronic structure in Cd21 conforms to the same description. We shall proceed now with cations, displayed in figures 3-4. The smaller clusters exhibit a lower compactness degree and average coordination number when compared with the corresponding neutral structures. Cd+ 3is still a linear chain, and Cd+ 4has a planar structure in the form of a centered triangle. Similarly, Cd+ 5is an expanded tetrahedron with one atom in its center. Cd+ 6also shows a non compact structure with a singlycoordinated adatom attached to a trigonal bi-pyramid, while Cd+ 7has instead a dimer attached to the same reference bipyramid. Cd+ 8essentially features a single bond between two tetrahedra. For most of these sizes, the presence of low coordinated atoms generates several isomers which are nearly degenerate with the GM structures, suggesting these small cations will all present significant fluxionality. Cd+ 9does not show the TTP unit as GM structure but rather a capped square metaprism (the perfect antiprism with C4 symmetry is slightly less stable). Nevertheless, clusters with 10 to 16 atoms are based on the TTP unit as their neutral counterparts. The main difference between cations and neutrals in this size range is the marked preference of cations for dangling atoms with the lowest coordination, as exemplified by the GM structures of Cd+ 11 and Cd+ 12. Cd+ 17 represents again the borderline between structures with and without a core atom. Its most stable structure displays a quite spherical and highly symmetric shell around the central atom. Cd+ 18 and Cd+ 19 have also high symmetry and are based on decahedral growth. For N=19 the GM structure of the cation is thus different from that of the neutral. Finally, Cd+ 20 and Cd+ 21 are similar to the corresponding neutral clusters. The anion structures are displayed in Figures 5-6. We will focus on those sizes for which the anion GM differs from the 1–25 |3 neutral GM. Cd− 5shows a dangling atom attached to a tetrahedron, while Cd− 6is a pentagonal pyramid, both being less compact than the corresponding neutral structures in terms of average coordination number. For these very small anions, the effect of the additional Coulomb repulsion in the electron density cloud seems to be a dominant factor on the atomic packing. Cd− 8is a perfect square antiprism. For Cd− 9we recover the TTP unit like in the neutral situation, which becomes a substrate for the growth of the majority of anion structures with 10–16 atoms. An exception is the GM structure of Cd− 10 which is rather based on the square antiprism unit with two atoms capping its square facets. Although the 17-atom clusters display different GM structure for each charge state, they share the characteristic of being the smallest structure with an incipient core atom. Cd− 18 is quite amorphous as opposed to both the neutral and cation. Cd− 20, despite sharing the C3symmetry and looking quite similar visually to Cd20 and Cd+ 20, shows indeed some substantial differences: the atom added to the 19-atoms structure initiates now a clear more external atomic layer, while it was integrated in the atomic shell of neutral and cation structures. Consequently, the Cd20 structure is no longer very stable when negatively charged. For the same reason, Cd− 21 is not based on a privileged structure with a low coordinated adatom, and instead shows a competition between different structures, all of them quite amorphous. In the whole size range (with the only exception of Cd− 5), anions do not show a preference for singly-coordinated adatoms and are clearly more compact than cations. For several sizes, they are even more compact than neutral clusters of the same size. Some general conclusions can be outlined from this description: there is no single size for which the GM structures are the same for all charge states, hence the number of electrons has a great influence on the structures, suggesting that electronic effects might also be a dominant factor affecting cluster stabilities. A more specific observation is that singly-coordinated adatoms are most frequently found for cations, so the global cluster charge plays a significant role in stabilizing this strange structural feature. We believe that the electron deficient character of these clusters is the key property needed to rationalize this trend. The Cd clusters studied in this paper are considered electron deficient systems in the sense they have fewer valence electrons than needed to form a 2-center 2-electron (2c-2e) bond between each pair of atoms, so the few valence electrons available have to screen nuclear repulsion through multi-center bonding. The electron deficiency becomes more important the smaller the size because s−phybridization is then less well developed, and may be extreme precisely for cations which have one electron less, favouring less compactly packed structures with exotic adatoms. The larger the size of the cluster the more complete are the s-p hybridization and electron delocalization. Consequently, for larger clusters we expect that structural differences become less and less dependent on the global cluster charge, and that singly-coordinated adatoms tend to disappear, but in the size range considered here the cations continue to show low-coordinated adatoms up to N=21. To further analyze this issue, we have evaluated excess Bader charges, defined as Qi,exc =Qi−Qtot/N, for all the clusters displaying adatoms. Here, Qtot is the total cluster charge, Nthe number of atoms, and Qithe Bader charge of atom i. Obviously for neutral clusters the excess charges are the same as the Bader charges, but with this definition we can compare neutral and charged clusters on an equal foot as the excess charges always add up to zero. In other words, the excess charge of atom iquantifies the deviation of its atomic charge with respect to the average charge that each atom would have in case the total charge was equally shared by all atoms. The complete set of Bader charges is included as an additional column in the xyz-coordinate files provided in the ESI. Here we just discuss specific examples that suffice to explain trends. The analysis reveals that the positive charge of cations tends to accumulate preferentially on the adatoms. A very revealing example is Cd+ 5, which has as many as four adatoms, each of them singly coordinated to the central atom (the bond distance is 2.98 ˚ A, while the distance between two adatoms is 4.87 ˚ A). Here the central atom has Qexc =−0.05e, while each adatom has Qexc = +0.0125e. The available electrons tend to concentrate around the central atom to stabilize the four radial bonds, while the adatoms themselves carry a substantial excess positive charge, which favors a long distance between them. As another example, Cd+ 12 has two singly-coordinated adatoms, each with Qexc = +0.067e. There is also one threefold coordinated adatom with Qexc = +0.034e. The rest of atoms (all of them more coordinated) have much more similar charges around an average of Qexc =−0.019e, so once more a large fraction of the total positive charge is located on the adatoms, the excess charge increasing with decreasing coordination number. Finally, the only adatom in Cd+ 21 has Qexc = +0.075e; for this cluster there is a central (or core) atom, which has Qexc =−0.103e, and the rest of charges fluctuate much less around an average value of Qexc = +0.0015e. This last example reveals a new general trend for all the studied cations: if there is a central atom, it tends to be negatively charged, while the shell tends to be positively charged. Once more the adatom takes the largest fraction of the total positive charge. It is interesting to compare Cd+ 5with Cd− 5, the only anion displaying a clear adatom. Its GM structure is indeed similar to Cd+ 5, but now three out of the four external atoms manage to bind with each other, and just a single atom remains singly-coordinated. This agrees with the general view that anions are less electron deficient, and the increased number 4|1–25 of electrons can glue more atoms together. Now the central atom has Qexc = +0.071e, i.e. of the opposite sign as compared to cations. Each of the three basal atoms (equivalent under C3vsymmetry) have Qexc =−0.014e, while the singlycoordinated adatom has Qexc =−0.028e. The total charge (now negative) tends to concentrate more on the shell atoms so that the central atom has a positive excess charge (even if its absolute Bader charge is, of course, negative). So what is common to both cations and anions is that the total charge (either positive or negative) becomes more concentrated on the shell in order to minimize the Coulomb self-repulsion energy, and that the excess charge of a shell atom is larger in absolute value the lower its coordination number. As a consequence, an interesting ‘asymmetry” results that differenciates cations from anions: in cations, the available valence electrons concentrate on the central atom to provide four strong radial bonds while sacrifying the tangential bonds, which results in four adatoms; the anion has two more valence electrons than the cation, and they concentrate more on the shell so they can create bonds between shell atoms, in this case only between a subset of the four shell atoms, suggesting that some electron deficiency remains even in the anion pentamer. The symmetry breaking from perfect Tddown to C3vsubgroup can then be tentatively assigned to an electronic frustration effect, in which there are not enough electrons to bind together the four atoms in the shell. In between cations and anions we have the neutral clusters such as Cd21. Now the central atom has Qexc =−0.073eand the adatom Qexc = +0.046e, with the rest of atoms showing a much smaller fluctuation around Qexc = +0.0014e. In other words, there seems to be a tendency, even in neutral systems, towards a negative charge for the core atom and a positive excess charge on the shell atoms, which is most important again for the adatom. This trend is re-inforced in cations, but for anions it enters in competition with their natural tendency to concentrate the total negative charge on the shell, probably explaining why the coordination-dependence of charges is less marked in anions. We shall compare now our results with those of previous computational studies. Initial works on Cd clusters either assumed or predicted an icosahedral growth for sizes around 13 atoms, while our results clearly shows that icosahedra have very low stability in this size range. For example,3estimated the ionization potential of Cadmium clusters using a simple Tight Binding model, without a previous search for structures. They assumed clusters up to 13 atoms to follow an icosahedral growth pattern, while larger clusters were modeled as fragments of a fcc crystal. Their estimation of the ionization potential is hence not trustworthy as it is based on wrong structures. Yonezawa et al.4employed simulated annealing to search for GM structures up to 20 atoms; the simulations were based on DFT but considering only 2 electrons in the active valence set and freezing the clusters from 5000 K to 0 K in less than 1 ps. The extremely fast cooling rate and the more approximate nature of the calculations may explain their finding of an icosahedron for N=13 or their missing of the TTP structure for N=9. Similarly, Flad et al.5reported first-principles calculations aimed to determine the relative weight of covalent and Van de Waals contributions to bonding in Cadmium clusters up to 6 atoms. However, they assumed Cd6structure to be an octahedron instead of performing a detailed structural search. As a consequence, the reliability of their results is also in doubt. Zhao et al.6and Michaelian et al.7employed an EP/DFT scheme with a Gupta potential, so their strategy is in principle identical to ours, yet the results are quite different: for example, Zhao et al. predict an icosahedral growth in the size range N=13 −17, and Michaelian et al. report a GM structure for Cd13 which is not based on the tri-capped trigonal prism unit. As those authors also employed a GGA-DFT level of theory in the final step, the disagreement may only arise from the different Gupta potential employed in their work, which is parametrized to reproduce Cadmium bulk properties instead of properties of small clusters; the two potentials are notably different as shown in our previous work15. In other words, the bulk-fitted potential proposes structures that are very high in energy on the GGA-DFT energy surface. We can explicitly test this assertion by optimizing at the GGA-DFT level the GM structures proposed in those works. As an example, the 13-atom icosahedron is less stable by 0.7 eV than the GM reported here. Johansson and Pyykk¨ o8have suggested that perfect tetrahedral shapes should be very stable for cadmium clusters with N=4,10,20,35 and 40 atoms. Our calculations confirm that the tetrahedron is the GM only for Cd4; for Cd10 it remains true that the tetrahedron is a competitive isomer, but for Cd20 the perfect tetrahedron is around 0.4 eV less stable than the C3 GM structure reported here. Considering the strong similarity between Cd and Zn clusters, we speculate that perfect tetrahedra with 35 and 56 atoms will be even less stable because the optimal number of core atoms for those sizes is expected to be substantially larger than in perfect tetrahedra15. Mu˜ noz et al.9computed neutral Cd clusters with up to 10 atoms using a similar EP/DFT scheme, based however on an extended Lennard-Jones potential which does not have a many-body character and should be less suitable to model atomic interactions in a metallic material like Cadmium. The agreement with their results is very good, though we identify slightly more stable structures for N=7, 8 and 10. Very recently, Kohaut and Springborg10 have reported GM structures for neutral Cd clusters with up to 60 atoms, obtained from global optimization runs based on a genetic algorithm and a tight binding DFT model, followed by re-optimization of structures at the PBE-DFT level. Although the agreement with 1–25 |5 their results is quite fair concerning gross features of the identified structural motifs, we locate different GM structures for many sizes in the range N=8−21. As the authors provided coordinate files for their GM structures, we have explicitly checked that their proposed structures were already present in our Gupta-generated databank, and all of them are found to be less stable than the GM reported here at the PBE+dispersion level. Here the main reason for the discrepancy may be that dispersion was not included in their calculations, while our results in Table I show that dispersion is essential to accurately reproduce cadmium properties in both dimer and bulk limits. Apart from that, the authors do not specify how many isomers were reoptimized at the DFT level for each cluster size, and choosing just a few representatives might also contribute to the observed discrepancy. All in all, we can conclude that the EP/DFT scheme employed is quite advantageous when used properly. Overall, our results seem to outperform most previous theoretical works performed at a GGA-DFT level of theory. At least for those works the theoretical energies from different calculations can be meaningfully compared to decide which method identified the most stable structure. But only a comparison with experimental observables (which we offer in the next sections) can fully assess the accuracy of our putative GM structures. 3.2 Electronic Properties We display in figure 7 the ionization potential (IP), electron affinity (EA) and fundamental gap (defined as [IP(N)− EA(N)]) of neutral cadmium clusters. All the quantities shown are adiabatic, i.e. they include the contribution from the structural relaxation that occurs after the detachment/attachment of one electron. These properties usually present marked variations with size as a result of the delocalized electron shell structure established in most sp−bonded metal clusters. Clusters with a particularly stable electronic shell structure typically display an enhanced stability against dissociation as well, so they can correlate with the abundance maxima observed in experimental mass spectra, to be discussed in the next section. Electronic shell closings are usually associated with local maxima in the IP curve, followed by a sharp drop that reflects the opening of a new electronic shell. In monovalent metals such as the alkalis, the IP additionally shows odd-even alternations, but Cd being divalent, all CdNneutral clusters have an even electron count. The IP curve is in general agreement with the expected electron shell closings, which occur for Ne=8,18,20,34,40 electrons in a spherical jellium model30 and also for Ne=14,30 in spheroidal31 and ultimate32 jellium models due to shape distortions. The IP shows local maxima for N=7,9,15,17,20, i.e. for Ne=14,18,30,34,and 40 electrons. It also shows a sharp drop after N=4 (8 electrons), although this occurs in a region of an overall decreasing trend of the IP, and the drop after N=5 is just as big. Cd9and Cd10 are neighbouring sizes where a jellium shell closing is expected, and our results show that the IP for N=10 is not particularly high. There are two factors that explain this observation: first, Cd10 geometry is not spherical at all and in second place, these clusters are hollow so the 18-electron shell closing is more important, as predicted by hollow sphere jellium models33. Clusters with a closed electronic shell can additionally display a local minimum in the EA curve. If so, the high IP and low EA signatures will be reflected in local maxima in the fundamental gap. Only Cd17 and Cd20 satisfy both conditions so they show the most pronounced maxima in the fundamental gap. Otherwise the gap curve is quite smooth and suggests a steady closing of the gap as the sp−hybridization increases and clusters approach a more metallic behavior. An additional EA minimum is observed for Cd10 which is the only signature of the electron stability expected for Ne=20. We can analyze also the electronic stability of anionic clusters, if we remember that EA of neutral clusters corresponds to the IP of anions. In fact, marked local maxima in the EA curve identify clusters which would like to gain one electron, so they are expected to be very stable as anions. Marked maxima in the EA are obtained for N=9 and 19, so those two anions display an enhanced electronic stability. Similarly, the IP of neutral clusters corresponds to the EA of cations. In figure 7 we can see a marked local minimum at N=18, that suggests an enhanced electronic stability for Cd+ 18. Experimental results for the IP of Cadmium clusters are available in the literature34. The authors analysed the evolution of metallicity from small clusters to the bulk solid using the IP as indicator. They concluded that cadmium clusters start exhibiting a typically metallic trend only for N>20. For sizes smaller than 20 atoms, they found a non-metallic behaviour comprising 2 different regions. Clusters with N=1−10 show a sharp IP decrease from the atomic value of 9 eV to 5.9eV. From 10 to 20 atoms, the IP decreases at a much slower pace from 5.9eV to 5.5eV . Our results agree with the experimental trends. In particular, we reproduce the 2 size intervals with different slopes, with quantitative accuracy. For example, IP experimental values for N=9 and N=10 (6.3 eV and 5.9 eV, respectively) accord quite well with our computed values. Thus, we get a further confirmation of the quality of our results. 3.3 Cluster Stabilities Figures 8-10 display different stability measures for neutral, cation and anion clusters, respectively. The first one is the binding energy per atom Ecoh =E1−EN N, which quantifies the total internal energy content, being therefore a global 6|1–25 stability measure. Secondly, we show the evaporation energy (Eevap = [E1+EN−1]−EN], which provides a more local measure of stability. Finally, the second energy difference (∆2=EN−1+EN+1−2EN=Eevap(N)−Eevap(N+1)) is shown. This last parameter compares the stability of a cluster against the average stability of its two neighbouring sizes, and is considered a more accurate indicator of experimental abundances in mass spectra. Enhanced cluster abundances in mass spectra can usually be explained in terms of two different factors. Some clusters might have an enhanced stability due to a geometric shell closing: these structures are particularly compact, without surface defects such as vacancies, and typically with a high average coordination number. Adding an atom to these structures opens a new geometric shell, the adatom being barely coordinated and so easily removed. On the other hand, clusters of metallic elements can be very stable due to an electronic shell closing. This situation occurs when the cluster behaves as a superatom with a delocalized electronic density describable by jellium models. For neutral clusters, our results predict an enhanced stability for N=4,9,10,15,17,20. As all the magic sizes coincide with expected electron shell closings, we conclude that electronic factors dominate the stabilities of Cd clusters in the small size regime. Nevertheless, geometric factors clearly contribute for some sizes: Cd4is a compact regular tetrahedron, with shorter bond lengths as compared to its neighboring sizes, resulting in an enhanced stability even if a clear IP maximum was not detected at this size; Cd9may be considered a doubly-magic cluster, displaying both electronic and geometric shell closings. In fact, the TTP unit is also very compact and with high average coordination, compared to Cd10 which has a clear adatom on top of the TTP unit. The evaporation energy decreases for Cd10 precisely because of such adatom. The unexpectedly high ∆2value for N=10 can thus be explained only by considering the marked instability of Cd11 coming from both electronic and geometrical factors. The anionic clusters show clear enhanced stabilities for N=3,4,9,15 and 19 (figure 9). Cd− 9and Cd− 19 owe their stability mainly to electronic factors as shown in figure 7. These clusters miss just one electron to complete their electronic shells, so electronic effects are expected to be quite important. Cd− 4and Cd− 15, on the contrary, must owe their stability to geometric effects, since they do not have a particularly high ionization potential. In effect, the Cd− 4tetrahedron is much more compact than either Cd− 3(which is 2-dimensional) or Cd− 5(which has a weakly bonded dangling atom); the same happens for Cd− 15, obtained by glueing together two compact TTP units, while Cd− 16 has an adatom with low coordination. Finally, the high stability of Cd− 3is driven both by the proximity of the 8-electron shell closing and the instability of the negatively-charged dimer. Cationic clusters, displayed in figure 10, are specially stable at sizes N=3,10,18,20, most of which can be explained by electronic factors. For example, figure 7 shows that the electron affinities of both Cd+ 9and Cd9are local maxima, demonstrating that Cd+ 10 is electronically more stable than Cd+ 9. Additionally, considering the low EA value of neutral Cd10, it becomes clear that Cd+ 10 is electronically more stable than Cd+ 11 as well. On the other hand, Cd+ 11 has a dangling adatom, easy to disassociate, promoting also for geometric reasons a higher population for Cd+ 10. Cd+ 20 can be explained in similar terms: on one hand, the EA of Cd19 is a marked maximum, making Cd+ 20 electronically more stable than Cd+ 19; on the other hand, Cd20 exhibits a marked minimum in the EA, so adding 2 extra electrons to Cd+ 20 is not so favourable. Furthermore, Cd+ 21 has a low coordinated adatom as well, reinforcing Cd+ 20 population. Cd+ 18 has an extra electron with respect to the 34 electron shell closure, with a clear electronic stability since this cation has the smallest EA of all. The GM structure of Cd+ 18 is also a geometric shell closing with high D5hsymmetry, so this cluster can be considered as doubly-magic. Finally, Cd+ 3shows a large stability not easily explained in terms of either electronic or geometric reasons. This size is far from an electronic shell closure, and its linear configuration does not suggest an enhanced geometric stability, implying that arguments valid for metallic clusters can not be extended to such small molecules. Once the energetic stability of the computed clusters had been properly analysed, we shall compare it with experimental results in the following section. 4 Interpretation of mass spectrometry experiments We found in the literature several independent measurements of mass spectra for cadmium clusters. For example, Katkuse and coworkers11 reported mass spectra of Cd+ Ncations with N=6−75, obtained using the Secondary Ion Mass Spectrometry (SIMS) technique. In this type of sputtering source, a highly energetic beam of Xe+ions is used to bombard a Cd sample, resulting in the ejection of hot charged clusters which have enough internal energy to evaporate atoms. The abundances measured in the mass spectrometer are established through those evaporation events and so clearly reflect the stabilities of cluster cations. For the range of sizes we are interested in, they observed enhanced abundances at the magic sizes N=10,15,18,20, which coincide very well with the local maxima in the theoretical evaporation energy curve. The experiments predict a high abundance also for N=6, which is in good agreement with theory as well, but unfortunately they could not perform reliable measurements for N<6 so 1–25 |7 a local maximum for N=6 is not completely demonstrated. The only slight disagreement is that our ∆2results predict a similar abundance for cations with 15 and 16 atoms, while no enhanced abundance is detected for N=16 in the experiments. Other than that, the agreement with experiment is very satisfactory. Two other experiments have tried to determine the abundances of neutral clusters. Ruppel and Rademann12 seeded the metal vapor of a high-temperature oven into a stream of noble gas which is then adiabatically expanded to form a supersonic molecular beam. Here the evaporation events occur on neutral clusters and so the populations established at this stage of the experiment should reflect the relative stabilities of neutral clusters. Photoionization of the beam is then employed to generate cations that can be mass selected in a time- of-flight spectrometer. For the photoionization process to be efficient, however, higher-than-threshold conditions were employed, so the resulting cation contains an excess internal energy (not quantified in the experiment) that can promote additional evaporation events during the flight of the cations towards the detector. The measured abundances will reflect the stabilities of neutral clusters only under the assumption that no evaporation occurs during the time of flight; otherwise they will reflect the stabilities of cluster cations, or something intermediate between the cation and neutral stabilities if equilibrium populations can not be established during the time of flight. Theoretical predictions can be very helpful to elucidate the experimental results in cases like this. Ruppel and Rademann report magic sizes at N=6,10,15,18,20 which perfectly coincide with the theoretical magic numbers of Cd+ N. Neutral clusters should show a clear maximum at N=17 instead of N=18, a much higher population for N=9, and finally a magic number at N=4 instead of N=6. Our calculations thus suggest that the mass spectra reported by Ruppel and Rademann correspond to cation and not neutral clusters. Additionally, their results coincide perfectly with the previous experiments by Katakuse for cations. More recently, Diederich et al.13 reported an improvement in the previous experimental setup, wherein the metal vapor is now seeded into a stream of ultracold big helium droplets. Each Cd atom picked up by the droplet travels towards the droplet center where the growth process occurs. The excess energy after addition of a new atom to the cluster is quickly redistributed over the whole system, and the droplet cools by evaporation of helium atoms, so each growth step occurs in an ultracold environment and finally produces cold neutral clusters. Photoionization is then employed to allow mass selection of the positively charged clusters. The authors state that no attached helium remains in the detected cations except for the monomer, so the excess energy introduced by photoionization is certainly enough to evaporate the whole droplet. It is not clear if it is enough to induce further evaporation during the flight of cations to the detector, and the authors explicitly state they have no means to rule out this possibility. The observation of some attached helium atoms to the monomer at least suggests that evaporation during the time of flight may be less probable the smaller the cluster size. Diederich et al. report magic sizes at N=4,6,10,15,18,20, i.e. essentially the same results as in the two previous experiments except for the additional magic number for the tetramer. Once more, our results suggest that even in this improved setup the measured abundances are representative of cation stabilities. The most convincing argument in our opinion is that for neutral clusters one should observe high abundances at N=9 and N=17, both of which are prominent magic numbers with a closed electronic shell. After photoionization, however, both Cd+ 9and Cd+ 17 have very low stability as shown in figure 10, so we expect their population in the cation beam to be easily depleted by evaporation. Similarly, Cd+ 19 has a local minimum in the evaporation energy and also in the ∆2 curve, so it will readily increase the population of Cd+ 18 during the flight time. However, according to our calculations the enhanced abundance observed at N=4 can only come from the original neutral populations inside the droplets, suggesting Diederich et al. approach succeeded in measuring neutral abundances only for the smallest clusters. Additional support for this claim is that the previous experiment by Ruppel and Rademann12 does not show any signature at N=4. Our results suggest that measuring neutral abundances through photoionization mass spectrometry is very difficult to accomplish, and we hope they can stimulate more experimental research on the subject. It seems that using even bigger helium droplets might allow to determine neutral stabilities up to bigger cluster sizes. Another obvious alternative would be to employ closer-to-threshold ionization, which would demand for more efficient detectors as then the cluster ion signal would be weaker. The abundances of N=9 and N=17 could be taken as good indicators in this process: they should steadily increase as the experiment is improved according to our calculations. 5 Emergence of Metallic Behaviour Strictly speaking, a metallic state involves a finite density of states at the Fermi energy. Obviously, finite systems always display a finite HOMO-LUMO gap and the closing of that gap can not be expected to complete until much bigger cluster sizes than those considered in this paper. But the small size regime is nonetheless interesting as it reveals the initial stages of a gradual evolution towards a metallic state. Here we will analyze the progression of sp−hybridization, electron delocalization, and jellium shell structure in the computed cadmium clusters, properties that can all be associated with a metallic-like bonding pattern. A similar analysis has been re- 8|1–25 cently reported for zinc clusters17 and proved to be successful in rationalizing the peculiar bonding properties found in that system. We will examine the electronic density of states (EDOS) as well as real-space visualizations of cluster orbitals. Figure 11 shows the eigenvalue spectra of neutral clusters for selected sizes. The EDOS of Cd4is fully compatible with a 1S21P6jellium-like electronic configuration, and shows a wide HOMO-LUMO gap of more than 3 eV. The high symmetry (Td) of this cluster does not lift the degeneracy of the 1P shell. The Tdisomer of Cd10 (right top panel) similarly displays a closed-shell electronic configuration (1S21P61D102S2) and the widest HOMO-LUMO gap for all isomers with this size. As a curiosity, notice that it additionally shows an accidental near-degeneracy between the 1Dand 2Ssuperatomic orbitals, and a negligible crystal-field splitting of the 1Dshell. These two structures are the most symmetric ones and can be described accurately in terms of a spherical jellium model. In the size range N=5−9 the 1Dshell is progressively filled. For example, Cd5has an electronic configuration of 1S21P61D2. Its prolate shape induces significant crystal-field splitting of the jellium shells, which opens a sizable HOMOLUMO gap between the stabilized 1D3z2−r2orbital and the rest of te 1Dshell, and makes also the 1Pzmore stable than the degenerate 1Px,1Pyorbitals. Cd9(the TTP unit) has a magic number of 18 electrons and the expected 1S21P61D10 configuration, but its slightly oblate distortion causes again fragmentation of 1Pand 1Dshells, with the 1D3z2−r2orbital being now destabilized and becoming the HOMO of the cluster. The GM for Cd10 with C3vsymmetry has a prolate distortion too, yet its closed shell jellium structure is easily identified. Cd17 and Cd20 also conform to a jellium description under a weak oblate distortion. Finally, in Cd15 the prolate distortion is high enough to open significant gaps at electron counts of Ne=4,14,30, as predicted by the ultimate Jellium model32. We can conclude that cadmium clusters generally show an electron shell structure in agreement with a superatom model of delocalized electrons, as most sp−bonded metallic clusters do. The HOMO-LUMO gap tends to decrease with increasing size since the energetic distance between superatom orbitals shrinks with their filling. Nevertheless, there are certain sizes following a spherical shell closure where the two added electrons do not seem to populate the next jellium orbital, hence those sizes do not display a tipically metallic behaviour. In particular, this happens for Cd11 and Cd21: both clusters display an intriguing EDOS, because their shell structure is very similar to that seen in the immediately previous size, apart from a new populated orbital deep below the Fermi level (identified with ’??’ in the figure as its superatomic label is not clear in advance). As a consequence, the HOMOLUMO gap after the expected shell closing can be even wider than at the shell closing (for example, this gap is bigger for Cd11 than for Cd10). This strange behaviour is reminiscent of that expected when an atomic impurity with highly localized orbitals is present within a metal cluster. Therefore, we observe an apparently insulating behavior for these sizes, as opposed to the expected jellium pattern. We shall get a deeper insight from visualization of the charge densities associated with each orbital. In figures 12- 14, we show 3D isodensity plots as well as 2D projections over some planes of interest. Additionally, we show plots of ∆ρ=ρ−ρ0,ρbeing the self-consistent electron density and ρ0the pro-molecular density (i.e. the superposition of atomic densities at the GM geometry). As a representative of the very small clusters, we show results for Cd4in figure 12. On one hand, they confirm the assignment of jellium labels to molecular orbitals; on the other hand, they demonstrate that sp−hybridization is still very incomplete. In fact, the 1S orbital is clearly a delocalized spherical orbital centered at the cluster center-of-mass, and has a four-center two-electron bonding character with its maximum located at the center of the tetrahedron. The charge density associated with the full 1P6shell is also spherical, although with a more fragmented topology, i.e. with lower connectivity between the four atomic basins (the isosurface chosen corresponds to 70% of the maximum density value). Contour pictures show that the 1Pelectron density has the expected node at the center of the structure, but additional local minima along the edges of the tetrahedron, precisely at the midpoint between each pair of atoms. 1P6orbitals hence show some anti-bonding character, indicating a poorly developed sp−hybridization (in absence of hybridization, a full s2atomic configuration will result in a filled s−band, and the top occupied levels will neccessarily show an anti-bonding character). Finally, the LUMO is clearly the 2S jellium orbital, due to its node in the radial direction and its spherical shape. The very small size of Cd4promotes the stability of 2Sover 1D, since the latter (not shown explicitly) is essentially located in an anti-bonding spatial region far away from the surface of the tetrahedron. A somewhat more quantitative picture can be obtained through projections of the EDOS onto atomic orbitals and from ∆ρ, which shows the combined effect of 1Sand 1Por- bitals. The positive values of ∆ρobserved inside the tetrahedron are mainly due to the bonding character of 1S; a shallow local minimum is however observed at the exact center of mass, reflecting the well-known contraction of atomic orbitals upon bond formation. Similarly, although the 1Pshell has a node at the midpoint of each edge, ∆ρis instead slightly positive at those points. Since the satomic orbitals contract during bond formation, some sp−hybridization must be invoked to explain that non-negative value. In fact, projecting the EDOS over atomic sand porbitals shows that the 1Por- bital has 87% s-character and 13% p-character, and that the total p−contribution is equally distributed into px,py,pxor- 1–25 |9 3 4 5 6 -- D 3h -- T d -- C 3v 10 D 3h 30 C 4v -- C 5v 10 C 2v 789 -- D 3 10 C 2v -- D 4d -- D 3h 10 C 4v 10 11 -- D 4d 10 C 3v -- C s 10 C s 10 C 1 Fig. 5 Putative GM structures, and approximate point group symmetries of Cd− Nanions with N=1−11 atoms. Competitive isomers are also shown, together with their energy difference with respect to the GM, expressed in meV. 16 |1–25 12 -- C s 10 C 3v -- C s 20 C 1 -- C 1 40 C 2 15 17 40 C 2 -- C 2v 30 C s -- C 2v 20 C s -- C s 40 C 2v 18 19 20 30 C s -- C 1 -- C 3 60 C s -- C 3 21 20 C 1 -- C s 4 C 1 20 C 1 16 13 14 Fig. 6 Putative GM structures, and approximate point group symmetries of Cd− Nanions with N=12−21 atoms. Competitive isomers are also shown, together with their energy difference with respect to the GM, expressed in meV. 1–25 |17 10 17 20 4 7 Fig. 7 Electronic properties of CdNclusters. From top to bottom, the adiabatic electron affinity, ionization potential and fundamental gap, are shown as a function of cluster size. Specific sizes discussed in the text are explicitly annotated in the plot. 18 |1–25 E coh (eV) E evap (eV) Fig. 8 Stability indicator of CdNclusters. From top to bottom, the binding energy per atom, the evaporation energy and the second energy difference, are shown as a function of cluster size. Specific sizes discussed in the text are explicitly annotated in the plot. 1–25 |19 4 Fig. 9 Stability indicator of Cd− Nclusters. From top to bottom, the binding energy per atom, the evaporation energy and the second energy difference, are shown as a function of cluster size. Specific sizes discussed in the text are explicitly annotated in the plot. 20 |1–25 18 Fig. 10 Stability indicator of Cd+ Nclusters. From top to bottom, the binding energy per atom, the evaporation energy and the second energy difference, are shown as a function of cluster size. Specific sizes discussed in the text are explicitly annotated in the plot. 1–25 |21 Fig. 11 Electronic density of states of selected neutral CdNclusters (notice two different isomers are shown for Cd10). The EDOS curves have been nornalizaed to the total number of electrons, and the HOMO (vertical dotted line) is shifted to zero energy, in such a way that positive energy values correspond to occupied orbitals (the first negative energy value is the LUMO). Standard superatomic jellium labels are attached to each peak or group of peaks. 22 |1–25 Fig. 12 Isosurface and 2D contour plots of several superatom orbitals and ∆ρfor Cd4. Contour plots employ a rainbow color scale, with reddish tones indicating lower values of the scalar field. In particular, red (purple) regions indicate local zones with depletion (accumulation) of electron charge in the ∆ρcontours. Fig. 13 Same as previous figure but showing results for Cd9. 1–25 |23 Fig. 14 Same as previous figure but showing results for Cd20 and Cd21. 24 |1–25 Fig. 15 TOC figure. Associated text: Putative Global Minimum structures and an analysis of the electronic structure of neutral and charged cadmium clusters are reported to gain insight into the gradual insulator-metal transition in the small-size regime. 1–25 |25