Multiscale Modeling of Au-Island Ripening on Au ( 100 )

We describe a multiscale modeling hierarchy for the particular case of Au-island ripening on Au(100). Starting at the microscopic scale, density functional theory was used to investigate a limited number of self-diffusion processes on perfect and imperfect Au(100) surfaces. The obtained structural and energetic information served as basis for optimizing a reactive forcefield (here ReaxFF), which afterwards was used to address the mesoscopic scale. Reactive force field simulations were performed to investigate more diffusion possibilities at a lower computational cost but with similar accuracy. Finally, we reached the macroscale by means of kinetic Monte Carlo (kMC) simulations. The reaction rates for the reaction process database used in the kMC simulations were generated using the reactive force field. Using this strategy, we simulated nucleation, aggregation, and fluctuation processes for monoatomic high islands on Au(100) and modeled their equilibrium shape structures. Finally, by calculating the step line tension at different temperatures, we were able to make a direct comparison with available experimental data.


Introduction
Surface diffusion plays a key role in dynamical processes occurring on metallic surfaces which involve physical or chemical displacements of atoms or molecules on a particular surface [1], such as adsorption-desorption, crystal growth, coarsening, chemical reactions, wetting, spreading, or capillarity.Hence, a profound understanding of this phenomenon has implications in a great variety of fields such as electrochemistry, microelectronics, tribology, and corrosion protection [2].Two main pathways can be distinguished for the diffusion of adatoms across a surface.The first one corresponds to the hopping of an adatom between different adsorption sites, whereas in the so-called exchange mechanism, adatoms replace atoms of the surface layer (or fall into vacancies) while surface atoms are simultaneously pushed onto adsorption sites [3].The latter pathway was proposed for Al(001) [4] and Au(100) [5] based on theoretical studies (Müller and Ibach stated that the hopping mechanism could be also competitive on Au(100)) [6].While for planar terraces these processes are rather well understood, [7][8][9][10][11][12][13][14] much less is known about diffusion processes on imperfect surfaces with lower coordinated adsorption sites (e.g., kink-and step-sites or vacancies).However, these types of sites often determine morphological changes of the surface structure, for instance, in Ostwald ripening, metal deposition and dissolution, or island and step fluctuation processes.It should be noted that clean Au(100) shows a quasi-hexagonal reconstruction in the first surface layer [15], however, this reconstruction can be lifted in an electrochemical cell by applying a sufficiently positive potential.Under these conditions, chemisorbed species (e.g., anions) and electric fields are present at the metal interface.Although these parameters are not taken into account in the present work; we reported the influence of these parameters in previous experimental and theoretical studies [16][17][18].Moreover, the lifting and the restoration of the reconstruction are slow processes, so that both types of surfaces can be assumed to be present over a sizable range of potentials.
In order to address diffusion phenomena at surfaces theoretically, different approaches are used, which range from experimentally based, semiempirical molecular dynamics up to first principles studies.In the present work, we describe a first principles-based approach for investigating surface dynamics on Au(100).While improvements in experimental techniques have led to considerable progress in the investigation of surface diffusion, there is still a lack on the experimental data for the self-diffusion on Au(100).The basis for the present study is a recent contribution from our group where the most reasonable processes for self-diffusion on Au(100) were studied using quantum mechanics (density functional theory, DFT).This study specifically focused on terraces and imperfect kink-and step-configurations [19].Since the rate of adatom migration depends mainly on its immediate environment, we concentrated on the nearest neighbor interactions only.Nevertheless, with this restriction, already more than a thousand surface configurations have to be considered [20].Given the high computational expense of DFT calculations, with the present computing resources, one is rather limited in the number of systems that can be studied from first principles.However, we can use the binding energies and self-diffusion processes previously studied by ab initio simulations to construct and optimize a reactive force field (here ReaxFF) that is able to reproduce this behavior.After the mandatory tests related to accuracy and transferability of the force field, ReaxFF has the advantage of being nearly as accurate as QM but with much lower computational demands [21].Using ReaxFF to investigate all possible diffusion processes (considering only first-nearest neighbors), the obtained process rates were used with kinetic Monte Carlo simulations (kMC) for large-scale simulations that allowed us to analyze island nucleation, aggregation, and fluctuation as well as the equilibrium shapes of adatom islands on Au(100). Figure 1 shows the overall computational procedure, where we started at the ab initio level of theory to obtain a reactive molecular dynamics forcefields, which was then used to provide the necessary process rates for macroscopic kinetic Monte Carlo simulations.Similar hierarchical simulation schemes have  been employed by the groups of Goddard, Neurock, and van Santen, for instance [22][23][24][25].

