Impurity diffusion in magic-size icosahedral clusters

Bimetallic nanoparticles, i.e., nanoalloys, have recently received considerable attention due to their possible applications in several different fields, such as optics, catalysis, energy storage, biosensing, and nanomedicine.^1–3 The physical and chemical properties of nanoalloys strongly depend on their chemical ordering, i.e., the pattern in which the two elements are arranged within the nanoparticle. The controlled production of different chemical ordering types is, therefore, required in sight of possible technological applications, and a clear knowledge of the thermodynamic stability of the produced phases is desirable as well. Different out-of-equilibrium chemical ordering types can be obtained by specifically tailored synthesis techniques⁴ or as the result of kinetic trapping phenomena taking place during the nanoalloys formation.^5,6 Atoms of the two species are then expected to rearrange within the volume of the nanoparticle, until reaching the most favorable configuration. Transformations of the chemical ordering have, indeed, been observed during annealing experiments for several bimetallic systems,^7–14 as well as in molecular dynamics (MD)^15–17 and Monte Carlo simulations.^18,19

The underlying process of chemical ordering transformations is the diffusion of atoms within the nanoalloy. Understanding the diffusion mechanisms at the atomic level is therefore of great importance as it can provide a substantial contribution to the study of the thermodynamic behavior of these systems. In recent years, some experimental studies have, indeed, focused on this point, analyzing the thermally induced evolution from core–shell to intermixed ordering at atomic resolution and in real time.^20–23 

In general, diffusion is enhanced at the nanoscale. Diffusion processes are typically activated at lower temperatures than those needed for the correspondent bulk system, and much faster rates are measured.^24–27 A significant increase in the diffusion rate when decreasing the nanoparticle size has been observed in both experiments²⁸ and molecular dynamics simulations.^29–31 Regarding the diffusion mechanisms identified so far, the role of the surface is always of primary importance. In most cases, the whole process begins with the disordering of the surface layer: surface atoms diffuse easily, and the motion is eventually transmitted to the atoms in the inner layers, therefore allowing for the rearrangement of the two atomic species.^29,30,32 A second mechanism is possible, which requires the presence of vacancies within the cluster volume; in this case, atoms can diffuse by jumping to unoccupied neighboring sites.^31–34 Again, vacancies first appear on the surface, where defects are created more easily, and then migrate inward through subsequent filling steps;^31,32,34 in Ref. 31, it is pointed out that in large nanoparticles, the vacancy-mediated process is restricted to the surface region, whereas in the inner part, where no vacancies are present, diffusion takes place with the same mechanisms observed for bulk systems.

In this article, we focus on the very fundamental process of the diffusion in nanoalloys, i.e., the displacement of one atom of species A within a volume occupied by atoms of species B. To this end, we use molecular dynamics techniques to simulate the diffusion of a single-atom impurity of either Ag or Au within otherwise pure nanoparticles of Cu and Co, respectively. We consider magic-size icosahedral clusters made of 147 and 309 atoms, corresponding to sizes in the range 1.8–2.2 nm. Non-crystalline icosahedral structures are frequently observed at the nanoscale, and perfect icosahedra have been proved to be particularly stable for several mono- and bi-metallic systems in the considered size range.^35,36

Although diffusion processes are enhanced at the nanoscale, we find that the diffusion of the impurity is not activated at room temperature in a time affordable to our simulations. To further accelerate the evolution, we employ two different simulation techniques, i.e., standard molecular dynamics performed at higher temperatures and metadynamics,³⁷ a well-known molecular dynamics-based technique that is able to accelerate rare events by acting on the free energy profile of the system. Metadynamics is widely used in the computational study of several condensed matter systems. In the field of nanoparticles and nanoalloys, it has been employed to investigate structural transformations of mono- and bi-metallic systems,^38–41 coalescence,⁴² and liquid-to-solid transitions.⁴³ 

Our simulations show that the displacement of the impurity within the icosahedral matrix takes place through some unexpected diffusion mechanims, in which vacancies are created in a completely different way compared to what is described in the literature. We demonstrate that this new type of diffusion process is due to the peculiar properties of the icosahedral structure.

The interaction between atoms in the nanoparticles is modeled by an atomistic potential, developed by Gupta⁴⁴ and by Rosato et al.⁴⁵ on the basis of the second-moment approximation to the tight-binding model.⁴⁶ The form and details of the potential are given in Ref. 17, whereas parameters for AuCo and AgCu are taken from Refs. 47–49 and Refs. 50 and 51, respectively. Structure and chemical ordering of AuCo and AgCu nanoalloys predicted by this model are in good agreement with both density functional theory (DFT) calculations^36,52–55 and experimental results.^11,56–59

Molecular dynamics (MD) simulations are performed using our own code. The equations of motion are solved by the velocity Verlet algorithm with a time step of 5 fs, with the addition of an Andersen thermostat that keeps the temperature constant during the simulation.^60,61 For each system and initial configuration, we consider different temperatures in the range 300–800 K, making at least two independent simulations per temperature. Simulations are stopped once the process under consideration is completed up to a maximum duration of 10 µs.

