Modeling of epitaxial film growth of ${\mathrm{C}}_{60}$ revisited

Epitaxial films evolve on time and length scales that are inaccessible to atomistic computer simulation methods like molecular dynamics (MD). To numerically predict properties for such systems, a common strategy is to employ kinetic Monte Carlo simulations, for which one needs to know the transition rates of the involved elementary steps. The main challenge is thus to formulate a consistent model for the set of transition rates and to determine its parameters. Here, we revisit a well-studied model system, the epitaxial film growth of the fullerene ${\mathrm{C}}_{60}$ on an ordered ${\mathrm{C}}_{60}$ substrate (111). We implement a systematic multiscale approach in which we determine transition rates through MD simulations of specifically designed initial configurations. These rates follow Arrhenius's law, from which we extract energy barriers and attempt rates. We discuss the issue of detailed balance for the resulting rates. Finally, we study the morphology of subatomic and multilayer film growth and compare simulation results to experiments. Our model enables further studies on multilayer growth processes of ${\mathrm{C}}_{60}$ on other substrates.

Figure 1

Example MD simulation setup. Right: Top view. Left: Side view. The dark gray particles are fixed in place to ensure an fcc crystal structure with a (111) surface. The three layers of light gray particles are freely evolving during the NVE simulation but are under the effect of a Langevin thermostat or velocity rescaling at the beginning of the simulation. The black particles are set in a stable configuration where they are unlikely to move and set up the environment for the tagged orange particle, which is set in a state where the transitions of interest can be observed. Both black and orange particles are never under the effect of a thermostat or velocity rescaling. In this example setup, the transitions of interest are the descension and diffusion of the orange particle.

Figure 2

Geometry and involved transition types. (a) Visualization of the sublattices. The unit cell is a parallelogram with two sublattice positions. A freely diffusing particle always jumps from one sublattice to the other, and particles that belong to the same cluster always occupy the same sublattice. Periodic boundary conditions are applied to the $x$ and $y$ directions, while the $z$ direction consists of two layers. (b) Visualization of the transition types in the simulated system. (c) Difference between A and B step edge diffusion. While on A step edges the orientation of the base layer facilitates the edge diffusion transition, on B step edges it does not. As a result, dissociation and reassociation can be the more probable trajectory at B step edges.

Figure 3

Example Arrhenius plots. The temperature range is varied depending on the configuration in order to observe as many transitions as possible. Left: Configuration (2NB) in Table 1 with a temperature range of $T\in [600,750]\mathrm{K}$, sampling a total of $\approx \mathrm{120 000}$ trajectories. Right: Configuration (2NA) with a temperature range of $T\in [500,750]\mathrm{K}$, sampling a total of $\approx \mathrm{23 00 000}$ trajectories. Because in configuration (2NA) transition I had a much higher transition rate than transitions II and III, a much larger number of trajectories was needed to sufficiently sample all transitions.

Figure 4

Energy barriers of the different transition types plotted against the number of initial neighbors. The dissociation barriers are taken from (1N)IV, (2NA)II, and (3NA)III, which are the longer dissociation paths compared to the ones in B step directions and are therefore the energy barriers for a complete dissociation. For the energy barriers of ascension there is no significant difference between the A step and B step cases, so we took the average of the two respective measurements for this plot. The data for edge diffusion along an A step edge is composed of (1N)II, (2NA)I, (3NA)I, and (4NA)I, while edge diffusion along a B step edge is given by (1N)I, (2NB)I,(3NB)I, and (4NB)I. All energy barriers are well described by linear fits (lines).

Figure 5

Comparison of the average entropy production between the RawMD and simple models plotted against time $t$. The ensemble average is calculated from 100 trajectories at a temperature of $T=318\mathrm{K}$.

Figure 6

KMC simulation snapshots of the RawMD (left) and simple (right) models at several temperatures (increasing from top to bottom). Parameters for the simple model are $E_B=235\mathrm{meV}$ and ${\nu }_0=0.25\mathrm{THz}$. The snapshots show an area of $1000\times 866{\text{nm}}^2$, and each contains 100 000 molecules, corresponding to $10\%$ coverage.

Figure 7

Average cluster features of the two KMC models in a temperature range from 225 to $360\mathrm{K}$. Left: Comparison of the RawMD model to three simple models with varying effective bond strength $E_B$ and a constant attempt rate of ${\nu }_0=0.25\mathrm{THz}$. Right: Comparison of the RawMD model to three simple models with varying attempt rate ${\nu }_0$ and a constant effective bond strength of $E_B=235\mathrm{meV}$. As suggested by the shape descriptors, the resulting cluster morphologies are very similar.

Figure 8

Cluster density as a function of molecular exposure (solid symbols) at $T=333\mathrm{K}$. Also plotted are results from the optimal simple KMC model: unmodified model for growth on a clean ${\mathrm{C}}_{60}$ substrate (light line) and with modified rates to reproduce the cluster density in the first layer (dark line).

Figure 9

Snapshots at different stages of the multilayer growth simulations at a temperature of $T=333\mathrm{K}$ and a flux of $F=0.1\text{ML/min}$. (a) Simplified KMC model from Ref. [13]. (b) Simple model derived in Sec. 3c. As a result of the two distinct sublattices, grain boundaries form between clusters and persist in the following layers. (c) Using the simple model for the second and higher layers. In the first layer, the diffusion barrier and attempt rate have been altered to reproduce the experimental first-layer cluster densities.

Figure 10

Peak cluster densities of the second ${\mathrm{C}}_{60}$ layer at the three experimentally observed temperatures (disks). Also shown are the KMC prediction of the optimal simple model (diamonds) and the upper bound equation (10).

Figure 11

Top: Density plot of the potential energy of a sample molecule (dark green circle) in a given configuration of fixed cluster molecules (black disks, viewed from above) interacting via the Girifalco potential (the substrate is omitted). The dotted dark green line represents a possible reaction path connecting the four states $S_1,\cdots ,S_4$ passing through the transitional states $T_1,T_2$, and $T_3$. Bottom: Plot of the potential energy along the dotted reaction path. The energy barriers $\mathrm{\Delta }{E}_{ij}$ ($i$ is the number of initial neighbors, $j$ is the number of final neighbors) are only weakly affected by $j$.