(Meta-)stability and Core-Shell Dynamics of Gold Nanoclusters at Finite Temperature

Gold nanoclusters have been the focus of numerous computational studies but an atomistic understanding of their structural and dynamical properties at finite temperature is far from satisfactory. To address this deficiency, we investigate gold nanoclusters via ab initio molecular dynamics, in a range of sizes where a core-shell morphology is observed. We analyze their structure and dynamics using of state-of-the-art techniques, including unsupervised machine-learning nonlinear dimensionality reduction (sketch-map) for describing the similarities and differences among the range of sampled configurations. We find that, whereas the gold nanoclusters exhibit continuous structural rearrangement, they clearly show persistent motifs: a cationic core of one to five atoms is loosely bound to a shell which typically displays a substructure resulting from the competition between locally spherical vs planar fragments. Besides illuminating the properties of core-shell gold nanoclusters, the present study proposes a set of useful tools for understanding their nature in operando.

T ransition-metal nanoclusters have been the focus of a large number of experimental and computational studies in physics and chemistry because of their expected technological applications in diverse areas, including catalysis, 1 optics, 2 and biomedicine. 3 Their physical and chemical properties are dependent on many characteristics such as size, 4 shape, 5 charge state, 6 ligand effects, 7 and temperature. 8 Transition-metal nanoclusters typically exhibit catalytic activity at high temperature, 9 but gold nanoclusters have recently received special attention because of the discovery of their catalytic activity at low temperature. 10 Experimental and computational studies of both colloidal 11−15 and bare 16−23 gold nanoclusters have achieved important advances in the last decades. Most of the reported ab initio studies have been dedicated to searching for their putative global (energy) minimum configurations (pGMC) in the gas phase at zero temperature. 14 An empirical-potential molecular dynamics (MD) study of Au 75 , Au 146 , and Au 457 showed that their melting temperatures were well below the melting temperature of bulk gold (1090 K, at standard pressure) but also that structural solid-to-solid transitions were present well below their melting temperatures. 24 For instance, for Au 75 , the melting temperature is around 550 K, but around 350 K there is a structural transition from truncated-decahedral to icosahedral.
The present study focuses on the gas-phase gold nanoclusters Au 25 , Au 38 , and Au 40 . Au 25 was chosen because it is expected to be the smallest cluster that has a core, in this case a single Au atom, in the low-lying isomer population. 18 Au 38 and Au 40 were chosen because they are both expected to have tetrahedral Au 4 cores and low-symmetry structures have been proposed as their pGMCs. 22,23 We use ab initio Born− Oppenheimer MD, an important computational technique that allows for the sampling of the configurational space of systems with many atoms, together with yielding reliable information on the dynamics of the system, for example, in terms of timecorrelation functions. When the Born−Oppenheimer MD equations of motion are extended by adding a thermostat, the canonical ensemble (at temperature T > 0 K) is sampled. In the specific case of density-functional theory (DFT)-based MD, the energy of the system and forces between atoms are determined by solving the Kohn−Sham equations. Hence, the ground-state electronic structure is accessible along the simulated trajectories.
DFT-based MD studies on small Au clusters have been performed in the past, aimed at describing the thermal evolution of the ground-state structures and the competition between structural isomers at a given temperature, 25−27 as well as at accessing the infrared (vibrational) spectra for comparison with experiments. 28 We analyze DFT-based MD trajectories of Au 25 , Au 38 , and Au 40 at room and higher temperatures in order to address important questions regarding the structural and dynamical aspects specific to this range of cluster sizes: (i) Is the core−shell morphology persistent at finite temperature? (ii) Do the atoms constituting the core exchange their positions with those constituting the shell? (iii) Does the shell of different isomers have recurrent structural motifs?
We ran DFT-based MD within the Perdew−Burke− Ernzerhof 29 (PBE) approximation for the electronic exchange−correlation functional, corrected for the many-body dispersion (MBD) interactions, 30 in the all-electron fullpotential implementation of the FHI-AIMS package. 31 The time step for the numerical integration of the MD equations of motion was set to 10 fs, and we applied a Gaussian broadening of 1 meV for the occupation at the Fermi level. Canonical sampling was imposed via the stochastic velocity-rescaling thermostat. 32 Further MD details are described in the Supporting Information. These settings have been already adopted in refs 27, 28, and 33 and were shown to yield observable quantities in good agreement with experiments in the size range between 3 and 20 gold atoms. Gold nanoclusters are strongly influenced by relativistic effects, 34,35 as the contraction of the 6s orbital and expansion of the 5d results in a smaller s−d gap and an increased directionality in the Au− Au bonding. In the present study, relativistic effects are included via the "atomic" zeroth-order regular approximation to the Dirac equation. 36 The initial configurations for the MD trajectory were ( Figure  1) the pGMC structures described in ref 19 for Au 38 and Au 40 , both containing a tetrahedral Au 4 core, and a higher-energy isomer for Au 40 (also from ref 19), containing a Au 3 core. Furthermore, to increase diversity in the starting point of the MD trajectories of Au 38 and Au 40 , we also selected the respective pGMC of a potential-energy surface described  , and cutoff distances (solid vertical lines) with respect to the cluster's center of mass defining the core, shell, and protrusion regions for the Au 38 nanocluster at 600 and 300 K. (b) Evolution of the Hirshfeld charges for core and shell atoms and (c) evolution of number of atoms belonging to protrusions (green line), shell (orange line), and core (blue line) regions at 600 and 300 K for Au 38 . Panels d−f show the same quantities as in panels a−c, respectively, but for Au 40 .

