A coarse-grained Monte Carlo approach to diffusion processes in metallic nanoparticles

A kinetic Monte Carlo approach on a coarse-grained lattice is developed for the simulation of surface diffusion processes of Ni, Pd and Au structures with diameters in the range of a few nanometers. Intensity information obtained via standard two-dimensional transmission electron microscopy imaging techniques is used to create three-dimensional structure models as input for a cellular automaton. A series of update rules based on reaction kinetics is defined to allow for a stepwise evolution in time with the aim to simulate surface diffusion phenomena such as Rayleigh breakup and surface wetting. The material flow, in our case represented by the hopping of discrete portions of metal on a given grid, is driven by the attempt to minimize the surface energy, which can be achieved by maximizing the number of filled neighbor cells.


Introduction
Metallic particles in the range of a few nanometers straddle the line between microscopic and quasi-classical systems. Many-body effects, which determine the properties of the corresponding bulk phase, are still strongly dependent on the particle size [1]. Characteristic features of the bulk metal are often not yet developed, while other properties such as the chemical behavior can overcome bulk values by orders of magnitude.
Regarding the shape and structure of metallic nanoparticles, modern imaging techniques of electron microscopy yielded significant advances in the experimental study of these transitions with particle size [2], and opened the possibility for the in-situ analysis of structural changes caused by chemical interactions or temperaturemediated diffusion processes [3]. A detailed understanding of the dynamic behavior in these systems requires a multi-scale approach towards material simulation which ideally covers the whole range from the macroscopic to the atomic scale, with the latter often being dominated by less-intuitive quantum effects. Therefore, numerous methods and program packages have been developed, each of them filling a certain niche in the broad spectrum of computational materials research and design. Particularly interesting with respect to multiscale approaches is the MBN Explorer suite of programs [4], which offers the ability to group atoms into larger units in order to reduce the computational effort in molecular dynamics or Monte Carlo simulations.
Contribution to the Topical Issue: "Dynamics of Systems at the Nanoscale", edited by Andrey Solov'yov and Andrei Korol. a e-mail: andreas.w.hauser@gmail.com In this article we present a new approach towards the simulation of diffusion behavior for metallic nanoparticles with diameters in the range of a few nanometers, aiming at finite many-particle systems consisting of several thousands of metal atoms. Special focus will be set on wireshaped nanostructures due to their potential applications in circuits or nanodevices [5,6], as membranes [7], biosensors [8][9][10], transparent electrodes [11], waveguides [12] or antennae in photochemistry [13].
Our approach is based on the concept of a cellular automaton (CA), a well-studied discrete model of computability theory with broad applications ranging from physics to economy [14]. The basic principle, originally conceived by Ulam and von Neumann [15,16], and introduced to a broader audience by Conway in the form of a zero-player game called "Game of Life" [17], is a discretization of the environment in the form of cells. The status of each cell is periodically updated based only on the status of its neighbor cells. This principle has been adapted and implemented here to describe the diffusion of metal atoms on the surface of metallic nanoparticles residing on a carbon substrate. A first, preliminary application of this model can be found in reference [18]. Instead of describing the motion of each atom separately, which introduces clear limitations to the number of particles involved in standard molecular dynamics simulations, we employ a coarse-grained model which handles material transfer in larger groups of atoms, defined by the grid size of the CA. This reduces, together with the intrinsic limitation to interactions with next-neighbor cells, the computational effort significantly, and allows the simulation of temperature-induced effects such as droplet formation, surface tension and Rayleigh-breakup for even nanometer-sized structures. In comparison to the comprehensive MBN package our method can be considered as a "top-down" approach in the sense of accessing nanostructures from the macroscopic side based on microscopy image data, while comparable stochastic Monte Carlo based algorithms, as implemented in the MBN package, are starting at the atomic scale but employ a grouping of atoms into clusters in order to address larger systems. Recently, the latter type of approach has been successfully applied to the description of surface diffusion of Ag 800 nanoparticles on graphite surfaces [19,20], and has also been extended to three dimensions via a stacking of the metallic nanoparticles [21]. As will be shown below, our treatment of surface kinetics via Monte Carlo sampling is very similar to these earlier studies on silver clusters, but with the conceptional difference that in our case fictitious, finite portions of metal, with their volume defined by the resolution of input TEM images, are subject to a random repositioning, and not well-defined concrete metal clusters. In this sense, we are not studying the migration of real nanoparticles along fractal structures, but model the deformation of larger nanowire with lengths of several hundred nanometers.
The basic principle and implementation of our method is presented in the next section. In Section 3 the model is then used to study surface diffusion processes of Ni, Pd and Au nanoparticles and nanowires on amorphous carbon. Our predictions are compared with recent experimental findings of our group obtained by temperatureprogrammed TEM imaging. Section 3.1 is dedicated to the phenomenon of wire breakup and discusses the stability of the selected metallic nanowires at various temperatures. Section 3.2 introduces interactions with the support and aims at a correct description of wire deformation by surface wetting effects. The latter occur after the deposition of inert wires on TEM grids of amorphous carbon. An attempt is made to reconstruct the actual shapes from two-dimensional TEM images and compare them to the model predictions.