Metadynamics (MetaD) simulations are performed using our own MD code interfaced with the PLUMED plugin.⁶² In metadynamics, the evolution of the system is affected by the addition of a history-dependent repulsive potential, which is built as a sum of Gaussian functions deposited along the trajectory in the space of a few properly chosen collective variables (CVs).⁶³ The effect of the bias potential is to enhance the exploration of the CV space since the system is forced to move away from CV regions already visited. In addition, the system is able to escape from local energy minima corresponding to metastable configurations, in which it would be trapped for a long time in a standard MD simulation. Thanks to these features, metadynamics can be used to efficiently sample the CV space, to reconstruct free energy profiles, and to accelerate rare events such as transitions between adjacent local energy minima. In this paper, MetaD simulations are used only to accelerate the kinetics of rare events that would not occur in standard MD simulations at room temperature. On the other hand, the reconstruction of free energy profiles has not been attempted because our MetaD simulations were unable to reach the diffusive regime. This is due to two reasons: (i) the initial configurations used in the following are largely energetically unfavorable so that returning back is highly unlikely and (ii) there are no constraints forbidding arbitrary deformations of the clusters so that the metadynamics bias, which forces the system to go away from the already visited regions of the CV space, easily leads to unphysical cluster shapes.

The CVs chosen are the path collective variables (CVs) s and z introduced by Branduardi et al.⁶⁴ combined with the permutation invariant vector (PIV) metric,⁶⁵ as proposed by Pipolo et al.⁶⁶ This scheme seems to be particularly appropriate in our case since, in principle, it directs the evolution toward a desired final configuration, and it is invariant under the permutation of identical atoms. Until now, it has been successfully employed in the computational study of several different systems,^66–68 and here, it is applied for the first time to nanoparticles. Further details on the CVs and on the choice of the CV parameters employed in our simulations can be found in the  Appendix. For each system and initial configuration, we perform at least eight independent simulations. The MetaD parameters are as follows: Gaussian height δ = 0.01 or 0.10 eV, Gaussian width (σ_s, σ_z) = (0.02, 0.1), and deposition stride = 10 or 100 ps. All MetaD simulations are performed at 300 K.

Structures obtained in MD and MetaD simulations are relaxed to find the closest local minimum structure of the energy landscape. All energy and pressure profiles shown in the following are calculated using such locally minimized structures.

The local pressure on a selected atom is calculated from its atomic stress tensor. The atomic stress tensor on atom i is defined as⁶⁹ 

in which a, b = x, y, z, E_i is the atomic energy of the ith atom calculated according to the Gupta potential, V_i is the volume of the ith atom in its bulk crystal structure, and r_ij is the distance between atoms i and j. The eigenvalues of this symmetric 3 × 3 matrix correspond to the principal stress components p_i,k (k = 1, 2, 3). The isotropic pressure on atom i is calculated as the sum of the three principal stress components.⁷⁰ 

We consider perfect icosahedral clusters of sizes 147 and 309. The clusters are completely made of either Cu or Co, with the only exception of a single-atom impurity of Ag and Au, respectively. First, we consider all possible positions of the impurity within the icosahedral matrix, and we calculate the total energy of the cluster in each case. The Mackay icosahedron⁷¹ is composed of concentric atomic layers, i.e., icosahedral shells. In the two icosahedra of size 147 and 309 atoms, four and five shells are completed, respectively. The impurity is placed in the different icosahedral shells, and for each one, all the possible inequivalent sites are evaluated, as shown in Fig. 1. The first shell is made of only one atom, corresponding to the central site of the icosahedron; then, the shell size and the number of inequivalent sites progressively increase (vertex, edge, and facet sites).

The energy of the four considered clusters as a function of the impurity position is shown in Fig. 2. All the systems exhibit the same behavior, i.e., an overall decrease in the total energy as the impurity occupies more and more external shells. Surface sites are most energetically favorable for the placement of the impurity as a consequence of the lower surface energy of Ag and Au compared to Cu and Co.^50,52,72,73 Regarding the inner sites, vertex sites, i.e., sites belonging to the fivefold symmetry axis, are always the least favorable.

This behavior is due to atomic stress effects. In icosahedral nanoparticles, interatomic distances are contracted in the radial direction,⁷⁴ resulting in compressive stress on the atoms placed in the inner sites.¹ The degree of contraction, and thus of compressive stress, is larger in the shells closer to the center,⁷⁵ the pressure on the central atom being particularly high.^75,76 The icosahedral matrix is therefore a non-homogeneous environment, and the placement on the impurity in different sites can strongly affect the total energy of the cluster. In general, we expect the energy to improve as the local pressure on the impurity goes closer to zero. Here, the atomic radius of the impurity is larger compared to the atoms building up the icosahedral matrix, the lattice mismatch being ∼13% for both Ag–Cu and Au–Co. The compressive stress on the impurity inside the cluster is therefore even higher, resulting in a more striking energy difference between different sites.

In Fig. 3, we show the local pressure on the impurity for the four considered clusters as a function of the impurity position. The behavior of the pressure and of the energy (see Fig. 2) look very much alike, showing that atomic stress effects are, indeed, the cause of the energy difference between different inner icosahedral sites. As expected, the highest value of the atomic pressure is found in the center in all cases; then, the pressure on the impurity gradually decreases as it occupies more and more external shells. In each shell, the highest value of the pressure is found when the impurity is placed in the vertex site. This behavior is in agreement with previous stress calculations on icosahedral nanoparticles, in which relatively high pressures are found for all inner sites along the fivefold axis,¹ and it implies that vertex sites are the least favorable for the placement of the impurity inside the cluster. Finally, we note that the pressure on the impurity is negative for all sites on the cluster surface. Negative values of the pressure are associated with tensile stress, which is always present on the nanoparticle surface for all nanoparticle geometric shapes,¹ being particularly strong in the case of the icosahedral geometry.⁷⁷ 

