Microscopic Origin of the Apparent Activation ... - ACS Publications

Article pubs.acs.org/JPCC

Microscopic Origin of the Apparent Activation Energy in DiffusionMediated Monolayer Growth of Two-Dimensional Materials Miguel A. Gosálvez*,†,‡,§ and Joseba Alberdi-Rodriguez§ †

Department of Materials Physics, University of the Basque Country UPV/EHU, 20018 Donostia-San Sebastian, Spain Centro de Física de Materiales CFM-Materials Physics Center MPC, centro mixto CSIC−UPV/EHU, 20018 Donostia-San Sebastian, Spain § Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastian, Spain ‡

S Supporting Information *

ABSTRACT: Current trends indicate that mass production of highquality, two-dimensional materials will likely be based on the use of onsurface synthesis (OSS) and/or chemical vapor deposition (CVD) on selected substrates. However, the success of these techniques heavily relies on a deeper understanding of terrace and perimeter diffusion of the adspecies across and around many self-generated, mobile clusters, and/or motionless islands with a variety of shapes (dendritic, compact, irregular, polygonal, ...). We show that, for typical monolayer growth conditions at constant deposition rate, the total square distance traveled by all adsorbed particles departs from the total number of diffusion hops due to the onset of correlations between subsequent hops along the perimeters of a growing density of obstacles. As a result, we propose a new expression to determine the tracer diffusivity of the adparticles, directly providing a simple understanding of the temperature and coverage dependence of this observable. Most importantly, we describe the microscopic origin of the apparent diffusion barrier extracted from a typical Arrhenius plot at any given coverage.

■

INTRODUCTION Graphene and related compounds are ideal candidates to achieve superior energy efficiency by enhancing the storage and release of charge and fuel.1,2 However, this requires the production of high-quality, large amounts of these twodimensional (2D) materials. On-surface synthesis (OSS) by evaporation of molecular precursors at moderate temperatures has recently received much attention as it enables atomic precision molecular assembly, generating defectless nanoribbons and complex 2D molecular networks on selected substrates.3−7 Another alternative is chemical vapor deposition (CVD), which is considered suitable for batch production of large amounts of extended monolayer samples with sufficient quality at reasonable cost.8−14 A common feature for OSS and CVD is the prominent role of diffusion. In a typical CVD process, a metal substrate is heated up to about 1000 °C in the presence of a carbon source vapor. As a result, the vapor is thermally decomposed and adparticles are adsorbed on the metal surface, where they diffuse randomly, eventually forming small clusters at random locations (nucleation). The clusters gradually evolve into larger islands by accretion of the diffusing adparticles (growth), and the islands eventually merge to form a single layer (coalescence). The quality of the resulting monolayer depends markedly on the density and structure of the grain boundaries formed during coalescence, which depends on the actual shape of the islands (dendritic, compact, irregular, polygonal, ...), © 2017 American Chemical Society

which is ultimately a function of the relative occurrence and strength of perimeter diffusion.15−17 At lower temperatures and with larger diffusing particles, the same picture is valid for OSS. In fact, both dendritic and compact islands can be obtained in OSS, depending on the details of perimeter diffusion.5 Thus, one cannot overstate the importance of diffusion for the future production of 2D materials. Previous work on surface diffusion has focused on such nucleation, growth and coalescence aspects, analyzing the time evolution of the monomer and island densities, characterizing the island size distributions, and describing the origin of the island capture numbers and capture zones.18−22 Similarly, much effort has been invested on clarifying the differences and similarities between the tracer and collective diffusion coefficients, and describing them for various systems.23−29 The tracer diffusivity, DT, which is the relevant transport coefficient when the particles can be followed, can be determined using three alternative formulations, e.g., eqs 1, 3, and 19 in ref 28. In this study, we present a fourth expression, describing DT in terms of the total hop rate, Rh, defined as the sum of all distinct hop rates accessible to the adparticles, each multiplied by the corresponding hop multiplicity (the number of hops that share a given rate). The proposed formula Received: June 13, 2017 Revised: August 4, 2017 Published: August 21, 2017 20315

DOI: 10.1021/acs.jpcc.7b05794 J. Phys. Chem. C 2017, 121, 20315−20322

Article

The Journal of Physical Chemistry C

to jump to a neighbor site of type j; and Mij = nimij is the corresponding total multiplicity, i.e., the number of hops with rate νij that can be performed. Here, ni is the number of particles of type i and mij is the local multiplicity, i.e., the number of sites of type j surrounding a particle of type i. For triangular and square lattices the maximum coordination number S is 6 and 4 (7 and 5 particle types), respectively. In general, any suitable variable−or set of variables−may be used to uniquely describe the different particle types. Based on transition state theory (TST),31 the hop rates in eq 2 are written as νij = νae−Eij/kBT, where the Boltzmann factor e−Eij/kBT indicates the probability to perform the hop at temperature T if the energy barrier is Eij, and the attemptfrequency νa (≈ 1 × 1013 Hz) describes the frequency with which the substrate phonons and the adparticle vibrations couple with each other.27 Regarding eq 3, Fi is the adsorption rate at a site of type i and nϕi is the number of empty sites of type i. A typical assumption (adopted in the right-hand-side of eq 3) is to take all adsorption N rates equal to a single deposition flux (Fi = F). θ = L La