Computational details
Surface diffusion processes can be seen as paradigmatic cases of multi-scale dynamics. Although atomic in essence, an atomistic approximation is rarely used to describe this group of phenomena due to their macroscopic character, which only emerges for sufficiently large system sizes and longer time evolutions. Typical examples are surface wetting, the intermixing of liquids, capillarity, Ostwald ripening or Rayleigh breakup. On the large scale, these phenomena are well described by continuous models such as surface tension or Fick's law of diffusion, but a fully macroscopic view has its clear limitations. For example, in the case of solid-liquid phase transitions for metals, the concrete behavior of an arbitrarily shaped, solid object as a function of temperature can only be captured by a spatial discretization in combination with a set of assumptions on local diffusion based on the shape of the material at a given point in time and space.
In some cases, models of reduced dimensionality are able to deliver very accurate results, as has been shown in the case of silver cluster fractals on graphite [19,20]. Experimental studies revealed that silver clusters on graphite arrange themselves in dendritic structures [22,23]. This behavior could be well described by a two-dimensional kinetic Monte Carlo model [4,24] based on the diffusion limited aggregation method [25].
Another interesting attempt to reduce the complexity of the problem has been published in reference [26]. It is based on Mullin's equation [27], originally introduced to simulate the thermal grooving at metallic grain boundaries. It relates the flux of surface atoms to the gradient of the local surface curvature, and derives an artificial force on the respective surface element. This is achieved by a discretization of the contourlines of the actual structure, followed by a propagation in time based on Mullin's equation. Surface diffusion is assumed to be the main form of mass transport and dominating over evaporationcondensation processes. A partial differential equation of the form can be derived, with z denoting the local transversal displacement of a contour line s. A change of the local displacement in time t is assumed proportional to the inverse curvature of the contour line, δ 2 z/δs 2 . All unknown material-specific properties are adsorbed into a constant B, which had been set to one in our previous work for the qualitative description of the degradation of the silver structures [3]. While the discretization in s had been determined by the resolution of the TEM images, the discretization in time was done in arbitrary units. Very fine steps had to be chosen to avoid instabilities in the propagation, which turned out to be only one of several downsides of our original implementation. Nevertheless, this ansatz has been successfully applied to the description of silver nanostructures despite its lack of computational practicality [3]. However, for the efficient description of larger and more complex structures a different route has to be taken, preferably also bearing the potential to capture adsorption effects on different substrates. For this reason an alternative computational concept will be employed: the direct simulation of material transfer via a CA. As will be shown below, this model intrinsically accounts for the locality of the diffusion process via the limitation to interactions between nearest neighbor cells, which keeps the computational effort much lower than techniques based on contour lines. The transfer or "hopping" of small portions of metal will be described by a kinetic Monte Carlo procedure. At first, the actual shape of the material, in our case metallic nanowires deposited on a weakly interacting substrate, needs to be translated into a pattern of cell occupation in a discretized environment (see Fig. 1). The useful cell size is limited to the resolution of the transmission microscopy (TEM) images, in our case to pixels of (0.95 nm) 2 , corresponding to approximately 55-85 metal atoms per cubic cell 1 . The lack of tomographical data in two-dimensional TEM images makes it necessary to introduce additional assumptions on the wire shape due to the loss of information in the projection. Assuming the height of the particles to be identical with the width, i.e. an approximately cylindric shape for wires and a spherical shape for particles, a three-dimensional structure can be derived from the pixel positions and intensities of the TEM image. Following the fundamental principle of CA models, each cell is either thought empty or completely filled with metal atoms. A three-dimensional, discretized model shape derived from a TEM image is shown in Figure 1a. For better illustration, an arbitrary cell and its local environment have been picked and enlarged in Figure 1b. The occupied or filled cell shown in red is surrounded by 5 occupied (golden) cells; the remaining 21 other cells in planes 0, 1 and 2 are empty.
In the second step we introduce a series of update rules to model the time evolution at a given temperature: 1. An empty cell can become filled if there is a filled neighbor cell and the corresponding hopping condition is fulfilled.