The results shown in Figs. 2 and 3 represent a reference point for the study of the diffusion pathways of the impurity within the icosahedral matrix, which we will present in the following. Essentially, if we initially place the impurity in any of the inner sites, we expect that it will diffuse toward the cluster surface in order to reach the most energetically favorable position. The driving force for the diffusion is the relaxation of the local stress in the volume around the impurity, which progressively decreases as the impurity moves toward the external layers.

The impurity is initially placed in the central site of the icosahedron. This configuration is most energetically unfavorable for all the considered systems, and it is associated with the highest degree of compressive stress on the impurity, as discussed before. The evolution from the initial structure is studied by MD simulations at constant temperature. First, two independent simulations at room temperature (T = 300 K) with a duration of 10 µs are run for each system; in these simulations, no significant changes are observed, neither in the impurity position within the cluster nor in the overall geometry. Perfect icosahedra with the impurity occupying the central site are therefore mestastable configurations, whose lifetime at room temperature exceeds the maximum timescale we are able to reach by MD simulations for the considered nanoparticle sizes (i.e., tens of μs).

To study the impurity diffusion pathways, we need to accelerate the diffusion process so that it can be completed within a time affordable to our simulations. Here, we exploit two different approaches: we increase the temperature of the MD simulations and we perform MetaD simulations at room temperature. In general, rising the temperature modifies the free energy landscape of the system and therefore may produce evolution pathways that are unlikely at lower temperatures. For this reason, we perform MetaD simulations at room temperature, as well. The combination of the two methods allows us to identify a more complete collection of possible evolution pathways, thus providing an overall view of the impurity diffusion process from the center of the icosahedral matrix.

MD simulations are performed at different temperatures, from the room temperature up to a maximum value of 600 K for Au₁Co₁₄₆, 700 K for Ag₁Cu₁₄₆ and Au₁Co₃₀₈, and 800 K for Ag₁Cu₃₀₈. In all cases, the chosen maximum temperature is below the melting point of the corresponding system, as we are interested in the impurity diffusion process within the solid icosahedral matrix. The lowest temperature at which the diffusion is activated on our timescale depends on the system: it is 600 K for Ag₁Cu₁₄₆, 400 K for Au₁Co₁₄₆, 800 K for Ag₁Cu₃₀₈, and 500 K for Au₁Co₃₀₈. We can note that for both sizes, the activation temperature is higher in the case of AgCu, denoting a higher rigidity of the Cu icosahedral matrix compared to Co. In MetaD simulations at room temperature, we always observe the diffusion of the impurity, which begins after few simulation steps. This confirms the efficiency of the employed metadynamics scheme (see the  Appendix for further details on this point).

The diffusion process is analyzed at the atomic level. Interestingly, we find that in both MD and MetaD simulations, the migration of the impurity to the second shell is coupled with the formation of a vacancy in the central site and with an adatom emerging on the surface. The clusters before and after the impurity displacement are shown in Figs. 4(a) and 4(c) and Figs. 4(e) and 4(g), respectively. Here, we show the migration of the Au impurity within the Co icosahedra, but the same process is, indeed, observed in the case of AgCu systems. In Figs. 4(e) and 4(g), we can see that the impurity has moved to one of the sites in the second shell; the central site left by the impurity is not occupied by any Co atom so that a vacancy is created.

Icosahedral structures with a central vacancy have been predicted by global optimization searches for Lennard-Jones clusters,^78–81 for Cu,⁷⁵ Ag,⁷⁵ Au,^82,83 Al,^83–85 Pt,⁸³ and Pb⁸⁶ clusters modeled by atomistic potentials, and also for Ag clusters at the DFT level.⁵³ The presence of the vacancy has been proved to release the strong compressive stress in the central site of the icosahedron, therefore stabilizing the structure.⁷⁵ Icosahedra missing the central atom have also been observed in freezing molecular dynamics simulations of Au nanoparticles.^87,88 On the experimental side, mass spectra obtained for methane clusters synthetized in helium nanodroplets exhibit peaks at icosahedral magic sizes, but the peak corresponding to the three-shell icosahedron is shifted from size 55 to 54;⁸⁹ this is consistent with the presence of a vacancy in the icosahedral structure, which is expected to be at the center on the basis of energy considerations.

Here, we can see the dynamic process of creation of the central vacancy as a consequence of the impurity diffusion. A previous study on single-atom impurity diffusion in nanoparticles has, indeed, described the formation of inner vacancies in icosahedral clusters during the diffusion process.³⁴ In Ref. 34, an Ag atom was initially placed in the central site of an otherwise pure Cu icosahedron of size 147 atoms, and the diffusion of the impurity was studied by means of atomistic MD simulation at T = 700 K. A detailed description of the impurity diffusion process was provided: first, the surface of the cluster becomes partially disordered on one side, and a vacancy appears in the subsurface layer; the vacancy migrates to the second shell and comes into contact with the Ag atom, which is still in the central site; the Ag atom jumps to the empty site in the second shell so that the vacancy is displaced in the center; finally, the structure reverts to the perfect icosahedral geometry, with ordered surface. In the described mechanism, the vacancy first appears in the subsurface layer and then moves inward until it reaches the central site; the whole process involves three steps, taking place at different times. The initial formation of the vacancy requires the partial disordering of the cluster surface, and therefore, many atoms are involved during this vacancy-mediated diffusion process.