accounts for small-cluster diffusion across the terraces (monomers, dimers, trimers, ...) as well as perimeter diffusion at mobile clusters and motionless islands of any size. In particular, our approach provides an exact formulation of an empirical procedure,23 where DT is described as the product of an average hop rate, Γ, and a correlation factor, f T, which accounts for all memory effects between consecutive hops at finite coverages.23−27 Such correlations refer to the fact that an adparticle hopping from site i to site j leaves site i empty and, thus, at finite coverage the adparticle has a higher chance of returning to i. We show that Γ must be understood as a triple average of the total hop rate, Γ = ⟨R h⟩/Na, where ⟨X⟩ is the ensemble average of X, X̅ =

∫ X dt Σ X Δt = Σk Δk t k is the time average ∫ dt k k

of X, and X/Na is the “per-particle” average of X, with Na the number of adparticles. Similarly, we stress that f T should be understood as the proportionality factor between the total square distance traveled by all adparticles, ⟨R2⟩, and the total number of diffusion hops, ⟨Nh⟩, namely: fT =

⟨R2⟩ l ⟨Nh⟩

x y

2

designates the total coverage, where LxLy is the total number of adsorption sites (before adsorption of any particle). Since the increase in coverage with time satisfies the relation dθ = F(1 − θ ), direct integration gives θ = 1 − e−Ft. Thus, dt initially θ increases linearly (θ ≈ Ft) and later saturates at value n 1. Defining θi = L Li as the density of particles of type i, the

(1)

with l the hop distance between adjacent adsorption sites. At very low coverage the particles perform completely independent random walks and ⟨R2⟩ = l2⟨Nh⟩.24,27 Thus, f T = 1. At finite coverages, however, each particle is affected by the presence of (and the interaction with) the other particles, and thus, their random walks become correlated.26,27 As a result, f T departs N from 1. Here, R2 = Σi=1a |ri(t) − ri(0)|2, with ri(t) denoting the position of particle i at time t. Although some studies have been devoted to understanding the behavior of the prefactor and the apparent diffusion barrier extracted from an Arrhenius plot of the temperature dependence of the diffusivity,23,24,26,30 the actual microscopic origin of these quantities is still poorly understood. As shown below, this issue becomes clear once DT is described in terms of the distinct hop rates and their multiplicities, as proposed here.

x y

total coverage satisfies the relation: θ = Σiθi. Similarly, Na = Σini. Regarding eq 4, Gi is the desorption rate for a particle of type i. For completeness, the desorption term is shown in our equations, although it can be neglected in both OSS and CVD. Note that we consider perimeter detachment events ending at terrace sites. Now, we propose describing the tracer diffusivity as follows: DT = g

■

THEORY Let us consider a two-dimensional lattice-gas in order to model a collection of adparticles on a substrate. The substrate is treated as a two-dimensional lattice of adsorption sites, where particles from the surrounding environment are deposited randomly (on empty sites) while the already existing adparticles are able to either hop diffusively to a neighboring (empty) site or desorb (back to the environment). Since adsorption, desorption, and diffusion events may occur at any given time, let us consider the corresponding total event rate Re = Rh + Ra + Rd, with Rh, Ra, and Rd the total hop rate, total adsorption rate, and total desorption rate, respectively: Rh = Ra =

∑ Fni iϕ = F(1 − θ)LxLy i

Rd =

i,j

∑ Gini i

(5)

⟨Nh⟩ t

(6)

= gl 2fT ⟨R h⟩

(7)

= gl 2fT

= gl 2fT

∑ ⟨Mij⟩νij (8)

i,j

= gl 2fT

∑ ⟨ni⟩⟨mij⟩νij (9)

i,j

Here, t is time and g =

∑ Mijνij = ∑ nimijνij i,j

⟨R2⟩ t

1 , 2dNa

where d is the dimensionality

(=2). eq 5 is based on previous definitions.27,28,32 eq 6 is equal to eq 5 based on the fact that ⟨R2⟩ = f Tl2⟨Nh⟩ (eq 1). In turn, the equality between eqs 6 and 7 is based on the observation that dNh/dt is nothing else but the number of hops per unit time, i.e., the total hop rate R h . Correspondingly,

(2)

(3)

(4)

Rh =

In eq 2, i, j = 0, 1, 2, ..., S denote the particle/site type, taken as the coordination number, i.e., the number of in-plane occupied nearest neighbor sites; νij is the hop rate for a particle of type i

ΣkR h , k Δtk ΣkΔtk

=

Σk

ΔNh , k Δtk Δt k

t

=

Nh t

is an exact expression. Finally,

the identity between eqs 7, 8 and 9 is based on the exact decomposition proposed in eq 2. Note that, defining the correlation-f ree dif f usivity as 20316


Article


Figure 1. A log−log plot of the tracer diffusivity (DT = g⟨R2⟩/t), correlation-free diffusivity (D̂ T = gl2⟨Nh⟩/t), correlation factor ( f T = DT/D̂ T), total coverage (θ), density of islands (Nisl), and density of particles of type i (θi) as a function of θ for the triangular lattice at (a) 150, (b) 250, and (c) 750 K. In all frames θ4+ = θ4 + θ5 + θ6, F = 5 × 104 ML/s, and LxLy = 800 × 924 sites. Insets: Typical morphologies at θ = 0.30 (subregion of 200 × 231 sites). Error bars for g⟨R2⟩/t and gl2⟨Nh⟩/t are smaller than the marker size for θ ≳ 10−3 and about double the marker size at the constant-diffusivity regions 10−5 ≲ θ ≲ 10−3.

D̂T = gl 2

Ea = E f +

⟨Nh⟩ t

(10)

α

then eqs 5 and 6 allow one to express the correlation factor (eq 1) equivalently as D fT = T D̂T

ωα =

■

(11)