First Principles Calculations
We used density functional theory to investigate different self-diffusion mechanisms in the presence or absence of a variety of surface defects on Au(100), for example, kink-and step-edges.Usually, DFT calculations in the gradient generalized approximations (such as the PBE density functional used in this work) are able to achieve a chemical accuracy of ≤0.1 eV for reactions [13].Since the details of these studies are described in [19], here we only give a brief summary.From these first principle studies, we were able to obtain binding energies of the stable intermediates and transition states and to extract rate constants for the diffusion events.An overview of the different pathways that have been studied is given in Figure 2. On terraced Au(100) surfaces, we found that the lowest energy mechanism was the exchange diffusion pathway from position A to C (E act = 0.55 eV), while hopping diffusion from A to B over a bridge site was less favorable (E act = 0.70 eV).The last evaluated reaction pathway at terraces was the atom hopping over the top position from A to C, being the least favorable with a rather high energy barrier of 1.39 eV.
The presence of step-edges has a tremendous influence on the Au-adatom diffusion, since the energy barrier is even lower for the diffusion along the step-edge (0.38 eV).As expected, diffusion perpendicular to or from the step edge has a much higher barrier of 0.84 eV (detachment) and 0.66 eV (attachment).Comparing the barriers shows that rearrangement of adatoms along the step would be faster than island growth or decay, which would lead to islands that quickly assume their equilibrium shapes before changing their size.
Afterwards, systems involving corners at step-edges were analyzed (see Figure 2).For the inner-corner system, we found that diffusion along long edges is preferred over migration along short edges.This suggests that Au atoms adsorbed at terraces first reach kink-sites or inner-corner positions at islands by exchange diffusion, followed by a hopping diffusion to the step-edges where they rapidly diffuse along the step-edges.As far as the outer-corner systems are concerned, adatoms, which are located directly next to an island corner, have essentially the same binding energy as an adatom at a step-edge.In addition, we found that detachment from the outer-corner was preferred over detachment away from the step-edge.
Finally, we calculated diffusion pathways for forming Au-dimers, which can be considered the first step of island formation, and for extraction of an adatom out of a stepedge forming a step-vacancy.While the latter process has a relatively high energy barrier of 1.07 eV, dimer diffusion is comparable to the case of a single atom.
In addition to reaction barriers, we calculated vibrational frequencies factors for each reaction.This was accomplished by fitting a frequency to the energy profile of each reactant state to within the harmonic approximation.Using these frequencies, we calculated the rates for the different diffusion events at variable temperatures (see Table 1) by means of the transition-state theory, as given by equation ( 1): where k 0 is a prefactor, E act is the barrier for diffusion (activation), k B is the Boltzmann constant, and T is the temperature (all of our mechanisms are first-order kinetic reactions, and so these terms are simply the vibrational frequency of the self-diffusion process in units of Hz).One of the main limitations of this approach is that it does not account for anharmonicites and thermal expansion, which cause deviations at higher temperatures.Moreover, it also ignores memory effects arising from the surface excitations, which change the effective Arrhenius barrier.Nevertheless, the Arrhenius law's validity within the limit of E act /k B T ≥ 4 provides reasonable temperature ranges, for which our results are valid [26].
The obtained diffusion rate constants (Table 1) clearly illustrate the importance of the surface structure on cluster formation and island growth.The fastest rate corresponds to the filling of an inner-corner by a diffusing atom (IV-2 pathway), although the same system also yields the slowest rate with pathway IV-3 (not considering the very unfavorable A → C hopping diffusion at terraces).Concerning the terraces, the diffusion over the top position is practically negligible while the exchange diffusion rate is around 300 times larger than the bridge diffusion rate.Our calculated activation energies for the A → B hopping process and for the A → C exchange process on the terrace are 0.70 and 0.55 eV, respectively.These values are in reasonable agreement with the approximated experimental activation barrier of 0.5 eV for the adatom hopping on Au(100) [27].In addition, they also agree well with previous theoretical studies.For example, Müller and Ibach reported energy barriers of 0.64 eV and 0.60 eV for the hopping and the exchange processes on an Au(100) [6].
Although ab initio studies only allow the investigation of a limited set of diffusion possibilities, the studies presented so far already provide useful insights into diffusion events on Au(100) with important implications for Ostwald ripening and island formation.However, the next step is to account for more possibilities, as described subsequently in the list of possible diffusion events section, by means of larger scale reactive forcefield calculations.Therefore, in the next section we will summarize the generation of the corresponding Auforcefield, which was then used to evaluate the rates for many more diffusion scenarios.