Here, we observe a rather different process. The vacancy is created directly in the center of the icosahedron at the same time as the impurity migrates to the second shell. In contrast to what described in Ref. 34, the vacancy does not play the role of promoter of the diffusion of the impurity, but it is a product of the diffusion itself. The whole process can be described as a push-like mechanism, in which the migration of the impurity to the second shell induces a concerted displacement of a row of atoms in the outward direction. The outer atom of the row, i.e., the one belonging to the cluster surface prior to the diffusion, is expelled from the icosahedral matrix and becomes an adatom; in some cases, it is again incorporated into the surface through a well-known surface reconstruction process (rosette reconstruction⁹⁰). In Figs. 4(b) and 4(d) and Figs. 4(f) and 4(h), the clusters before and after the impurity displacement, respectively, are again shown, highlighting all the atoms involved in the diffusion. The remaining Co atoms (displayed as small spheres in Fig. 4) do not change their position during the process; in particular, no surface disordering is observed. We do not observe a preferential direction for the displacement, for example, along a fivefold symmetry axis; instead, moving atoms are often rather misaligned, as in the case displayed in Figs. 4(b) and 4(f). We remark that the processes of Fig. 4 are only representative examples of the diffusion mechanism, which may also take place through more complex displacements. Specifically, it frequently happens that the adatom pops up on the surface in positions that are far away from the emerging subsurface atom. In this case, there, a concerted movement of many surface atoms is observed.

Sometimes, the creation of surface defects precedes the diffusion of the impurity: one or more vertex atoms are detached and diffuse on the cluster surface or are incorporated into it through some reconstruction process. In that case, the collective displacement described before takes place in the direction of a missing vertex so that the empty site on the surface is filled by an atom coming from the subsurface shell. This phenomenon is observed in Ag₁Cu₁₄₆, Ag₁Cu₃₀₈, and Au₁Co₃₀₈ in all MD simulations in which the push-like diffusion mechanism actually takes place, i.e., those performed at relatively high temperatures. Instead, it is never observed in MetaD simulations, suggesting that the creation of surface defects is not activated at room temperature; actually, the detachment of the surface vertex is never observed during the 10 µs-long molecular dynamics simulations at T = 300 K.

Even though the stability of the central vacancy in icosahedral clusters has been extensively discussed in the literature, the one-step creation process described here, in which the vacancy directly appears in the central site, is somewhat unexpected. By comparison, the formation of the vacancy in the subsurface layer as described in Ref. 34 is less impressive since defects frequently appear in the surface of nanoparticles due to the low coordination and enhanced mobility of surface atoms. The formation of defects inside the nanoparticle is much less common. As a matter of fact, the combined process of impurity diffusion and vacancy formation inside the cluster can be easily explained on the basis of stress considerations alone. The atom initiating the whole process is the impurity, which moves to the second shell in order to release its initially strong compressive stress; as a consequence, a vacancy appears in the central site. Normally, one would expect an inner vacancy to be almost instantaneously filled by one of its neighboring atoms. Here, this does not happen due to its favorable position. The relative stability of the vacancy, indeed, prevents atoms in the second shell to move to the central site. Due to the suppression of atomic movement toward the center of the cluster, the collective displacement triggered by the impurity can only proceed outward so that one atom is eventually expelled, as previously described. Actually, each atom involved in such displacement reduces its local pressure; this kind of motion is more favorable than inner circular motions needed to fill the central vacancy, in which atoms move within the same shell or inward, i.e., toward sites with similar or higher degree of compression. The filling of the central vacancy actually takes place in our simulations, but on much longer timescales. This point will be further discussed in Sec. III D.

Alternative diffusion pathways are also observed, in which the nanoparticle geometric structure is significantly deformed; the cluster temporarily takes a liquid-like configuration, at least on one side, and the impurity quickly diffuses toward the surface; the perfect icosahedral geometry is then restored. This disorder-mediated diffusion process was described in Ref. 34. Here, it is observed in the case of the smaller icosahedron, in simulations performed at the highest considered temperatures (700 K for Ag₁Cu₁₄₆ and 600 K for Au₁Co₁₄₆). It is never observed for the icosahedron of size 309 atoms. We do not observe the formation of the vacancy in the subsurface layer, which in Ref. 34 is said to be competitive with the disordered-mediated diffusion process for Ag₁Cu₁₄₆ at 700 K. This may be due to the lack of sufficient statistics at this temperature. However, here, we are mainly interested in the impurity diffusion pathways in which the icosahedral matrix is not deformed, and therefore, we decided not to further investigate the evolution at high temperatures.

The diffusion process is analyzed in terms of energy and atomic pressure on the impurity. In Fig. 5, we show the energy and pressure profiles calculated during the first step of the evolution, i.e., the displacement of the impurity from the first to the second shell. Here, we consider only the simulations in which the cluster keeps the perfect icosahedral geometry, with the only exceptions of adatoms or minor surface reconstructions. Energy and pressure levels corresponding to perfect icosahedral configurations (without the vacancy) are displayed as continuous lines, whereas for configurations with the central vacancy, dashed lines are used. The energy difference between the structures prior and after the diffusion strongly depends on the considered system: for Ag₁Cu₁₄₆, the impurity diffusion slightly increases the total energy [see Fig. 5(a)], whereas in the case of Ag₁Cu₃₀₈, the energy slightly decreases [see Fig. 5(c)]; for both AuCo systems, a significant decrease in the energy is observed [see Figs. 5(e) and 5(g)]. In all cases, the configuration obtained after the diffusion (impurity in the second shell, central vacancy) is less energetically favorable than its perfectly icosahedral counterpart (impurity in the second shell, no central vacancy). On the other hand, we always observe a major improvement of the atomic pressure on the impurity, which decreases by tens of GPa in all cases [see Figs. 5(b), 5(d), 5(f), and 5(h)].