d(fT ∑i , j ⟨Mij⟩νij) 1 dβ fT ∑i , j ⟨Mij⟩νij

where l and g =

1 2dNa

(12)

in eq 8 are considered independent of T

(and β). [Na = LxLyθ = LxLy(1 − e−Ft) only depends on the deposition flux.] Then, applying the chain rule to eq 12 and M

f

writing νij ∝ e−Eijβ, ⟨M ij⟩ ∝ e−Eij β , and f T ∝ e−E β easily leads to Ea = E f +

∑ ωij(Eij + EijM ) i,j

(13)

where the weights ωij describe the relative contribution of the hops with rate νij to the total hop rate, i.e., the probability of observing a hop with rate νij: ωij =

⟨Mij⟩νij Σij⟨Mij⟩νij

⟨Mα⟩να Σα⟨Mα⟩να

(15)

(16)

COMPUTATIONAL METHOD To illustrate our theoretical formulation we present computational simulations of monolayer growth using the Kinetic Monte Carlo (KMC) method.16,17,21,27,33−39 We use a typical lattice-gas model and the rejection-free, time-dependent implementation.17,27,34,35 Although the algorithmic details are provided in the Supporting Information, here we stress the fact that the time increment is calculated as Δt = −log(e)/Re, where e ∈ (0, 1] is a uniform random number and Re is the total event rate (see eqs 2−4). With some exceptions,37 self-learning KMC studies usually indicate that only a few key barriers are important in order to describe surface diffusion, even if the KMC simulations automatically generate many dozens of barriers.38,39 Thus, in this study the activation energies are chosen in order to generate realistic morphologies for both triangular and square lattices while keeping the resulting models (i) simple enough (to enable a foundational analysis of the apparent activation energy of the tracer diffusivity) and (ii) disengaged from any specific material (to maintain a general discussion for both OSS and CVD). Although the particular activation energies are listed in Table 1S of the Supporting Information, we indicate here that (i) they fulfill detailed balance27,33 and (ii) a few barriers for the triangular lattice are so large that they can be regarded as ∞ for the discussions below. This way there are only four (six) different hop rates for the triangular (square) lattice, referred to as να = νae−Eα/kBT, with α = 0, 1, 2, 3 (α = 0, a, 1, b, c, 2). To keep the discussion general, we use a wide range of temperatures and deposition fluxes (triangular lattice, T = 50 to 1050 K and F = 5 × 102 to 5 × 106 ML/s; square lattice, T = 120 to 600 K and F = 3.5 × 10−4 to 3.5 × 104), and the typical hop distance is taken as l = 1. This choice of l scales all presented diffusivity results, leaving the overall analysis unaffected.

To our best knowledge this is the first study where the tracer diffusivity (eq 5) is exactly linked to the total hop rate (eqs 7−9). Mathematically, the apparent activation energy (Ea) of DT is defined as Ea = −d(log DT)/dβ = −(1/DT) dDT/dβ, where β = 1/kBT. Thus, using eq 8 gives Ea = −

∑ ϵα , with ϵα = ωα(Eα + EαM )

(14)

Note that ωij depends not only on the value of the hop rate itself, νij, but also on the number of hops with that rate at a given time, i.e., the hop multiplicity Mij. In this manner, eq 13 includes the dependence on both coverage and morphology (through their effect on the actual values of the multiplicities), in addition to the temperature dependence (through the activation energies assigned to the correlation factor, Ef, the multiplicities, EM ij , and the hop rates, Eij). For the triangular (square) lattice used here (see next section), the notation can be simplified in terms of the four (six) distinct hop rates, να, and the corresponding hop multiplicities, Mα:

■

RESULTS For typical monolayer growth at constant deposition flux, the square distance traveled by the particles will grow linearly with 20317


Article

The Journal of Physical Chemistry C time only at the very beginning, before the density of deposited particles becomes too large and collisions start affecting the linearity through the formation of dimers, trimers and larger clusters/islands. In addition to the gradual reduction in the free area for terrace diffusion, the perimeter hops occurring in dimers, trimers and larger islands will result in limited contributions to the overall traveled distance per unit time, because the perimeter hop rates are usually slower than the terrace hop rates. Thus, the slope of ⟨R2⟩ determined through DT should decrease with time (and coverage), as shown in Figure 1 for the triangular lattice (see Figure S1 in Supporting Information for the square lattice). Together with the density of monomers (θ0) the density of islands containing two or more atoms (Nisl) is shown in Figure 1 to emphasize (i) the similarity of our results to previous studies,18−22 and (ii) the correlation between increasing obstacle density and decreasing diffusivity. At low coverage (θ ≈ 0) an early stage of constant diffusivity is observed. Here, the density of islands (Nisl) remains significantly smaller than the density of monomers (θ0). Similarly, the densities of the particles with one or more neighbors (θi, with i > 0) remain negligible with respect to θ0. Thus, the value of the tracer diffusivity essentially corresponds to the monomer contribution alone: DTθ≈ 0 =

l2 ⟨R h⟩ 2dNa

Figure 2. Universal Arrhenius plot of the tracer diffusivity (DT = g⟨R2⟩/t) and correlation-free diffusivity (D̂ T = gl2⟨Nh⟩/t) vs inverse temperature for the triangular lattice at θ = 0.30 (F = 5 times 102 to 5 times 106 ML/s, 400 × 462 sites). In addition, the correlation factor ( f T = DT/D̂ T) is plotted on the right-hand-side axis. Values for the apparent activation energy of DT and f T are given (in eV) for the middle flux and the three identified regions (I, II, and III). Inset: Same data for the square lattice (F = 3.5 × 10−4 to 3.5 × 104 ML/s, for 400 × 400 sites).