The Journal of Physical Chemistry Letters
Letter through the embedded-atom method (EAM), globally optimized via the revised basin-hopping Monte Carlo (RBHMC) algorithm, 37 both containing a Au 5 core. For Au 25 , we adopted the embedded-atom-RBHMC pGMC, which has a one-atom core. More details on the RBHMC search are given in Figure S1 of the Supporting Information. We note that close-packed structures, such as icosahedra and decahedra, which can be important for larger gold nanoclusters, are not relevant (i.e., energetically very unfavored) for the smaller sizes considered in the present Letter. 38 For each starting point, one or more ab initio MD trajectories of 100 ps is run, where the first interval of 50 ps is thermostated at a temperature between 400 and 600 K (see Table S1 for details on the single trajectories) and the last interval of 50 ps is at "room temperature", i.e., 300 K.
As detailed below, the statistical analysis of the structure and dynamics of the studied clusters is performed by means of radial-distribution functions, bond-lifetime analysis, nonlinear dimensionality reduction (in order to represent, by exploiting "structure similarities", the high-dimensional configurational space in two dimensions), and state-of-the-art Boltzmann reweighting (in order to use the MD sampling at all temperatures in the most efficient way).
The (center-of-mass centered) radial-distribution function g(r) is the histogram in which bin i counts how many atoms are found between distance r i−1 and r i from the center of mass of the cluster, averaged over all configurations along a trajectory. Here, r 0 = 0 and the radial increment r i+1 − r i = δr is constant and set to 0.1 Å. For extended systems, the histogram is typically normalized by the volume of the spherical shell between r i+1 and r i , but for finite systems, this normalization is not relevant and omitted here. This descriptor, invariant with respect to rigid rotations and translations of the whole cluster and to permutations of atom identities, gives direct and statistical access to the radially layered structure of the clusters and allows for a parameter-free definition of the core and shell regions. Figure 2a shows representative g(r) plots for one trajectory of Au 38 and one trajectory of Au 40 . At both cluster sizes, three clear regions can be distinguished. Some gold atoms are found in the range 1.0− 2.5 Å. These are the atoms constituting the core. Another group of atoms is found in the range 2.5−5.5 Å. These are the atoms constituting the shell. A third group of atoms appears at larger distances. A minimum in g(r) at r ≈ 5.1 Å (Au 38 ) and r ≈ 5.5 Å (Au 40 ) suggests that atoms contributing to the outer part of the distribution form a well-separated group from the shell atoms. Visual inspection of the cluster structures (e.g., green-colored spheres in Figure 4) reveals that these are lowercoordinated atoms "protruding" from an otherwise triangularlattice shell, the latter typically composed of 6−7 fold coordinated atoms. We therefore name this outer region protrusions. The g(r) plots for all trajectories, classified also by temperature, are shown in Figures S2−S4. For Au 38 and Au 40 , the three-region feature is observed for all trajectories and temperatures, with cutoff radii at 2.5 Å for core−shell and in the range of 5.0−5.5 Å for shell-protrusions. The average numbers of gold atoms in the core/shell/protrusion regions (obtained by integrating the g(r) in panels a and d, which is the same as taking the time average of the trajectories in panels c and f) are, for Au 38 at 300 K (600), 4.0/31.0/3.0 (4.1/30.2/ 3.7); for Au 40 at 300 K (600), 4.0/31.8/4.2 (4.1/32.2/3.6). We notice that for Au 38 , atoms move on average from protrusion to shell region when going from 600 to 300 K, while the opposite trend is observed for Au 40 . The numbers of atoms in the core tends to remain constant. These temperature effects are accompanied by a less clear boundary between shell and protrusion regions for both cluster sizes at 600 K, compared to 300 K, signaled by a less marked minimum of g(r) at r ≈ 5.5 Å. Another observation, from the trajectories in panels c and f, is that both shell and protrusion regions for Au 38 are rather stiff at 300 K while more dynamical at 600 K and for Au 40 at both temperatures. For Au 25 , only a (one-atom) core region and a shell region (the latter in the range of 2.5−5.5 Å, as for the larger clusters) is observed.
We note here that atoms are observed to change region along the MD trajectories, so their assignment to a region is dynamical and based on solely the distance from the center of mass. This flexible definition is crucial for our subsequent analysis. In Figure 2, we also show the time evolution of the Hirshfeld charges 39 of the core and shell+protrusion regions (panels b and e), and of the number of atoms belonging to these regions (panels c and f), for Au 38  The qualitative behavior is remarkably consistent at all cluster sizes and temperatures: the core is always fractionally cationic, while the shell+protrusions is always fractionally anionic.
The number of atoms in the core is very stable. For Au 25 , it is always one, and for Au 38 , it is almost always four, with infrequent fluctuation to five (in the trajectories that started from the structure with a five-atom core, one atom quickly, in a few picoseconds, moves to the shell region). For Au 40 , the behavior is more complex. Trajectories starting with four atoms in the core maintain a four-atom core (with sporadic fluctuations to five) throughout the time evolution, even at the higher temperature (600 K). However, the trajectory starting from a three-atom core (structure "J2" in Figure 1) fluctuates between three and four atoms in the core for the 50 ps at 600 K and remains with a three-atom core during the subsequent 50 ps time interval at 300 K, indicating a kinetic trapping in a metastable state. When starting with a five-atom core (EAM-RBHMC pGMC structure), it takes 500 K and more than 10 ps to see the trigonal-bipyramidal five-atom core transform into a tetrahedral four-atom core. At 600 K, the transformation is faster (fraction of a picosecond, but followed by frequent fluctuation back to a five-atom core), while at 400 K, the fiveatom core remains stable. The subsequent parts of the trajectories at 300 K maintain the number of atoms in the core, as at the higher temperature. In this case, we have also the signature of a kinetic trapping that is easily overcome at 600 K.
The number of atoms belonging to either shell or protrusions regions (for Au 38 and Au 40 ) is quite stable, but with larger fluctuations (±2−3 atoms around the average value along the trajectories). Average sizes of the protrusions regions range from 2 to 7 atoms, depending on the cluster size and the temperature.
To further analyze the nature of the different regions in the clusters, we measured the average bond lifetimes, according to whether the bonded pair belongs to one region or is across regions. To this purpose, we first define a bonded pair (at a given time step) as any pair of atoms at a distance smaller than 3.6 Å. This cutoff value is assigned as the first minimum of the pair-distribution function. It is a concept similar to the g(r) defined above, but here the distance between any pair of atoms is binned. The position of the first minimum (see Figure S6) is relatively independent of cluster size, temperature, and region The Journal of Physical Chemistry Letters Letter to which the atoms belong. Next, if a bond is found at any time step, it is monitored until it is first broken (i.e., until the distance between the two atoms exceeds the threshold): The elapsed time is the bond lifetime, and it is binned considering whether both atoms are in the core, or one in the core and one in the shell, etc. In Figure 3, we show τ, the region-dependent average value of the bond lifetime, at each cluster size and temperature (if more trajectories sampled the same temperature, we averaged over them). We also show the sample standard deviation of each bond-lifetime histogram as an error bar in order to quantify the spread of this quantity under the different conditions. More details on the bond lifetimes, broken down to single trajectories, are given in Table S2. We observe a striking difference in behavior between bonds shared by atoms both belonging to the core, or both belonging to the shell, and bonds across core and shell. The latter are systematically at least 3−4 times shorter lived than the former. The shell−protrusion bonds behave similarly to shell−shell bonds, suggesting that the two regions are structurally separated but dynamically very similar. Protrusion−protrusion bonds were too few to yield robust statistics at all sizes and temperatures. The dynamics of the core−shell bonds suggests a picture of a shell that is loosely bound to the core. Visual inspection of the trajectories of Au 40 and Au 38 confirms that the core atoms behave like a inner cluster "confined", and almost freely rotating, inside the shell. In the case of Au 25 , the single-atom core "rattles" around the center of mass.
Consistently, g(r) for Au25 ( Figure S3) shows that the oneatom core is positioned on average at about 0.25 Å from the center of mass.
In order to describe the statistical distribution of structures at each cluster size, as sampled along the MD trajectories, we adopted an unsupervised machine-learning approach, namely, the nonlinear dimensionality-reduction algorithm termed sketch-map, as introduced in ref 40. First, similarly to ref 41, we describe each configuration along a MD trajectory with a coordination histogram, where the jth bin is the number of atoms that have coordination j, i.e., that have j atoms at a distance within the same threshold defined above for the bond breaking. A graphical example of how a cluster structure is mapped into a coordination histogram is given in Figure S7. The histogram is normalized; that is, each bin is divided by the number of atoms in the cluster. For the clusters considered, we found that the maximum coordination is 12, and each histogram was therefore set to have dimension 13, starting from 0-fold coordinated atoms (i.e., atoms evaporated from the clusters) which were, however, not observed in the present study. This 13-dimensional representation is already a dimensionality reduction from the 3N coordinates for an Natom cluster. Such representation is roto-translational and permutation invariant, similar to g(r). In this preliminary dimensionality reduction, some information is necessarily lost, as two different clusters may accidentally have the same coordination histogram. Nonetheless, the coordination histogram captures essential information regarding the bonding network for the clusters. Next, we define a distance in the coordination-histogram representation, which is simply the Euclidean distance between histograms, considered as vectors. The sketch-map algorithm 40 then embeds, in two dimensions, the high-dimensional points (here, in the coordinationhistogram representation), by representing points that are close (distant) in high dimension as close (distant) in the twodimensional (2-D) map. The axes in the plots do not have any physical meaning and cannot even be labeled explicitly; what matters are the distances. Still, each high-dimensional point is mapped into a specific point on the 2-D map. Therefore, by discretizing the 2-D coordinates, one can count how many points land in each discrete (2-D) bin and color-code the map according to the number of points in each bin. For better statistics, we adopted a Boltzmann reweighting approach, the multistate Bennet acceptance-ratio (MBAR) algorithm, 42 to estimate the probability that a structure falls in a given bin, by using the trajectories at all temperatures, while the shown plot refers to 300 K only. In this way, the above-described kinetictrapping issues at 300 K are mitigated by simultaneous use of trajectories at higher temperatures, where no kinetic trapping is observed. For a converged sampling, the logarithm of these probabilities is proportional to the negative of the free-energy of the given Au N cluster. In our case, however, the samplings are most likely not converged, and the probabilities are therefore only qualitative. The results for Au 38 and Au 40 are shown in Figure 4 (results for Au 25 are in Figure S8). Absolute errors for the same sketch-maps, also evaluated via MBAR, are in Figures S9−S11.
Because the positioning of the data points on the sketchmaps is determined only by their 13-dimensional coordinationhistogram representation, the criterion of the color coding can be varied at fixed sketch-map, to highlight different aspects of the phase-space sampling. In Figures S12−S14, we assign one color in the sketch-maps for data points belonging to the same For Au 38 (Figure 4, left), we observe three rather distinct configurational basins. The structures belonging to the left and center basins have a majority of 7-fold coordinated atoms, while for the right basin, structures have a similarly high population of 6-and 7-fold coordinated atoms. The difference between the left and center basins is that structures in the center basin are more spherically shaped, while structures in the left basin are more distorted, with locally planar portions of the shell. Incidentally, both the RBHMC and "J" initial structures belong to the center basin. Furthermore, structures similar to the J structure are visited during the MD trajectories started from the RBHMC structure, suggesting a rather broad sampling of the configurational space along the MD trajectories. The connections between the three basins are narrow funnels, which means that, during the MD trajectories, while structures diffuse freely inside the basins, they need specific, coordinated, and infrequent structural rearrangements in order to pass from one basin to another. The power of the sketch-map is to reveal this latent order, invisible to human eye in the high-dimensional representation, even though there might not be a simple description in terms of geometrical features of the structures.
For Au 40 (Figure 4, right), we also observe three basins, albeit less well-defined than for Au 38 . The center-bottom basin contains structures with a trigonal-bipyramidal five-atom core and a preference for 6-and 7-fold coordinated atoms. The left basin has structures with a tetrahedral four-atom core, while in the right basins there is a combination of structures with threeand four-atom cores, where the common structural motif (of the whole cluster) is the twisted trigonal pyramid. Structures in both the left-and right basins have a majority of 7-fold coordinated atoms.
For Au 25 (see Figure S6), we found tubular-like structures (always with one atom in the core) with one side distorted, an aspect reminiscent of the known Au 24 − pGMC. 20 In a different basin, Au 25 shows a further distortion of the tubular structure, yielding locally planar motifs.
The sketch-maps in Figure 4 (and Figure S6) report also the position of the pGMC and other local minima. We note that these structures are not located in the main basins, which means that these structures are not in the regions of the configurational space mostly visited by the MD sampling. Keeping in mind the caveats expressed above, we refer to these "darker", most visited regions as low-free-energy regions in the following. To rationalize the finding that low-energy, relaxed structures deviate significantly from the low-free-energy configurations, we selected from the MD trajectories several structures in the lowest-free-energy areas of the sketch-maps and geometrically relaxed them to find their local minimum energy structures. This is shown in Figure 4 (and Figure S6) with differently colored dots (selected structures from the MD sampling) and stars (relaxed structures), where dot−star pairs with the same color depict initial and final structures for one relaxation. In all cases, the relaxed structures are in higher-freeenergy regions compared to the selected initial structures. Interestingly, one of the selected structures from the highprobability regions, labeled as 40d, relaxes into the pGMC, further supporting the observation that high-probability structures at finite temperature are significantly different, based on their coordination histograms, from the pGMC as well as from the other lower-energy local minima. We inspected the change upon relaxation of (i) gyration radius, (ii) number of core−shell bonds, and (iii) the average effective . The selected structures distinguish core (blue), shell (yellow), and protrusion (green) regions. For each selected structure, also the corresponding 13-dimensional coordination-histogram, the input for the sketch-map dimensionality reduction, is shown. The dots mark structures selected along the MD runs, positioned in the areas of highest probability, and the stars locate the corresponding geometrically relaxed structures (dot/star of the same color indicate initial/final structures for a geometry relaxation). For Au 38 , the pink star is the pGMC found by us (see its structure in Figure S15). The pGMC identified in ref 19 is marked by a diamond. For Au 40 , the red star is the pGMC found by us, which is also the pGMC identified in ref 19.