A filled cell can become empty if there is an empty
neighbor cell available and the hopping condition is fulfilled. 3. The reaction rate k for a hopping from cell a to b is given by the Eyring equation [28] as described below.
As soon as the rates for all possible hoppings of a filled cell are evaluated, a random number is drawn between 0 and 1. The normalized hopping rates and the random number are then used to determine the actual hopping event.
with k B as Boltzmann constant, h as Planck's constant, T as temperature, R as gas constant and ΔF as a free energy barrier height for the hopping process in kJ/mol. With the definition of the free energy as F = E − T S this equation above can be rewritten as where we have introduced the temperature-dependent entropy. Note that this ansatz is borrowed from reaction chemistry and therefore only valid in the atomic picture, where the relative energy of a saddle point on an energy surface for a given chemical reaction can be interpreted as a reaction barrier height. In the concrete case, where the migration of a whole group of atoms (the atoms contained in a filled cell) is simulated, neither a reaction coordinate nor a barrier height can be clearly defined. However, even such a collective migration can be reduced to a series of atomic transition states in principle, which justifies the analogous mathematical treatment, but leads to the question of suitable expressions for ΔF . Here we introduce a strong assumption known as the Bell-Evans-Polanyi principle in reaction chemistry. It is based on the observation that, for similar chemical reactions, the differences in activation energy are often proportional to the differences in reaction enthalpy [29,30]. In other words, for the same class of chemical reaction -a condition trivially fulfilled for the given problem of surface migration, where the "difference" of a reaction simply refers to different localization of the migrating portion of metal -reaction kinetics can be derived from thermodynamics. We approximate ΔF by a difference of sums over n metal-metal bonding energy terms E mm , with n as the number of bonds formed in the final and broken in the initial position. The corresponding values are taken from the literature [31]. The unknown proportionality constants as well as the temperature dependence of the entropy can be adsorbed into a temperature-dependent pre-exponential factor. Its concrete value can not be derived within this simplified model, but can be determined via a direct comparison of update cycles of the CA to the actual time evolution of an experimentally observed nanostructure at constant time (if available). Details of this approach can be found in a preliminary study with stronger focus on the experimental setup [18]. Note that condition 3 introduces not only the tendency to minimize the energy of movable cells by migration, but also provides the random character which is intrinsic to the temperature-driven process of surface diffusion. In Figure 1c we provide a symbolic representation of the hopping options for the selected cell mentioned above. A cell migration from its current position to any empty cell of plane 2 will be lower in energy due to the higher number of filled neighbors in these future positions. Therefore, a jump into one of these places is most likely to happen.
Double occupations of a given cell must be forbidden if a constant metal density is assumed. This is achieved by tabulating all possible jumps for a given time step, followed by a random permutation of the jumping events before the actual migration takes place, in order to eliminate any bias introduced by the serial processing of the cells.

Simulating surface diffusion in monometallic nanoparticles
Recently, our group focused on the experimental study of surface diffusion processes in nanowires formed of coinage metals and nearby elements [3,18], using superfluid He nanodroplets for the production of one-dimensional metallic structures [32,33]. The dotation of helium droplets, a technique [34][35][36] originally established for spectroscopy investigations, has recently been extended by new facet, the controlled synthesis and structure preserving soft deposition of metal clusters [37][38][39][40][41]. A peculiarity of rotating superfluid He droplets with diameters of a few μm, the formation of quantized vortices [42,43] within the liquid [44][45][46], turns out to emphasize a preferably one-dimensional growth of metal structures [47][48][49] in the liquid helium environment due to a pressure gradient effect around the vortex core [50,51]. The advantage of this unusual and experimentally challenging technique is the ability to create metallic structures in a fully inert environment, without the additional stabilization through templates [52][53][54], which introduce a large bias in stability measurements.
These recent experiments are the basis for the computational studies performed in the following sections.