3

= 4 ν0l 2 .

Here, the total number of deposited atoms is taken as Na ≈ n0 and the total hop rate as (eq 2): Rh = Σjn0m0jν0j = n0Σjm0jν0j = n06ν0. In addition, n0 ≈ Na = LxLyθ ≈ LxLyFt, which increases l i n e a r l y w i t h t i m e . T h u s , ⟨n0⟩ = n 0 / 2 a n d ⟨R h⟩ = ⟨n0⟩6ν0 = n03ν0 . For the square lattice, Rh = n04ν0, 1 ⟨R h⟩= n02ν0 and DTθ≈ 0 = 2 ν0l 2 . After the constant stage the tracer diffusivity decays rapidly with coverage. Typically, DT = g⟨R2⟩/t departs from D̂ T = gl2⟨Nh⟩/t with increasing θ, following the underlying separation of ⟨R2⟩ from l2⟨Nh⟩. This is best shown by plotting the correlation factor f T = ⟨R2⟩/l2⟨Nh⟩ = DT/D̂ T as a function of θ on a log−linear diagram for a wide range of temperatures (see Figure S2 in the Supporting Information). At vanishing coverage, the particles perform completely independent random walks, and f T fluctuates around the theoretically correct value of 1. However, the fluctuations in f T give way to a systematic variation with θ, as many hops take place along the contours of the clusters and islands, and thus each particle’s random walk is not absolutely random any more. This is traditionally referred to as the onset of correlations or memory effects.26,27 Near full coverage, for instance, the hops are extremely correlated since they take place back and forth between pairs of neighboring sites. Nevertheless, f T usually takes values between 0.5 and 2, and thus, it remains essentially constant when compared to the strong decay of 2−3 decades exhibited by DT (Figure 1). This highlights the fact that the coverage dependence of DT is essentially contained in the coverage dependence of ⟨R h⟩ (eq 7). The mainframe of Figure 2 demonstrates that DT = g⟨R2⟩/t and D̂ T = gl2⟨Nh⟩/t are essentially identical over at least 6 orders of magnitude and a wide range of temperature and flux (three and five decades, respectively) for the triangular lattice. The inset shows similar results for the square lattice. As before, the difference between the two variables is shown by the correlation factor, f T, which now exhibits a systematic dependence on T, while the variations are small (≈ 0.5−2) compared to the many decades of variation exhibited by DT. Note that, systematically, f T can be both larger and smaller than

1. In these Arrhenius plots, both the abscissa and ordinate have been multiplied by powers of F, collapsing all DT data into a single curve (above a threshold temperature). Although only one coverage value is shown (θ = 0.30), similar curves are obtained for other values (see Figure S3 in Supporting Information). Accordingly, the diffusivity DT = g⟨R2⟩/t in Figure 2 can be split into three regions (I, II, and III), corresponding to low, medium, and high temperature. Within each region, DT exhibits a seemingly linear behavior, with an increasing slope (or apparent activation energy, Ea) as we go from the low- to the high-temperature region. In region I, DT converges to the straight behavior as F decreases. Since T is very low here, only monomer diffusion is active and diffusion stops as soon as the monomers collide with other monomers or islands. Thus, the diffusivity is literally collision-quenched at high flux while the hops are increasingly allowed for longer distances as the flux is reduced, converging to the straight line. In region II, in turn, the temperature is high enough to activate the lowest-energy perimeter hops, typically occurring between one-coordinated sites. Thus, the diffusivity should contain a contribution from type-1 particles, in addition to the hops by the type-0 particles (monomers). Similarly, in region III the next lowest-energy hops have been activated, typically between two-coordinated sites, and now the diffusivity should also contain a contribution from type-2 particles. Because of the universal behavior (Figure 2) we focus on analyzing the activation energy for one particular flux (F = 5 × 104 and 3.5 × 100 ML/s for the triangular and square lattices, respectively). First, we examine one coverage value (θ = 0.30), considering all coverage values afterwords. Figure 2 shows typical values (in eV) for the apparent activation energy (Ea) of DT = g⟨R2⟩/t in regions I, II, and III, directly obtained by the shown linear fits. Similarly, Figure 2 also shows values for the activation energies (Ef) assigned by similar linear fits of the correlation factor f T = DT/D̂ T = ⟨R2⟩/l2⟨Nh⟩ in the three regions. Additionally, Figure 3a shows a typical Arrhenius plot of the multiplicities ⟨Mα⟩ and the corresponding values for their activation energies (EM α ) obtained by analogous linear fits. Note 20318


Article


Figure 4. Apparent activation energy (Ea) of the tracer diffusivity (DT = g⟨R2⟩/t) vs coverage (θ) for the triangular lattice in the regions I, II, and III defined in Figure 2. Ea is described well by Ef + Σαϵα, where ϵα f = ωα(Eα + EM α ). The relative error |1 − (E + Σαϵα)/Ea| is plotted using 4 the right-hand-side axis. (F = 5 × 10 , 400 × 462 sites). Inset: Same information for the square lattice. (F = 3.5 × 100, 400 × 400 sites).

Figure 3. Arrhenius plot of (a) the time-averaged hop multiplicities ⟨Mα⟩, and (b) the relative contributions ωα to the total hop rate (triangular lattice, θ = 0.30, F = 5 × 104, 400 × 462 sites). Inset in part b: ωα in linear scale, stressing the reduction observed in ω0 at high T.