Previously, we have shown that the driving force for the impurity diffusion within the icosahedral matrix is the release of the compressive stress on the impurity itself. As shown in Fig. 5, configurations with the central vacancy are, indeed, effective in this sense, as the pressure on the impurity is much lower compared to the initial value. In the two-step process observed at T > 300 K, the initial surface reconstruction corresponds to a rather large increase in the cluster energy and has small effects on the local pressure on the impurity [see Figs. 5(a)–5(d), 5(g), and 5(h)]; however, if one surface vertex is removed, the atomic row that has to be “pushed” outward by the impurity is shorter by one atom, and the energy barrier of the process is expected to be lower as well. Therefore, at temperatures in which such minor surface reconstruction is activated, it always precedes the migration of the impurity, actually accelerating the diffusion process.

Here, we evaluate the energy barrier for the process of impurity diffusion to the second shell with vacancy and adatom formation. First, we consider the case of Au₁Co₁₄₆, in which the process takes place through a single step at all temperatures. We run 30 independent MD simulations at four different temperatures (450, 500, 550, and 600 K), we evaluate the average time for the occurrence of the event at each temperature, and we make the Arrhenius plot of the rate r_vf of the process of vacancy formation,

where $rvf0$ is the frequency prefactor and E_a is the activation barrier. The activation barrier turns out to be E_a = (0.97 ± 0.03) eV, and $rvf0∼1017s−1$⁠. The reasons for such a large prefactor are likely to be as follows: (i) E_a is an effective barrier that takes into account several similar diffusion processes and (ii) each diffusion process involves the concerted displacement of many atoms.

We tried the same kind of activation barrier evaluation for Au₁Co₃₀₈, Ag₁Cu₁₄₆, and Ag₁Cu₃₀₈. In these cases, one should keep in mind that the process does not take place through a single step, but the displacement of the impurity is preceded by the detachment of a vertex atom, as discussed before. The estimated energy barrier refers to such a combined process; that estimate should reflect the energy difference between the highest saddle point encountered during the whole process and the starting configuration. The highest barrier is that of the rate limiting process, which we observe to be the one in which the impurity moves from the center to the second shell by a collective displacement. Given these caveats, for Au₁Co₃₀₈ and Ag₁Cu₁₄₆, we obtain E_a = (1.16 ± 0.04) eV and (2.1 ± 0.2) eV, respectively. In the case of Ag₁Cu₃₀₈, the temperature range in which it is possible to observe the process in our MD time scale is very limited so that we are unable to produce a reliable estimate of the barrier.

To study the following steps of the diffusion process, we perform new simulations, in which the impurity is initially placed in the second and third shell of the perfect icosahedral matrix, without central vacancies. The initial placement in the third shell is considered only for the larger icosahedron of size 309 atoms; in that case, both inequivalent starting sites are considered (see Fig. 1). Again, in 10 µs-long MD simulations at room temperature, the systems do not display any significant evolution so that MD simulations at higher temperatures and MetaD simulations are needed to study the diffusion process.

For all systems and initial conditions, we observe diffusion pathways similar to the one previously described, in which the migration of the impurity to the outer shell takes place together with the formation of an internal vacancy. The vacancy appears in the site originally occupied by the impurity. The position of the vacancy suggests that its formation is again due to a push-like mechanism initiated by the impurity, which is, indeed, confirmed by the short-timescale analysis of the simulations. Again, small surface reconstructions are sometimes observed before the diffusion, especially at high temperatures. Due to the presence of non-equivalent sites in the third and fourth shell, multiple diffusion pathways are possible. The combination of MD and MetaD simulations allows us to observe and analyze many of them; they are shown in Figs. 6 and 7, respectively, for the diffusion from the second and third icosahedral shell.

The stability of vacancies in the different icosahedral shells was discussed in Ref. 75, in which the authors demonstrated that only central vacancies are favorable, as the vacancy formation energy is positive elsewhere. Energy profiles in Figs. 6 and 7, indeed, show that the process of impurity diffusion, with formation of the vacancy in the second and third shell, always increases significantly the energy of the system. However, as observed for the impurity migration from the central site, the local pressure on the impurity decreases during the process (see again Figs. 6 and 7). In addition, if we compare the two possible configurations (with and without the vacancy) for any given impurity position, we find that the presence of the vacancy is always associated with lower values of the atomic pressure on the impurity. This effect is able to temporarily stabilize the vacancy, which is not immediately filled and, therefore, can actually be observed in our simulations. Sometimes, we observe the subsequent migration of the vacancy from the second shell to the more favorable central site. This process is observed only for the larger icosahedron; it is frequent in MD simulations at high temperatures, but it is observed only in few cases in MetaD simulations at T = 300 K. We never observe the vacancy migration from the third shell toward the inside.

As in the case of impurity migration from the central site, alternative diffusion pathways are observed, in which the nanoparticle structure is greatly deformed. Such pathways are, indeed, much more frequent. In general, the probability of observing structural deformations is higher for the following: higher temperatures compared to room temperature; a smaller icosahedron of size 147 atoms compared to a larger icosahedron of size 309 atoms; AuCo systems compared to AgCu ones; and the initial position of the impurity closer to the cluster surface. Unlike what happens during the impurity diffusion from the central site, deformation-mediated diffusion pathways are now observed even at room temperature; however, the probability of such pathways is rather low in all cases, with the only exception of the diffusion of Au from the second shell on Au₁Co₁₄₆, in which significant structural deformations are actually observed in the majority of our simulations. The vacancy-mediated diffusion described in Ref. 34 is observed as well, but solely during the migration of the impurity from the second to the third shell of the larger icosahedron, and in few cases.