Rayleigh breakup of metallic nanowires
We start with the simulation of the breakup behavior observed via TEM imaging of Au nanowires on a heated substrate. At first, the TEM image of freshly deposited nanowires at room temperature is discretized on the grid. The obtained structure is shown in the upper left image of Figure 2. Note that the actual grid of the CA is threedimensional, but only a 2D projection is shown for easier comparison to the TEM images. In the very first step a slight smoothing of the structure appears as the naturally occurring noise of the image translation is removed instantly. In the following series of cell update steps or "cycles" in CA terminology the shape shifting effects of surface diffusion become evident. After 20 cycles a thinning of certain regions becomes apparent, which gets more and more emphasized (cycle 85) until the first breakup occurs (cycle 100), here at the lower end of the wire structure. A series of additional breakups follows (cycles 200-500). After that, the qualitative changes in the structure between cycles become marginal. For very long propagation in time (cycle 3000) all remaining fragments reach their local minimum of spherical shape, i.e. a circle in the twodimensional projection.
As can be seen in Figure 2, the simulation agrees well with the experimental outcome and with the dynamics observed in the TEM study. Breakup sequences and positions are fully reproduced by the simple CA model. However, a direct comparison of different metals with respect to stability is not straightforward as material properties only enter via hopping probabilities. Therefore, the Rayleigh breakup patterns are material-independent in the sense that infinite time propagation will lead to the same results for all metals. However, the process itself will be faster in cases where metallic binding energies are higher since they enter the hopping probability in the exponent. This puts the cycles needed for a breakup in a given structure into a direct relation with the material and allows to rank different metals with respect to their stability at a given constant temperature. Details of this investigation are published in a separate article [18].
The accelerating effect of higher temperatures is also captured by the ansatz given in equation (2). However, in order to be able to predict breakup events in time, it becomes necessary to estimate the amount of material displaced in each cycle. In principle, this can be achieved by a direct comparison of the number of cycles needed for a breakup to the real breakup time observed with the TEM at a given temperature.