Reactive Forcefield (ReaxFF) Simulations
The next step in the multiscale approach presented in this work corresponds to the generation and application of a reactive forcefield within the ReaxFF framework [28][29][30].ReaxFF is a reactive molecular dynamics method that uses a bond-order-dependent potential energy formulation, similar to the Tersoff [31] or Brenner [32] potentials.These potentials are all based on Pauling's idea of mapping bond distances onto bond orders to enable the determination of different quantum chemical states of a molecular structure [33].
The dependence of the energy contributions on the bond order also means that such energy terms implicitly contain multibody contributions.The Au forcefield employed in the present study contains the following three energy terms (2): where E bond is the energy corresponding to interatomic bonds, E over is a penalty energy that corrects atomic over-coordinations, and E vdW accounts for van der Waals interactions and interatomic repulsions when atoms are too close to each other.The full expressions for these terms as well as the included parameters can be found in more detail in a previous publication from our group [21].
The formulation of the system energy is more sophisticated than those of nonreactive potentials, for example, EAM [34][35][36][37], MEAM [38], UFF [39,40], CHARMM [41], OPLS [42], or AMBER [43,44], and therefore computationally more demanding.Reactive forcefields represent a useful tool for overcoming limitations of these empirical potentials as it allows for the description of chemical reactions, that is, bond formation and dissociation, with almost QM accuracy.In contrast to most other forcefields, the ReaxFF potential used in the present work is completely trained against ab initio DFT calculations.
The goal of this training process is to optimize the parameters of the forcefield such that the DFT energies and structures are reproduced as accurately as possible.In this specific case, first principles DFT-PBE calculations on the equations of state for several gold bulk phases, as well as binding energies and self-diffusion processes on Au(100) surfaces as described in the previous section, were used in the training procedure.The validity and the transferability of the obtained interaction potential was successfully tested for additional diffusion processes, surface free energies of Au low index surfaces, and cohesive energies of molecular Au clusters.Details concerning these tests as well as a first application to study the morphology of larger Au nanoparticles can be found in [21].The purpose was to obtain a large database of diffusion processes on Au(100) under a wide range of defects that one might expect to exist on real surfaces, in order to obtain the diffusion barriers which are needed to run kMC simulations.As already mentioned before, the simulation of all the diffusion processes that could occur is a task that cannot be achieved in a reasonable time using DFT or other ab initio methods.
Figure 3 shows exemplarily the performance of the ReaxFF Au-forcefield after the optimization procedure.Besides simple terrace (bridge) diffusion, diffusion along step edges, around kink sites, and diffusion forming molecular dimers are compared.The comparison shows that the forcefield is indeed a useful tool for calculating all the remaining diffusion processes required to evaluate the database of diffusion rates required for the kMC approach, and thus for bridging the gap between QM and the macroscopic regime.