also be included if the fine details of Ea are desired with utmost accuracy. The relative error |1 − (Ef + Σαϵα)/Ea| in Figure 4 reaches a maximum of 5.59% and 4.93% for the triangular and square lattices, respectively. The major source of error is the constant approximation used for the weights ωα within each region, as shown by the horizontal fits in Figure 3b. This is not surprising, since the three regions with linear behavior highlighted in Figure 2 are an oversimplification, used here to introduce the main ideas in simple terms. Truly speaking, the behavior of the diffusivity is not purely linear in each region and, thus, the slope (or apparent activation energy) changes continuously from one temperature to the next. Thus, instead of three regions, one should consider as many regions as there are measured temperatures and make use of every value of ωα and the actual (local) slope of ⟨Mα⟩ at every temperature. As an example, using nine regions (instead of three) reduces the maximum error to about 2.82% and 3.48%, respectively (see Figure S4 in Supporting Information). Although a similar analysis for each temperature will reduce the relative error to a negligible level, we do not feel the need to present it, since the picture is already clear. From a general perspective, the overall similarity between the results given in Figure 4 for the triangular and square lattices indicates that the proposed theoretical model for the interpretation of the apparent activation energy is rather robust. Although there are differences in the detailed behavior of the two lattices in Figure 4, every aspect of their behavior can be described well by the proposed theoretical framework.

that both f T and ⟨Mα⟩ may decrease with temperature, thus resulting in negative activation energies. Finally, Figure 3b shows the relative contribution ωα to the total hop rate. Figure 3 is for the triangular lattice, obtaining similar plots for the square lattice (not shown). Regarding region I of Figure 3b, the weights satisfy ω0 ≈ 1 and ωα ≈ 0 for α > 0. Thus, only the monomer hops have a significant contribution to the diffusivity, as expected. In region II, however, ω1 ≈ 1−3% and ω2 ≈ 0−5%. Thus, the hops from single-coordinated sites (α = 1) and doubly coordinated sites (α = 2) need to be taken into account if the value of the diffusivity must be accounted for accurately. This is unavoidable in region III, where ω1 ≈ 4−25% and ω2 ≈ 5−2%. By performing similar linear fits for all coverage values between 1 and 100% we are able to explain any particular value of the apparent activation energy (Ea) of DT = g⟨R2⟩/t in terms of the contributions Ef and ϵα = ωα(Eα + EM α ) defined in eq 15. This is shown in Figure 4 for both the triangular and square lattices. In region I, the monomer contribution alone (α = 0) accurately describes the strong reduction in Ea at small coverage and the constant behavior at higher coverage. Since Ef ≈ 0 and ω0 ≈ 1 in this region, eq 15 gives Ea ≈ E0 + EM 0 , and the observed behavior is due to the change in the slope of ⟨M 0⟩ with coverage, i.e., the θ-dependence of EM 0 . In regions II and III, explaining Ea requires considering more contributions for both lattices. For instance, in region III, four contributions are needed (α = 0, 1, 2, 3 for the triangular and α = 0, a, 2, b for the square lattice). Note that here Ef is negative and, thus, Ea = Ef + Σαϵα is smaller than Σαϵα. In addition, the contributions ϵα of the perimeter hops with increasingly higher energy (i.e., between sites with increasingly higher coordination numbers) become more relevant as full coverage is approached. In this manner, it becomes clear that the apparent activation energy cannot be explained by simply averaging the corresponding atomistic hop barriers, (ΣαEα)/S, nor by considering a weighted average of only the hop barriers, ΣαωαEα. The contribution from the hop multiplicities, ΣαωαEM α , must be included, because they depend markedly on both temperature and coverage. In addition, the contribution from the correlation factor, Ef, must

■

DISCUSSION We stress that the present study has a strong theoretical foundation, with eq 13 as the primary result. On the basis of knowledge about the relative activity of the processes (ωij) and the temperature dependence of the multiplicities (EM ij ) and correlation factor (Ef), eq 13 provides a theoretical link between the apparent activation energy of the tracer diffusivity (Ea) and the microscopic activation energies of the elementary processes (Eij). Within a purely conceptual lattice-gas model (in which the KMC simulations need not be included), eqs 1 through 11 provide exact mathematical equalities. Note that the 20319


Article

The Journal of Physical Chemistry C neighborhood-dependent hop rates (νij) appearing in eq 2 do not depend on time (nor coverage, other than from the neighborhood dependence itself), and thus, they appear without modifications in eqs 8 and 9. On the other hand, eq 12 makes the assumption that the lattice parameter of the substrate essentially remains constant with temperature. This ensures that both the hop distance (l) and the dimensions of the substrate (Lx and Ly) also remain constant. This is justified considering that typical variations of the lattice parameter of most materials are limited to only a few percent in a wide range of temperatures, in comparison to the corresponding large variations by orders of magnitude observed in the tracer diffusivity itself. Finally, beyond the standard TST formulation f of νij, we allow Eij, EM ij , and E to vary with temperature, thus

where ⟨Wi ⟩ =

M

f

e−E β, respectively, for the derivation of eq 13 from eq 12. The variation of Eij with temperature is not used here, but it might be useful in future studies. As a result, the only assumptions made in order to derive eq 13 are that the lattice parameter of f