The analysis of MD and MetaD simulations with different initial placements of the impurity gives us a general idea of the whole diffusion process. By putting together the simulation results, we can describe overall diffusion pathways from the center to the subsurface shell of the icosahedral matrix. At each step of the diffusion, a vacancy can appear, diffusion processes coupled to the formation of vacancies being by far the dominant ones at room temperature. Since in each simulation set the impurity is placed in a perfect icosahedron (without inner vacancies), so far, we have actually assumed that the vacancy created at each step is filled before the following step begins. However, due to the already discussed stability of inner vacancies in icosahedral clusters (particularly high in the case of the central vacancy), we cannot exclude the possibility that they survive during the whole diffusion process. This is, indeed, observed in our simulations, as we will briefly discuss in the following.

We initially place the impurity in the central site, and we let the system evolve until the diffusion is completed, i.e., until the impurity reaches the nanoparticle surface. Different diffusion pathways are observed. The central vacancy created during the first step of the diffusion is often filled, while the impurity is in the second shell; however, sometimes, it is still present when the impurity migrates to the third shell, and sometimes, it even survives until the impurity arrives on the surface. Such behavior can be observed for all the considered systems.

In Fig. 8, we show two possible diffusion pathways, obtained for Au₁Co₁₄₆ in two independent MetaD simulations at room temperature. For both pathways, we show the evolution in the space of MetaD CVs, i.e., in the s–z plane, and the behavior of the total energy of the cluster [in Figs. 8(a)–8(d), respectively]. In Figs. 8(e) and 8(f), snapshots of the fundamental steps of the evolution are shown as well, in which we can see the Au impurity gradually migrating from the center to the surface of the cluster. In the first case (left panel of Fig. 8), the central vacancy created during the initial diffusion step is filled, while Au is still in the second shell; as Au migrates to the third shell, a second vacancy appears in the second shell, which is again filled before the subsequent diffusion step takes place. At the end of the evolution, Au reaches the cluster surface, but the lowest energy configuration is not obtained, as some surface defects (rosette reconstruction) are still present. In the second case (right panel of Fig. 8), the central vacancy is maintained during the whole diffusion process so that again the obtained final structure does not correspond to the best configuration for this system. Here, the migration of Au from the second to the third shell takes place without the formation of a further vacancy.

Now, we briefly analyze the progression of the two MetaD simulations in the s–z plane, as displayed in Figs. 8(a) and 8(b). In both cases, we can easily identify four well-separated regions through which the system passes during the evolution. Such regions in the CV space are associated with the placement of the Au impurity in the four different icosahedral shells, and therefore, the sudden jumps between adjacent CV regions correspond to the migration of Au between adjacent shells. The ability to discriminate between configurations in which the impurity occupies different icosahedral shells is a crucial property of the chosen CVs, which, indeed, guarantees the efficiency of our MetaD simulations in accelerating the process of impurity diffusion toward the surface. We note that for any given placement of the impurity, different possible configurations are not so easily distinguishable. Configurations with or without the vacancy, even though rather different regarding the total energy of the cluster, are placed in the same region of the s–z plane so that the vacancy-filling process is not easily identified by solely looking at the behavior of the CVs during the evolution. In short, the CVs used in this work are strongly sensitive to the position of the impurity within the icosahedral matrix but less sensitive to the presence of a vacancy and to other fine details of the cluster structure, such as some surface reconstructions. This point is analyzed in more details in the  Appendix.

In some evolution pathways, we observe the contemporary presence of multiple vacancies inside the cluster, each one created during a different diffusion step. This last scenario happens frequently in MetaD simulations at room temperature; configurations with two or three vacancies are also observed in MD simulations at high temperatures for short times. Multiple vacancies are typically associated with high energy values; therefore, they are expected not to frequently appear during the evolution and to be poorly stable. We believe that the probability of creation of multiple vacancies during the diffusion is somewhat overestimated by metadynamics: probably, in MetaD simulations, the migration of the impurity toward the surface is accelerated to a higher degree than the vacancy-filling process, and therefore, a new vacancy appears before the one created in the previous diffusion step is filled.

In this article, we have described a new mechanism for the diffusion in metallic nanoparticles, in which the displacement of a single-atom impurity between adjacent internal sites is coupled with the formation of a vacancy. In contrast to the vacancy-mediated diffusion process identified so far, the vacancy appears directly in the inner part of the cluster in sites quite far from the surface. The role of the vacancy in the diffusion process is rather different as well. In this case, the vacancy does not promote the diffusion of the impurity, but, on the contrary, it is the result of the diffusion itself, as it appears in the site originally occupied by the impurity prior to the displacement. The combined processes of impurity diffusion and vacancy formation create defects in the surface of the cluster due to collective displacement mechanisms, which involve several atoms. This type of mechanism has been observed for both the considered bimetallic systems and icosahedral sizes and for all initial placements of the impurity within the cluster, being by far the dominant impurity displacement mechanism when the impurity is far from the surface.