The Journal of Physical Chemistry Letters
Letter coordination number 43,44 (ECN). The ECN is an indicator of the average coordination of atoms in a cluster. It adopts a selfconsistent, parameter-free, smoothly decaying definition of bond and therefore coordination number for each atom. These values are reported in Table S3. As a general trend, the gyration radius decreases and the number of core−shell bonds increases upon relaxation. At the same time, the ECN always increases. This indicates that there is a contraction of the clusters upon relaxation and, in particular, a tightening of the core−shell binding. Together with the bond lifetime analysis (Figure 3), this observation further strengthens the picture of typical gold cluster structures at 300 K and higher temperatures with weak coupling between core and shell and therefore larger configurational freedom (hence lower free energy) than relaxed structures.
An analysis of the dynamics of the single trajectories projected onto the 2-D maps reveals that each trajectory of Au 38 and Au 40 visits along its 50 ps all the main basins in the map (i.e., those with an exemplary structure shown in the plots), swapping between neighboring basins between 10 and 100 times. Much slower is the diffusion for Au 25 , where at 300 K the trajectories never leave the basin they start in, but in all three cases, during the initial 50 ps at higher temperature, the trajectory swaps few times between the two main basins. We can therefore conclude that, although full ergodicity of our trajectories cannot be ensured and indeed we notice signatures of some kinetic trapping, single trajectories do not merely (anharmonically) vibrate in one basin but their diffusion across the 2-D maps is significant.
We note, incidentally, that the relaxation of the low-freeenergy structures let us identify pGMC for Au 25 and Au 38 , which are previously unreported in the literature (pink stars in Figures 4 (left) and S6). The corresponding structures are shown in Figure S15. For Au 40 , the pGMC we found (red star in Figure 4 (right)) corresponds with the pGMC in ref 19.
Finally, we investigated the reason for the stability of threedimensional four-and five-atom cores, while isolated neutral Au 4 and Au 5 are known to be planar. 28 We measured the amount of positive charge needed to stabilize the 3-D structure versus the planar one for both Au 4 and Au 5 . In Figure 5, we show the energy of the 3-D structure, relative to the planar one for both Au 4 and Au 5 , where both structures are geometrically unrelaxed for each charge. Some structures like the neutral tetrahedron or trigonal-bipyramidal are mechanically unstable, i.e., have imaginary vibrational frequencies (see Figures S16 and S17). In Figure S16, we show the relaxation of the structures by fixing their symmetry and in Figure S17 without fixing symmetry. We find that tetrahedral Au 4 becomes more stable than planar Au 4 with a charge of 1.5e, while trigonalbipyramidal Au 5 needs 2.0e to become more stable than planar Au 5 . These charge values, even without symmetry constraints, are much larger than the small positive charges observed for the core atoms. Therefore, charge alone is not sufficient to stabilize the 3-D structures, leaving confinement and bonding interactions with the shell as likely reasons for the stabilization.
In summary, we have presented a theoretical study of the stability and dynamics of Au 25 , Au 38 , and Au 40 nanoclusters whose configurational space is sampled via Born−Oppenheimer molecular dynamics at various temperatures. We showed that in the considered range of temperatures (between 300 and 600 K) the gold clusters are not amorphous. Instead, they exhibit a dynamical core−shell structure, with a cationic core (single Au for Au 25 and mainly tetrahedral Au 4 for Au 38 and Au 40 ) loosely bound to an outer anionic flexible shell. We also identify, as a persistent feature, a substructure in the cluster shells-one triangular-lattice region with 6-and 7-fold coordinated atoms and an outer region with fewer lowercoordinated atoms that we term protrusions. By using sketch-mapstwo-dimensional embeddings of the coordinationhistogram representation of the clustersand Boltzmann reweighting, we built approximate free-energy landscapes, at room temperature, of the studied nanoclusters. These show that the global and local minimum energy structures differ significantly from the structures in low free-energy areas, i.e., the most probable at finite temperature. Obtained by using a suite of statistical tools that can be applied readily to similar systems, our results are able to highlight some typical characteristics of the finite-temperature populations of Au 25 , Au 38 , and Au 40 , providing insights into the stability and dynamics of their core−shell regions. Our results (especially those for Au 38 and Au 40 ) suggest a possible mechanism of (equilibrium) cluster growth. If gold atoms are added to a nanocluster in the size range that we have considered, they will tend to remain in the shell and protrusion regions until a critical value is reached when one of the gold atoms stably transfers to the core.

The Journal of Physical Chemistry Letters
Letter