the substrate and the prefactors of e−Eijβ, e−Eij β and e−E β remain essentially constant with temperature. Should any of these variables turn out to be strongly dependent on temperature, then their actual contribution can easily be added into eq 13. In this study, we simply regard their contribution to eq 13 as a second order effect and, thus, negligible. The computational aspects of the study (i.e., the KMC simulations) are used to test the validity of eq 13 for particular examples. For this purpose, the main approximation consists on considering a number of regions in the overall temperature f range, so that EM ij , E , and, especially, ωij are estimated as accurately as possible. We conclude that, from a theoretical perspective, eq 13 provides an excellent model to describe the origin of the apparent activation energy and, from a computational viewpoint, the estimates for EijM, Ef, and, especially, ωij become progressively more accurate as the number of regions is increased. Overall, we conclude that, even when a single process dominates, the macroscopic activation energy is a modified version of the microscopic activation energy of that process, due to the temperature-dependence of the corresponding multiplicity.

where θiϕ =

niϕ LxLy

(18)

is the density of empty sites of type i. Eq 18

θi Σj ≠ imijνij − Σj ≠ iθjmjiνji = Fθiϕ −

where Γ = ⟨R h⟩/Na. Thus, this study establishes that the average single-particle hop rate, Γ, loosely introduced in ref 23, should be understood as the time average of the total hop rate per particle. Although Γ has been determined in various ways,24−27 the exact formulation based on the time average has not been reported previously. In fact, our exact formulation provides a fundamental relation between the effective barrier Ea and the microscopic activation barriers (eq 13), unavailable before. Regarding the correlation factor itself, eq 1 provides a more fundamental formulation of f T, which is equivalent to the empirical definition24 (f T = 4DT/l2Γ), provided Γ is defined properly (e.g., as above). Should part of (or all) hop distances lij be different, eq 9 is generalized as follows:

dθi dt

(19)

Thus, by directly constraining the search space, these equations should be useful to improve the efficiency and accuracy of evolutionary procedures aimed at determining the hop rates.17 In this manner, extracting the reactivity directly from microscopy images might be possible in the future for both OSS and CVD processes. This might be useful as an automated alternative to existing approaches for extracting elementary hop rates.18,21,29,36

■

CONCLUSIONS In summary, we show that the tracer diffusivity, traditionally defined as DT = g⟨R2⟩/t = gl2f T⟨Nh⟩/t, is exactly described by the less familiar expression DT = gl2f T⟨R h⟩, where ⟨R h⟩ = ⟨Nh⟩/t is the time average of the total hop rate. On the basis of this, we are able to describe accurately the apparent activation energy of DT in terms of the microscopic activation barriers (eq 13). As a

∑ ⟨Wi ⟩⟨μij ⟩lij 2 i,j

is the occupancy, understood as

means that the density of adparticles of type i increases due to direct adsorption from the gas phase (Fiθϕi ) and/or due to hops by particles of type j ≠ i ending at sites of type i (Σj≠iθjmjiνji), while it decreases due to direct desorption to the gas phase (−Giθi) and/or due to hops by atoms of type i ending at sites of type j ≠ i (−θiΣj≠imijνij). Now, we note that during typical monolayer growth (at constant deposition and negligible desorption), the morphology of the system (e.g., the perimeter shapes) at any given coverage directly reflect a particular density of empty and occupied sites (θϕi and θi), as well as a particular set of values for the local multiplicities mij (i.e., the number of sites of type j surrounding the particles of type i). This opens the possibility for determining the hop rates (νij) by directly extracting the values of θi, θϕi and mij from two images separated by small increments in θi and solving

l2

1 2d

θi θ

dθi = Fiθiϕ + Σj ≠ iθjmjiνji − Giθi − θi Σj ≠ imijνij dt

Regarding eq 7, note that it can be rewritten as DT = fT 2d Γ

DT = fT

=

the probability of occupying a site of type i, and ⟨μij⟩ = ⟨mij⟩νij is the rateplicity (hop rate × local hop multiplicity). Both concepts have been defined in a recent study of the lowcoverage tracer diffusivity for complex energy landscapes with any number of distinct hop rates.32 It is remarkable that eq 17, which is valid out of equilibrium and contains all correlations between the hopping particles at any coverage, has the same form as eq 29 in ref 32 for noninteracting particles (θ ≈ 0) in equilibrium, where f T = 1. We also point out that, for the reverse process, i.e., anisotropic etching, we previously explained the origin of the apparent activation energy of the etch rate based on a similar formulation.40 Although future image processing techniques might enable analyzing sequences of briefly spaced, detailed microscopy images, eventually providing a reasonable description of the motion of all adparticles, we argue that a sequence of just two such images should contain essential information. For this purpose, let us consider the master equation solved in this study (i = 0, 1, ..., S)

describing νij, ⟨Mij⟩ and f T as proportional to e−Eijβ, e−Eij β, and

M

ni Na

l2

(17)

result, we show that the semiempirical expression DT = fT 2d Γ , 20320


Article

The Journal of Physical Chemistry C based on loosely defining Γ as an average single-particle hop rate in previous studies, becomes an exact formulation once Γ is understood as the time average of the total hop rate per particle, Γ = ⟨R h⟩/Na . Presumably, our formulation has eluded previous studies because (i) the tracer diffusivity is proportional to a slightly hidden version of the total hop rate (the timeaverage) and (ii) traditional studies have mainly focused on the density of islands, Nisl, and/or variables related to it, such as the θ−θ θ average island size s = N 0 ≈ N , the rate of growth of the isl

ds

average size ( dt ≈

F(1 − θ ) ), Nisl