This kind of diffusion process may be expected to be a general process for the diffusion of impurities within icosahedral clusters in systems that share the same types of features of AuCo and AgCu; some examples are AgCo, AgNi, AuNi, AgPt, AuPt, and AuRh, in which Au and Ag impurities are larger in size and tend to segregate to the surface. In order to check this point, we calculate the vacancy formation energy, i.e., the energy difference between a configuration in which the impurity is displaced on the surface in an adatom position, leaving a central vacancy, and the starting configuration in which the central site is filled by the impurity. This is done for icosahedral and fcc clusters, for those systems for which Gupta potential models are available in the literature.^91–94 The results are reported in Table I. In all icosahedral clusters, the energy gain in displacing the impurity from the central site to the surface is always large (strongly negative vacancy formation energy). This reinforces the expectation that all binary icosahedra are likely to behave in the same way with regard to the diffusion of the central impurity. We note, however, that for an elemental icosahedron, the vacancy formation energy is always positive or very slightly negative (data not reported in Table I) so that the strong compressive stress on the large impurity is an essential point. Moreover, we performed MetaD simulations for AuRh, which is the system with the smallest energy gain in displacing the impurity on the surface of the icosahedron. Specifically, we performed MetaD simulations at room temperature for Au₁Rh₁₄₆, finding the same mechanism of impurity diffusion as in AuCo and AgCu systems. For what concerns the dependence on cluster size, we expect that the presence of a vacancy at the center of the icosahedron would become more and more energetically favorable as the size increases because the pressure on the central site becomes larger and larger.⁷⁵ This is confirmed by the results in Table I for the icosahedra of sizes 147 and 309. However, it is difficult to predict whether in larger icosahedra the impurity diffusion and vacancy formation will take place by the same elementary mechanism observed for our small icosahedra.

We have shown that the combined process of impurity diffusion and vacancy formation is related to the presence of non-homogeneous compressive stress in the inner part of the icosahedral structure. Therefore, we expect that for other geometric structures (decahedra and crystalline clusters), these collective diffusion processes leading to the formation of internal vacancies are less likely than in icosahedra, as their internal pressure is much smaller and more homogeneous. In fact, we find that the displacement of the impurity from the center to an adatom position in fcc octahedral clusters causes an increase in the energy, instead of the decrease observed for the icosahedra (see Table I). This is a clear evidence that the formation of an internal vacancy in fcc clusters will be much more difficult than in icosahedra.

To simulate the impurity diffusion, we have employed molecular dynamics and metadynamics simulations; the combination of the two techniques has allowed us to identify several possible evolution pathways and, therefore, to achieve a good insight on the impurity diffusion process within the icosahedral matrix. The approach we have employed is, indeed, very general; therefore, it can be applied to the study of more complex phenomena of chemical rearrangement in nanoalloys in order to achieve a more complete understanding of the diffusion at the nanoscale.

D.N. and R.F. acknowledge support from the PRIN 2017 project UTFROM of the Italian MIUR and from the Progetto di Eccellenza of the Physics Department of the University of Genoa. D.N. and R.F. acknowledge networking support from the IRN Nanoalloys of CNRS. F.P. would like to thank Silvio Pipolo for assistance with the use of PIV-based collective variables.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Here, we discuss the form and main properties of the CVs employed in our metadynamics simulations, with particular attention to the choice of the optimal CV parameters for our study.

The path CVs are defined from a reference evolution pathway, consisting of the initial and final states of the studied transition and, optionally, of a discrete set of putative intermediate configurations. Here, we consider only two reference structures, i.e., the chosen starting configuration (A) and the lowest energy structure for the nanoparticle (B). Path CVs for running configuration X are defined in the following way:⁶⁴ 

where $DAX$ and $DBX$ are the distances between the running structure and the two reference structures according to some properly chosen metric. The s-variable tracks the progress of the evolution toward the selected final structure, as it increases as the system moves away from A and approaches B, whereas z is a measure of the distance from the reference pathway since it takes high values when the system explores configurations distant from both A and B. For the parameter λ, we use $λ=2.3/DAB$⁠, with $DAB$ being the distance between the two reference structures, which ensures smooth and well-resolved evolution pathways in the s–z plane.⁶⁶ 

In the framework of the path CVs, a metric for the calculation of distances between different configurations of the system is required. Here, we employ the metric introduced by Gallet and Pietrucci,⁶⁵ which we will briefly describe in the following and we will use the same notation as in Ref. 65, to which reference should be made for a more detailed explanation.

In the metric proposed by Gallet et al., each configuration of a system of size N is associated with a vector v of $12N(N−1)$ components, corresponding to the total number of possible atomic pairs; the vector is built in the following way:

where R_i (R_j) are the Cartesian coordinates of atom i (j) and C_ij is a monotonically decreasing function of the interatomic distance d_ij = |R_i − R_j|. The vector v includes all the elements from the upper triangular part of the symmetric matrix C. These elements are sorted in ascending order within each block characterized by the pairing of the same elements. The sorting operation ensures the invariance of v under the permutation of identical atoms; for this reason, v is called permutation invariant vector (PIV). Finally, distances between configurations are computed as squared Euclidean distances between their PIVs.

The matrix C from which the PIV is built [see Eq. (A2)] is called adjacency matrix; each entry refers to a specific pair of atoms and provides its level of connection, which is here evaluated by a function C(d_ij) of the interatomic distance. This function, also called switching function, monotonically decreases from one (maximum level of connection) to zero (no connection) as the spatial distance between the two atoms increases. Different choices for the switching function are possible; in this work, we have used the decreasing sigmoid function

The parameters r₀, m, and n can be adjusted to select the range of interatomic distances in which the switching function actually decays from one to zero; in this way, it is possible to focus the metric on the range of distances, which are relevant for the studied transition.