Surface wetting on the nanoscale
In this section we extend the simulations by the introduction of additional parameters which characterize the support material in more detail. We analyze the wetting behavior of a selection of metals (Ni, Pd, Au) on amorphous carbon, with the aim to reproduce the typical shapes obtained for originally spherical nanoparticles deposited on a substrate. Fig. 2. CA simulation of the surface diffusion for an Au nanowire structure derived from TEM images. For a given temperature, the propagation in time is simulated by a series of update cycles for the CA. A reshaping of the nanowire becomes evident, which leads to multiple Rayleigh breakup events, labelled with Greek letters in chronological order. Eventually, qualitative changes in the structure between cycles become marginal as the remaining fragments reach their local minimum of spherical shape. Two HAADF-TEM images have been added for direct comparison. At first, we characterize the wetting situation for a given metallic nanoparticle based on simple geometric arguments as illustrated in Figure 3. Assuming an incomplete sphere as the shape of the metallic nanoparticle after adsorption onto the substrate, we can define an angle α which determines the size of the missing spherical cap. In this simplified model, an analytical relation between α and the intensity profile I(x) for such an adsorbed and deformed particle can be derived, (4) with the additional relations h = r(1 − cos (α/2)) and s = 2r sin (α/2), where r denotes the particle radius, h the height of the cap, l(x) the thickness of the unmodified upper half of the sphere, s the radius of the cap, and H(x) the Heaviside step function. Equation (4) is based on the assumption of a linear relationship between the thickness of the particle (with thickness defined as particle width in the direction of observation) and the intensity in the twodimensionsal scanning-TEM image. This well-known feature holds as long as no saturation effects occur, i.e. for metallic particles with diameters of only a few nanometers, and is commonly used to obtain additional morphology and composition information in high-angle annular dark-field (HAADF) scanning transmission electron microscopy [55][56][57].
A fit of this simplified model to the real TEM data for small Ni, Pd and Au nanoparticles is shown in Figure 4, together with the corresponding images. We find a very small α value for Au (9 • ) when compared to Pd (31 • ) and Ni (53 • ), which indicates a significantly reduced wetting in the case of gold. This is also deducible from the slightly more abrupt change of TEM contrast at the edges of the Pd and Ni particles in comparison with Au.
In a second step, we extend the CA approach described in Section 2. The interaction with the substrate is included by means of added binding energies for the occupation of cells directly above the surface plane. For a particle of a given size and metal the outcome is only determined by the ratio of the metal-metal and metal-surface interaction energies. A series of simulations with various ratios is depicted in Figure 5, showing the converged result of the diffusion process for cases from minimal to significant surface wetting. From this, an approximately linear relation can be derived between the angle α (53 • , 31 • and 9 • for Ni, Pd and Au) which defines the size of the spherical cap in the geometric model and the actual ratio of metal-metal energies E mm to metal-surface energies E ms :  (4) to the experimental data with α as free parameter. This relation allows for a crude estimation of E ms values of Ni, Pd and Au with the amorphous carbon surface based on the corresponding metal-metal interaction energies.
Despite its simplicity, this method produces relative surface interaction energies which are in reasonable agreement with ab initio-derived adsorption energies for the selected (111) metal surfaces on graphene [58]. Table 1 contains the metal-metal interaction energies E mm which we used as input parameters for the CA, and compares the derived metal-carbon interaction energies E ms to surface adsorption energies from reference [58]. The approximate number of metal atoms per cell n v as well as the average number of bonds per filled cell neighbor n b can be derived from the ionic radii r i listed in the second column. From this, we can estimate the energy gain which occurs due to the interaction of the filled cell with its 9 "substrate" neighbor cells. The latter can be directly compared to the ab initio values from reference [58] for the interaction Table 1. Comparison of Emm, the metal-metal interaction energies [31] (the only input data of the CA) to derived Ems values and the corresponding ab initio surface adsorption energies ΔEC taken from reference [58]. Absolute energies are given in kJ/mol, ionic radii ri in nm. with graphene. An atomic density of 37 carbon atoms per (0.95 nm) 2 has been assumed for this comparison, derived from the well known molecular structure of graphene [59] with a C-C bond length of 1.42Å.

Emm
The largest discrepancy between our prediction and the ab initio value occurs for nickel, which might be related to a reduced image contrast for the lighter Ni atoms and due to more pronounced oxidation effects. Furthermore, nickel is known as a wetting electrode on graphene at higher temperatures, forming carbides at the contact surface [60][61][62].

Conclusion
A kinetic Monte Carlo method based on the cellular automaton principle has been introduced which allows the simulation of surface diffusion processes in metallic nanoparticles. For the simulations, realistic shapes are derived from TEM images and translated into discrete 3-D models. The grid size for the model structures is determined by the resolution of the TEM image and corresponds in our case to cells of (0.95 nm) 3 or approximately 55-85 metal atoms depending on the element. A series of cell update rules is then defined in order to simulate surface diffusion by cell hopping events. Hopping probabilities are estimated based on the difference between the number of occupied neighbor cells before and after the hopping. To introduce an element-specific hopping probability we derive a set of parameters from known binding energies between the elements involved. With this ansatz, even the behavior of larger metallic structures can be correctly predicted and reproduced. The reduction to nearest-neighbor interactions, a fundamental aspect of the cellular automaton model, allows for a highly efficient computation of surface diffusion processes.
We showed that this computational approach is able to predict Rayleigh breakup phenomena observed for metallic nanowires with diameters of a few nanometers. It reproduced the final shape after heating treatments and correctly predicts the chronological order and the location of wire breakup events for any given shape, as could be verified by in situ TEM-imaging of heated samples. Assuming a non-wetting behavior for all studied metals on amorphous carbon, the cellular automaton produces material distributions which are in good agreement with the particle shapes derived from the TEM images. We further extended the model by explicitly introducing interactions with the substrate, which allowed us to reproduce the deformations we observe for the Ni, Pd and Au nanostructures due to surface adsorption effects. The relative strength of metal-carbon interactions derived from our TEM-inspired CA approach agrees reasonably well with the available ab initio data for the binding energies of graphene on bulk (111) surfaces of the selected metals.