Dimensional Crystals, and Hybrid Systems for Energy Conversion and Storage. Science 2015, 347, 1246501. (3) Cai, J.; Pignedoli, C. A.; Talirz, L.; Ruffieux, P.; Söde, H.; Liang, L.; Meunier, V.; Berger, R.; Li, R.; Feng, X.; et al. Graphene Nanoribbon Heterojunctions. Nat. Nanotechnol. 2014, 9, 896−900. (4) Ruffieux, P.; Wang, S.; Yang, B.; Sánchez-Sánchez, C.; Liu, J.; Dienel, T.; Talirz, L.; Shinde, P.; Pignedoli, C. A.; Passerone, D.; et al. On-Surface Synthesis of Graphene Nanoribbons with Zigzag Edge Topology. Nature 2016, 531, 489−492. (5) Bieri, M.; Nguyen, M.-T.; Gröning, O.; Cai, J.; Treier, M.; AïtMansour, K.; Ruffieux, P.; Pignedoli, C. A.; Passerone, D.; Kastler, M.; et al. Two-Dimensional Polymer Formation on Surfaces: Insight into the Roles of Precursor Mobility and Reactivity. J. Am. Chem. Soc. 2010, 132, 16669−16676. (6) Bronner, C.; Björk, J.; Tegeder, P. Tracking and Removing Br during the On-Surface Synthesis of a Graphene Nanoribbon. J. Phys. Chem. C 2015, 119, 486−493. (7) Liu, J.; Dienel, T.; Liu, J.; Groening, O.; Cai, J.; Feng, X.; Müllen, K.; Ruffieux, P.; Fasel, R. Building Pentagons into Graphenic Structures by On-Surface Polymerization and Aromatic Cyclodehydrogenation of Phenyl-Substituted Polycyclic Aromatic Hydrocarbons. J. Phys. Chem. C 2016, 120, 17588−17593. (8) Li, X.; Cai, W.; An, J.; Kim, S.; Nah, J.; Yang, D.; Piner, R.; Velamakanni, A.; Jung, I.; Tutuc, E.; et al. Large-Area Synthesis of High-Quality and Uniform Graphene Films on Copper Foils. Science 2009, 324, 1312−1314. (9) Lin, L.; Li, J.; Ren, H.; Koh, A. L.; Kang, N.; Peng, H.; Xu, H. Q.; Liu, Z. Surface Engineering of Copper Foils for Growing CentimeterSized Single-Crystalline Graphene. ACS Nano 2016, 10, 2922−2929. (10) Chen, C.-C.; Kuo, C.-J.; Liao, C.-D.; Chang, C.-F.; Tseng, C.-A.; Liu, C.-R.; Chen, Y.-T. Growth of Large-Area Graphene Single Crystals in Confined Reaction Space with Diffusion-Driven Chemical Vapor Deposition. Chem. Mater. 2015, 27, 6249−6258. (11) Geng, D.; Wang, H.; Yu, G. Graphene Single Crystals: Size and Morphology Engineering. Adv. Mater. (Weinheim, Ger.) 2015, 27, 2821−2837. (12) Kobayashi, T.; Bando, M.; Kimura, N.; Shimizu, K.; Kadono, K.; Umezu, N.; Miyahara, K.; Hayazaki, S.; Nagai, S.; Mizuguchi, Y.; et al. Production of a 100-m-Long High-Quality Graphene Transparent Conductive Film by Roll-to-Roll Chemical Vapor Deposition and Transfer Process. Appl. Phys. Lett. 2013, 102, 023112. (13) Novoselov, K. S.; Fal’ko, V. I.; Colombo, L.; Gellert, P. R.; Schwab, M. G.; Kim, K. A Roadmap for Graphene. Nature 2012, 490, 192−200. (14) Gong, C.; Colombo, L.; Cho, K. Photon-Assisted CVD Growth of Graphene Using Metal Adatoms As Catalysts. J. Phys. Chem. C 2012, 116, 18263−18269. (15) Zhang, Z.; Chen, X.; Lagally, M. G. Bonding-Geometry Dependence of Fractal Growth on Metal Surfaces. Phys. Rev. Lett. 1994, 73, 1829−1832. (16) Cox, E.; Li, M.; Chung, P.-W.; Ghosh, C.; Rahman, T. S.; Jenks, C. J.; Evans, J. W.; Thiel, P. A. Temperature Dependence of Island Growth Shapes During Submonolayer Deposition of Ag on Ag(111). Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 71, 115414. (17) Ferrando, N.; Gosalvez, M. A.; Ayuela, A. Evolutionary Kinetic Monte Carlo: Atomistic Rates of Surface-Mediated Processes from Surface Morphologies. J. Phys. Chem. C 2014, 118, 11636−11648. (18) Venables, J. A.; Spiller, G. D. T.; Hanbucken, M. Nucleation and Growth of Thin Films. Rep. Prog. Phys. 1984, 47, 399. (19) Amar, J. G.; Family, F.; Lam, P.-M. Dynamic Scaling of the Island-Size Distribution and Percolation in a Model of Submonolayer Molecular-Beam Epitaxy. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 8781−8797. (20) Ratsch, C.; Venables, J. A. Nucleation Theory and the Early Stages of Thin Film Growth. J. Vac. Sci. Technol., A 2003, 21, S96− S109. (21) Evans, J.; Thiel, P.; Bartelt, M. Morphological Evolution During Epitaxial Thin Film Growth: Formation of 2D Islands and 3D Mounds. Surf. Sci. Rep. 2006, 61, 1−128.

isl

etc.. These variables are defined