The choice of the parameters of the switching function is, indeed, crucial to ensure the efficiency of MetaD simulations. In general, the transition pathway must be well-resolved in the CV space, i.e., the chosen metadynamics CVs must be able to distinguish between the most relevant configurations appearing during the evolution. The path CVs defined in Eq. (A1), together with the choice of the parameter λ already discussed, ensure that the two chosen reference structures are always well separated in the s–z plane; it is easy to verify that s_A ≃ 1.1 and s_B ≃ 1.9, regardless of the details of the chosen metric. However, we need to ensure that the intermediate configurations of the transition, for which we have 1.1 ≲ s_X ≲ 1.9, are sufficiently resolved as well. This is possible if distances between intermediate configurations are reasonably large according to the chosen metric. Strictly speaking, because of the choice of the parameter λ as a function of the distance $DAB$ between the two reference structures, what is actually relevant in the building of the path CVs is the relative distance: if we consider two intermediate configurations X and Y, the quantity we have to evaluate is the ratio $DXY/DAB$⁠. Optimal values of the parameters of the switching function are therefore those that provide large relative distances between intermediate structures.

In this study, the parameters m and n have been fixed to 6 and 3, respectively. Here, we will describe the protocol we have used for the choice of r₀; data refer to the system Au₁Co₃₀₈.

First, a set of putative intermediate configurations have to be selected. Here, in the first reference structure (A), the Au impurity is in the center of the icosahedron, whereas in the second one (B), Au is on the surface on the vertex site; it is therefore reasonable to choose configurations in which Au is placed in the different icosahedral shells, as those considered in Figs. 2(d) and 3(d). We select one configuration per shell, namely, X_k, with k = 1, 2, 3, 4, 5, referring to the number of the shell in which the impurity is placed. When multiple sites are present in the same shell, we always consider the vertex site (see Fig. 1). Of course, we have X₁ ≡ A and X₅ ≡ B.

We consider different values of r₀ in the range 1–12 Å, and for those values, we evaluate the relative distances between all pairs of adjacent structures $DXkXk+1/DAB$⁠, with k = 1, 2, 3, 4. The results are shown in Fig. 9.

Before analyzing the results of Fig. 9, it is useful to make some general remarks on the PIVs of the chosen structures. All configurations in the X_k set exhibit the same overall geometry, as they differ only in the position of the Au impurity within the icosahedral matrix. As a consequence, changes in PIV terms referring to the Co–Co pairs between the different configurations are small. On the other hand, we expect a quite large variation in the Au–Co terms, which, therefore, mostly affects the calculation of pair distances. To rationalize the behavior of the calculated relative distances, we thus have to focus on the effect of the parameter r₀ on the Au–Co PIV terms.

Two different trends are observed. In the case of A–X₂ and X₄–B, the relative distance takes its maximum value for r₀ = 1 Å and then decreases as r₀ is increased [see Figs. 9(a) and 9(d)]. Such behavior is consistent with short-range differences between the two structures of each pair; in particular, the structures can be well-distinguished if we focus on the first neighbor range of interatomic distances. This is clear for the X₄–B configuration pair: Au occupies an inner site when in X₄ and a surface site when in B; therefore, the number of its first neighbors is different in the two configurations. Regarding A–X₁, the number of Au first neighbors does not change since Au in placed in the inner part of the cluster in both cases. However, the values of the first neighbors distance are slightly different: when Au is in the compressed central site, all its twelve first neighbors are 2.45 Å apart, whereas when it is placed in the second shell, the average first neighbors’ distance is 2.53 Å. The difference is, indeed, small, but it can be detected by a properly chosen metric; here, values of r₀ in the range 1–3 Å are able to distinguish the two configurations to an acceptable degree.

Relative distances calculated for X₂–X₃ and X₃–X₄ exhibit a rather different behavior, as shown in Figs. 9(b) and 9(c). Small values are obtained when r₀ is small; then, the distance increases with r₀, and the maximum is reached for r₀ = 8 and 5 Å, respectively. In configurations X₂, X₃, and X₄, Au is placed in inner sites with similar first neighbors distances; therefore, short-range effects poorly affect the distance between structures. In this case, the largest relative distances are obtained for higher values of r₀, which are able to detect the different position of Au within the whole icosahedral matrix.

We calculate the product of all relative distances as a function of r₀; optimal values of the parameter r₀ are those that maximize this quantity, therefore ensuring sufficiently large adjacent pair distances between all the putative intermediate structures. The results are shown in Fig. 10(a), in which we can identify the maximum value in correspondence of r₀ = 4 Å. In Fig. 10(b), we show the results of the same calculation but without taking into account the relative distance between X₄ and B. For this configuration pair, quite large relative distances are obtained for all the considered values of r₀ [see the comparison between Figs. 9(d) and 9(a)–9(c)]; therefore, it is reasonable to focus on the contribution of the other pairs only. In this case, the maximum of the product of the relative distances is shifted to slightly larger values of r₀ and the peak is more broadened; good values of r₀ are in the range 4–7 Å.

Similar results for the optimal values of r₀ are found for all the considered systems, i.e., Ag₁Cu₁₄₆, Ag₁Cu₃₀₈, Au₁Co₁₄₆, and Au₁Co₃₀₈. In our MetaD simulations, we therefore decide to employ r₀ = 4 Å in all cases.

In the case of Au₁Co₃₀₈, the value r₀ = 10 Å is tested as well, but it always leads to non-physical trajectories, in which the cluster disintegrates after few simulation steps. This behavior is due to the inability of the metric to distinguish between different configurations in the early stage of the evolution [see Fig. 9(a)], which are, therefore, all located in a small region of the s–z plane. After some MetaD steps, a high bias potential is deposited in this small area; the resulting forces on the system are huge, and the system is eventually destroyed. The bad results obtained in this simulation set further demonstrate the importance of the correct setting of the CVs parameters for the efficiency of the metadynamics technique.