Kinetic Monte Carlo
Kinetic Monte Carlo (kMC) is a well-known technique for large-scale simulations in various areas of research such as chemical physics or materials design [45][46][47].In order to perform such kMC simulations on the dynamic behavior of two-dimensional gold islands on Au(100), we used the previously described ReaxFF forcefield to generate a list of independent (Markov chain) diffusion rates for all possible diffusion scenarios (taking the first-nearest neighbors into account) [48].Along the kMC simulations, the total rate R (cumulative function) of all diffusion events of a particular system is calculated by means of (3): ( In this function, all possible elementary diffusion processes are taken into account by their individual rates k i .Then, by using a random number, ρ 1 , between zero and one, we assign it to a rate out of all the processes described in the cumulative function (R): where l then specifies the particular process which will be executed.After each diffusion event, the simulation time is propagated by However, for a realistic value of Δt, we have to allow fluctuations, therefore, a second random number, ρ 2 , is chosen between zero and one and afterwards used in (5) in the following way: After executing the particular diffusion process and propagating the time the iterative simulation proceeds by calculating the cumulative function (R) of the new structure, and so forth.As the kMC approach is based on random numbers, reasonable statistics are required to allow for reliable qualitative and quantitative conclusions.

List of Possible Diffusion Events
Given the structure of the Au(100) surface, the number of possible diffusion events, (considering only next-nearest neighbors and neglecting symmetry considerations) is 2 10 = 1024, which is the number of reactions required for the rate catalog.The profile of the considered diffusion pathways is depicted in Figure 4.While the atom of interest diffuses from position 1 to 2, the nearest neighbor (NN) sites (labeled from 3 to 12) are either empty or occupied.
Although the classification described in Figure 4 only considers the next-nearest neighbors, when calculating the diffusion events with ReaxFF and afterwards the rates to be used in the following kMC simulations, a slightly modified system has been employed.For the ReaxFF studies, each system was characterized by a four-layer slab for the underlying Au(100) surface with a lateral extension of 10 × 10 atoms.When occupying the surrounding sites, not only the positions specified by the labels 3-12 have been filled, but always the entire outgoing rows of atoms.Thus, for instance, the diffusion process indicated in Figure 4(a) has been studied with the model shown in Figure 4(b).On the basis of this structure, we used ReaxFF to map the energy profile when moving the diffusion atom (red) successively from position 1 to 2. The energy profile provided the activation energy and vibrational frequencies, which together with transition state theory (TST) provided all diffusion rates to be included in the event database for the kMC studies.

Input Structures
Kinetic Monte Carlo was utilized to study different systems, in which the Au(100) substrate was modeled by a square lattice of 100 × 100 positions.Periodic boundary conditions were imposed in x-and y-directions (parallel to the substrate) and all atoms were assumed to be uncharged.In this work, we analyzed four different initial input configurations, their evolution are being followed at 850 K. Figure 5 shows (a) a random distribution of adatoms, (b) a two-island structure, (c) a spherical island, and (d) a square island.
It should be noted that although during the kMC simulations the systems were always free to evolve, for simplicity, the exchange mechanism at terraces was not included in the rate catalog.

Nucleation and Aggregation
Subsequently, we followed the evolution of the initial input structures (see Figure 5) at 850 K.The first structure (a) starts with a randomly distributed dispersion of adatoms on the Au(100) surface.In Figure 6, we show successive snapshots of the kMC simulation of this structure.Along the structural evolution, particle nucleation is observed, as, after several hundreds of kMC steps, the randomly dispersed adatoms begin forming nucleation seeds.These nucleation seeds afterwards form small islands that successively grow during the ∼4000 ns simulation period.Thus, the process indicated in Figure 6 indeed shows the first stages of island formation.As expected, we also found that the nucleation speed increases with temperature and adatom coverage.This can be explained by the fact that with randomly dispersed adatoms at the beginning, the highest rate corresponds to the free hopping of adatoms on the surface.After several simulation steps, adatoms meet each other forming slowly growing islands.This makes the step-edge diffusion possible and avoids the diffusion away from the step (the latter process has a higher energy barrier), therefore, continuously favoring the growth of the islands.With the second system (Figure 5(b)), we followed the evolution of a pair of disparately sized islands at 850 K (see Figure 7).After several kMC steps, the smaller island starts dissolving, while the detached Au atoms quickly migrate to the larger island where they are assimilated.Hence, this simulation provides insights on how island growth proceeds.

Island Shape (Equilibrium and Fluctuation)
Step-edge fluctuation is the dominating process in determining the dynamic shape of islands on the Au(100) surface.Several diffusion processes contribute to step-edge fluctuation: diffusion along the step-edge, along a kink-site, or along an outer-corner.As explained in the DFT section, the rate for step-edge-diffusion is the highest among these processes, followed by the diffusion along the outer-corner and by the diffusion along kink sites.Figure 8 shows the evolution of an initially spherical island (Figure 5(c)) at 850 K.The snapshots show that after several iterations, the roughness of the island border increases.Qualitatively, this is in good agreement with experimental results obtained by Giesen et al. by means of STM measurements in the case of two-dimensional Cu islands on Cu(100) [14,49].
Interestingly, starting with a square-shaped island and maintaining a temperature of 850 K (Figure 5(d)) showed no significant deviations from the behavior of the spherical island, and hence the corresponding snapshots are not shown.Every migration step of an adatom out of the square island produces a vacancy.This mechanism has the lowest reaction rate among all the considered ones, and therefore this occurs very rarely.This can be explained by the degrees of freedom.As noted already by Ibach [50] step-atoms can only move in one direction and therefore the binding order on the (100) surface is three.This means that three bonds need to be broken along this process, which results in an increased energy barrier for this diffusion process.
Finally, we analyzed the equilibrium shape of the islands for a fixed number of gold adatoms but at different temperatures.To obtain the equilibrium shape of islands at a given temperature, one needs to average over many individual samples [51].Therefore, we analyzed the island shapes during the ripening process starting with an initially random adatom distribution Figure 5(a).Islands already begin to show their equilibrium shape during growth, when the effects of island aggregation are averaged out [49,52,53].We Advances in Physical Chemistry studied the temperature dependence of equilibrium fluctuations of monatomic high islands on Au(100) by analyzing the obtained island shapes at three different temperatures, 500, 750, and 1000 K.As shown in Figure 9, the equilibrium island shape becomes rougher with increasing temperature, which was also suggested by Giesen et al. [49].This roughening was found in several STM-studies on Pb and Cu, Ag, Pt, [49,50,54], and also for the case of Au(100) [55], The origin of this temperature dependence of the island shape is connected with the step line tension, which becomes smaller as the temperature increases, causing the island borders to increase their roughness.The step line tension β of a monoatomic step on a solid surface is defined as the work per length required for generating the step.For uncharged surfaces, this work is equal to the Helmholtz step free energy.Bombis and Ibach reported a step line tension of 170 meV per atom length based on STM images in vacuum.This contrasted with the comparative smaller value of β ≈ 35-70 meV per atom length obtained by Giesen et al. for the unreconstructed Au(100) in an electrolyte solution (sulfuric acid) at higher potentials [56][57][58].Also, previous theoretical calculations have estimated step energies on Au(100) between 65 [59] and 180 meV [60] per step atom.Hence, in order to obtain a direct comparison between the available experimental data and our calculated values, we calculated the step line tension (β A ) for different temperatures ranging from 200-800 K.We calculated the step line tension according to the same procedure used by Bombis and Ibach [55] and described by Schlößer et al. [61] and hence do not describe it in detail here.In this approach, the expectation value G(R) is calculated with (7): where R is the middle radius of a circular island, N is the number of islands used in calculating the average value (N = 40-80), M is the number of angles for which r i (θ, j) is calculated (M = 360), and r i (θ, j) is the radius at a given moment.G is proportional to 1/β A , where β A depends on the direction of the step relative to the surface.In the particular case of gold clusters, G is obtained by (8): From ( 8), we can calculate β A at different temperatures.These results are summarized in Table 2.The corresponding values can also be expressed as a function of the gold atom radius since β A is proportional to the island radius (a = 288 pm).
The relation between the temperature T and β A follows a quadratic dependence.Based on this quadratic dependence and the data reported in Table 2, we obtained a value of 439 meV/nm at 353 K for the step line tension β A .As function of the gold atom radius (β A,a|| ), we obtained a value of 127 meV.At the same temperature (353 K), the value obtained by Bombis and Ibach of 170 meV can be compared to our value [55].The difference between both values might be due to the rate catalog (exchange processes with the surface have so far not been considered) or the difference in the quotient of the step line tensions (β B /β A ) along perpendicular steps B and A. This quotient provides information on how far away the islands are from reaching equilibrium conditions (at equilibrium β B /β A = 1).The value obtained by Bombis and Ibach for this quotient [55] was 0.96 compared with our value of 0.99 [62].
Ongoing work in our group seeks to include the exchange diffusion process and other processes which could influence the kMC description of the Ostwald ripening on Au(100) in future models.

Conclusions
In the present work, we have described a multiscale approach for studying the model case of Au-island ripening on Au(100).Starting with the atomistic level, we performed ab initio DFT calculations on a set of selected self-diffusion processes on perfect and imperfect Au(100) surfaces, obtaining both binding and activation energies.Afterwards, this data was used to generate a reactive ReaxFF force field for Au-Au interactions, which provided the basis for investigating many additional diffusion scenarios.Indeed, using ReaxFF, we were able to consider all possible processes involving nearest neighbors.The thus generated database of diffusion rates was finally used to study the system on the macroscopic level by means of kMC simulations.
Using kMC, we evaluated different initial surface configurations, by which we were able to investigate adatom nucleation as well as island aggregation and fluctuation processes on Au(100).Finally, by averaging over island shapes during the simulation runs at different temperatures, we studied the equilibrium shape of monoatomic high islands at different temperatures (here 500, 750, and 1000 K).In good agreement with experimental results for gold islands on Au(100), we found that increasing the temperature leads to island roughening.We explain this temperature dependence of island roughing in terms of the temperaturedependent behavior of the step line tension.Our step line tensions agree fairly well with the reported experimental values and validate the present multiscale approach, which should be easily adaptable to the study of other systems.

Figure 1 :
Figure1: Schematic of the hierarchy of computational chemistry used in the present work.Starting on the atomic scale, the information obtained via DFT was used to train a reactive force field (ReaxFF), which addresses the molecular scale.ReaxFF was then used to obtain the relevant data for kinetic Monte Carlo simulations to model the mesoscale.While training was performed upwards, each approach was tested against the next higher-level (more accurate) method.

Figure 2 :
Figure 2: Overview of the different diffusion pathways on Au(100) that have been investigated with DFT.

Figure 3 :
Figure 3: Selected Au diffusion pathways on an Au(100) surface (a) bridge diffusion on a terrace, (b) along a step, (c) kink diffusion, and (d) dimer diffusion.Insets in each figure illustrate the diffusion pathways, where filled circles depict adatoms on an Au(100) surface.

Figure 4 :
Figure 4: (a) Schematic representation of a diffusion event.The adatom of interest (blue) diffuses from position 1 to 2, while the next-nearest neighbor positions 3-12 can be either occupied (yellow) or empty.(b) Hard-sphere model of a possible diffusion process.The red atom would correspond to position 1, while the blue atoms are occupied adsorption sites (yellow spheres in (a)) and yellow atoms are the underlying Au(100) surface.

Figure 5 :
Figure 5: Initial input configurations considered in the kMC simulations: (a) random distribution of 3375 adatoms; (b) two islands with 52 and 508 atoms, respectively; (c) round island with 912 Atoms; (d) squared island with 256 Atoms.The Au adatoms are depicted in red while the Au(100) surface atoms are shown in yellow.

Figure 6 :
Figure 6: Snapshots of the nucleation process at 850 K, starting from randomly dispersed adatoms on the Au(100) surface.(a) Initial structure, (b) after 0.25 ns, (c) after 7.1 ns, and (d) after 4039 ns.After 10 9 iteration steps the average of the island size reaches 180 atoms.

Figure 7 :
Figure 7: Snapshots along the kMC simulation at 850 K of the aggregation of a two differently sized islands of 52 and 508 atoms, respectively.(a) Initial structure, (b) after 0.42 ns, (c) after 4.7 ns, and (d) after 95 ns.

Figure 8 :Figure 9 :
Figure 8: Snapshots along the kMC simulation of an initially spherically shaped island of 912 atoms at 850 K. (a) Initial structure, (b) after 0.43 ns, (c) after 4.65 ns, and (d) after 322 ns.

Table 1 :
Activation energies, prefactors, and calculated rate constants (at 300 K) for the different diffusion pathways.

Table 2 :
Calculated step line tension (β A ) and step line tension as function of the gold atom radius (β A,a|| ) at different temperatures.