only between island nucleation and island coalescence, thus reducing the coverage window where the studies can be performed. On the other hand, the total hop rateand the tracer diffusivityare properly defined for any coverage value, well before island nucleation and well after island coalescence. Thus, we conclude that the tracer diffusivity is the natural transport coefficient that describes monolayer growth at constant deposition flux. Finally, we stress the idea that the actual shapes of the clusters and islands at any given coverage for both OSS and CVD processes directly reflect particular values for the density of empty and occupied sites, as well as for the numbers of sites of type j surrounding the particles of type i. This should help extracting the hop rates directly from microscopy images in the future.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b05794. Details of the KMC simulations, Table S1 (atomistic activation energies used in the simulations), Figure S1 (equivalent to Figure 1, now for the square lattice), Figure S2 (coverage-dependence of the correlation factor for the triangular and square lattices), Figure S3 (complementary to Figure 2, now for additional coverage values), and Figure S4 (complementary to Figure 4, now for nine temperature subregions) (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*(M.A.G.) E-mail: [email protected]. ORCID

Miguel A. Gosálvez: 0000-0002-7621-0755 Joseba Alberdi-Rodriguez: 0000-0002-9078-8553 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge support by the Ramón y Cajal Fellowship Program by the Spanish Ministry of Science and Innovation, and the 2015/01 contract by the DIPC. The KMC calculations were performed on the ATLAS supercomputer at the DIPC.

■

REFERENCES

(1) Brownson, D. A. C.; Kampouris, D. K.; Banks, C. E. An Overview of Graphene in Energy Production and Storage Applications. J. Power Sources 2011, 196, 4873−4885. (2) Bonaccorso, F.; Colombo, L.; Yu, G.; Stoller, M.; Tozzini, V.; Ferrari, A. C.; Ruoff, R. S.; Pellegrini, V. Graphene, Related Two20321


Article

The Journal of Physical Chemistry C (22) Körner, M.; Einax, M.; Maass, P. Capture Numbers and Island Size Distributions in Models of Submonolayer Surface Growth. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 085403. (23) Reed, D. A.; Ehrlich, G. Surface Diffusion, Atomic Jump Rates and Thermodynamics. Surf. Sci. 1981, 102, 588−609. (24) Gomer, R. Diffusion of Adsorbates on Metal Surfaces. Rep. Prog. Phys. 1990, 53, 917−1002. (25) Hjelt, T.; Vattulainen, I.; Merikoski, J.; Ala-Nissila, T.; Ying, S. Dynamical Mean Field Theory: An Efficient Method to Study Surface Diffusion Coefficients. Surf. Sci. 1998, 402-404, 253−256. (26) Vattulainen, I.; Ying, S. C.; Ala-Nissila, T.; Merikoski, J. Memory Effects and Coverage Dependence of Surface Diffusion in a Model Adsorption System. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 7697−7707. (27) Ala-Nissila, T.; Ferrando, R.; Ying, S. Collective and Single Particle Diffusion on Surfaces. Adv. Phys. 2002, 51, 949−1078. (28) Danani, A.; Ferrando, R.; Scalas, E.; Torri, M. Lattice-Gas Theory of Collective Diffusion in Adsorbed Layers. Int. J. Mod. Phys. B 1997, 11, 2217−2279. (29) Tringides, M. C.; Hupalo, M. Surface Diffusion Experiments with STM: Equilibrium Correlations and Non-Equilibrium Low Temperature Growth. J. Phys.: Condens. Matter 2010, 22, 264002. (30) Uebing, C.; Gomer, R. The Diffusion of Oxygen on W(110) The Influence of the p(2 × 1) Ordering. Surf. Sci. 1997, 381, 33−48. (31) Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107−115. (32) Gosál vez, M. A.; Otrokov, M. M.; Ferrando, N.; Ryabishchenkova, A. G.; Ayuela, A.; Echenique, P. M.; Chulkov, E. V. Low-Coverage Surface Diffusion in Complex Periodic Energy Landscapes: Analytical Solution for Systems with Symmetric Hops and Application to Intercalation in Topological Insulators. Phys. Rev. B: Condens. Matter Mater. Phys. 2016, 93, 075429. (33) Voter, A. Radiation Effects in Solids. NATO Science Series 2007, 235, 1−23. (34) Chatterjee, A.; Vlachos, D. G. An Overview of Spatial Microscopic and Accelerated Kinetic Monte Carlo Methods. J. Comput.-Aided Mater. Des. 2007, 14, 253−308. (35) Jansen, A. P. J. An Introduction to Kinetic Monte Carlo Simulations of Surface Reactions. Lect. Notes Phys. 2012, 856, 1. (36) Bott, M.; Hohage, M.; Morgenstern, M.; Michely, T.; Comsa, G. New Approach for Determination of Diffusion Parameters of Adatoms. Phys. Rev. Lett. 1996, 76, 1304−1307. (37) Trushin, O.; Karim, A.; Kara, A.; Rahman, T. S. Self-Learning Kinetic Monte Carlo Method: Application to Cu(111). Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 115401. (38) Nandipati, G.; Kara, A.; Shah, S. I.; Rahman, T. S. Kinetically Driven Shape Changes in Early Stages of Two-Dimensional Island Coarsening: Ag/Ag(111). Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 115402. (39) Shah, S. I.; Nandipati, G.; Karim, A.; Rahman, T. S. SelfLearning Kinetic Monte Carlo Simulations of Self-Diffusion of Small Ag Islands on the Ag(111) Surface. J. Phys.: Condens. Matter 2016, 28, 025001. (40) Gosálvez, M. A.; Cheng, D.; Nieminen, R. M.; Sato, K. Apparent Activation Energy during Surface Evolution by Step Formation and Flow. New J. Phys. 2006, 8, 269.

20322


Microscopic Origin of the Apparent Activation ... - ACS Publications

Recommend Documents