Kinetic Monte Carlo Model for Gas Phase Diffusion in Nanoscopic

Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

pubs.acs.org/JPCC

Kinetic Monte Carlo Model for Gas Phase Diffusion in Nanoscopic Systems Elisabeth M. Dietze and Philipp N. Plessow* Institute of Catalysis Research and Technology (IKFT), Hermann-von-Helmholtz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany S Supporting Information *

ABSTRACT: Transport of atoms and molecules via the gas phase plays an important role in many processes in heterogeneous catalysis. Macroscopic diffusion, for example, in reactors, is typically modeled with continuum models. Much smaller length scales are involved if diffusion occurs between nanoparticles. One such example is a sintering mechanism, where volatile PtO2 mediates mass transfer between Pt particles. We developed a kinetic Monte Carlo model that explicitly simulates the kinetics of single atoms or molecules in the gas phase that result from collisions with a background gas. This model accurately reproduces ideal gas properties such as the diffusion constant. In model applications, we study gas-phase-mediated mass transfer as a function of the distance between the involved surfaces. If these distances are within the mean free path, typically a micrometer or lower, continuum models based on Fick’s laws deviate from the explicit simulation. This can be explained by the low number of collisions that occur if the length scale of diffusion is not significantly larger than the mean free path in the gas phase.

1. INTRODUCTION Most industrially relevant catalytic processes are heterogeneous. Catalytic reactions take place on the surfaces of the catalyst, for example, transition-metal nanoparticles or porous zeolites. Reactants and products are typically gases at the reaction conditions. This means that mass transport of the reactants to the catalyst and of the products away from the catalyst is crucial. Therefore, a lot of work is devoted to understanding and improving how the gas is transported through reactors by simulation. Simulations for reactors (computational fluid dynamics) typically employ continuum models that do not directly incorporate the collision processes in the gas phase that lead to convection and diffusion.1−5 A guide to an appropriate model to describe gas flow is given by the Knudsen number Kn.6 It is defined as Kn = λ/L, the ratio of the mean free path λ to the characteristic length L. For Kn < 0.1, continuum models are valid, and for Kn > 10, molecular models. In the case of sintering nanoparticles, Kn is in the transition regime.7 One example, in which diffusion on the nanometer scale may be important is the sintering of Pt nanoparticles in an oxidizing atmosphere.8 Besides particle migration9−11 and surfacemediated Ostwald ripening,12−19 another potential mechanism for sintering is ripening that is mediated by volatile species, PtO2.20−26 Experiments have shown how sintering depends on the gas flow and the geometrical setup of the reactor for the sintering of Pt particles in an oxidizing atmosphere,27−32 which is indicative of ripening via the gas phase. In recent simulations, we have shown that ripening through the gas phase is indeed a realistic mechanism.33 The employed mean-field model © XXXX American Chemical Society

altogether neglects mass-transport limitations in the gas phase. It is assumed that the sintering kinetics is controlled exclusively by the rate with which PtO2 is emitted from the particles. This is equivalent to assuming a constant background pressure of PtO2. Besides neglecting diffusion limitations, a main drawback of this model is that it cannot describe how much PtO2 is not redeposited onto Pt nanoparticles and is therefore lost. In the actual sintering process, PtO2 is emitted and (potentially) readsorbed by the supported nanoparticles. The distance between the particles is often on the order of a few ten nanometers. It is therefore not unlikely that, after being emitted, a PtO2 molecule will collide with another nanoparticle without being scattered in the gas phase even once. This effect will not be observed in a continuum simulation that is based only on the diffusion constant of PtO2(g). Being able to incorporate these effects is the motivation for the kinetic Monte Carlo (kMC) model described in this work. It is worth noting that the simulation of surface-mediated Ostwald ripening kinetics faces related problems. Usually, surface-mediated ripening is described by assuming a constant concentration of migrating species on the surface, equivalent to neglecting mass-transport limitations on the surface.18,34,35 This is similar to the mentioned mean-field model in the gas phase. It was shown, in various examples, that local effects such as defects and diffusion barriers lead to a breakdown17,36−38 of the Received: February 22, 2018 Revised: April 6, 2018

A

DOI: 10.1021/acs.jpcc.8b01816 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 1. Inverted probability function of (a) vb(ρvb(vp)) and (b) θ(ρθ(vs)) for two different example velocities vp and parameters vs, respectively.

simply referred to as particles) and a background gas. The particles do not interact with each other and do not change the properties of the background gas. This is expected to be a good approximation when the particles are present at a very low concentration. This will be the case in the applications of interest, where mass transport between solid particles is mediated by species with low vapor pressure while the gas phase is mainly composed of inert gas, air, or a reactive atmosphere. It is also a valid approach for the simulation of pure gas-phase properties, where only a single particle is simulated. This will be the case that we first study, the diffusion of 40Ar in 40Ar. The only kMC processes are collisions between particles and the background gas, in between which the particles move with a constant velocity. Initially, particles are supplied with certain velocities and from there on move until collisions change both direction and magnitude of the velocity. The probability for collision of the particle with an (implicitly present) particle of the background gas is proportional to the relative velocity of these two. For a particle with a given velocity vp, the probability to collide with a particle from the background gas thus depends both on its absolute velocity vb and the relative orientation. The probability function for the collision, integrated from 0 to vb, is a continuous, smooth function of vb and follows from the Maxwell−Boltzmann (MB) distribution f(v):

mean-field approximation. One possibility to take into account local effects is to extend the mean-field approximation to a nearest neighbor or local correlation approach.37 In this work, we develop a kMC model that describes the collision of an explicitly simulated molecule within a uniform background of ideal gas. Similar models were developed to describe high-power impulse magnetron sputtering and direct current sputtering.39 To calculate the flow of dilute gases, often the direct simulation Monte Carlo (DSMC) method after Bird7 is used. The DSMC is a probabilistic approach, where the number of simulated particles is not equal to the number of real particles. In consequence, it is not suitable to model the pathway of a single molecule influenced by a background gas. In our approach, the background gas is not explicitly simulated but interacts with the molecule through collisions that follow from its temperature, its mass distribution, collision cross section, and collective gas flow. We will first describe the kMC model of gas phase diffusion and test it by reproducing ideal gas parameters for 40Ar. As a second step, two model systems are investigated to compare the kMC model with analytical solutions of Fick’s laws.

2. METHODS To simulate the diffusion of a particle through an ideal gas, we use a kMC approach with a Bortz−Kalos−Lebowitz (BKL) algorithm.40 In a conventional BKL algorithm, all possible processes N of the system state are first identified and the rate constants ki are calculated for each process. The sum over all rate constants gives the total rate constant ktot:

v

ρv (v b) =

(4)

where

∑ ki

z(vp , v b) dv = f (v b)nbσ vr̅ dv

(1)

i=1

q

q−1

∑ ki ≥ ρ1k tot ≥

∑ ki

i=1

i=1

(2)

z(vp) =

The system time is advanced by

∫0

∞

z(vp , v b) dv b

(6)

This expression already includes all possible orientations of vp and vb. The inverted probability function is shown in Figure 1a, where vb is plotted against ρvb(vp). In contrast to usual kMC algorithms, which have a finite number of processes, here, we have an infinite number of possible processes. This is due to the continuum nature of the background gas, where at any time, all velocities follow the MB distribution. Figure 2 illustrates ρvb(vp),

ln ρ2 k tot

(5)

with the collision cross section σ, the mean relative velocity vr̅ averaged over all angles θ, and the molecular density of the background gas nb (for more details, see the Supporting Information). The collision frequency of the particle Mp is

The executed process q is chosen with a random number ρ1 ∈ ]0, 1] according to

Δt = −

z(vp)

p

N

k tot =

∫0 b z(vp , v b) dv

(3)

with ρ2 ∈ ]0, 1] being a second random number. In an ideal gas, particles interact only through collisions and otherwise move linearly at a constant velocity. Our simulation approach contains explicitly simulated particles (from now on B


Article

The Journal of Physical Chemistry C which is a smooth function (right side), whereas a finite number of rate constants would lead to a steplike curve (left side).

θ = θ(ρv = ρ1b )

(10b)

s

with ρ1a and ρ1b ∈ ]0, 1] being random numbers. The resulting particle velocity after the collision process is determined according to a hard-sphere collision. The change in time Δt is calculated according to eq 3, and the position of the particles in space up to the collision process is updated:

Δx = Δt ·vp

The precision of a kMC simulation depends on the interpolation procedure. In our case, the precision depends on the number of points Nvp to precompute the probability distribution of vb and Nint, the number of points to precompute the integral in eq 4. To investigate the precision of the simulations, gas-phase properties such as the diffusion coefficient D and the mean velocity v:=⟨|v|⟩ were studied. ̅ The kMC simulation, with NkMC steps, is repeated Ntraj times. The particle is initialized with a velocity which has the highest probability according to the Boltzmann distribution. To allow the system to equilibrate, Dsample and vs̅ ample are calculated as follows for each trajectory after 100 initial steps (n > 100):

Figure 2. Visualization of the transition from a finite number of processes (left side) to a continuum of processes (right side).

The direction of the velocity vb can be described with two angles θ and ϕ, using polar coordinates. The angle between vp and vb is given by θ (see Figure S1), whereas ϕ describes different orientations of vb with respect to the direction of vr. In the absence of an external field, all values of ϕ have an equal probability. Therefore, ϕ is chosen randomly in the range [0,2π]. The probability for collision with a certain angle θ is given by 3/2

(vs − cos θ )

ρv (θ ) =

3/2

(vs + 1)

s

vsample = ̅

− (vs − 1)

vp2 + v b 2 2vpv b

tkMC

(12)

∑i = 1 xn , i 2 n

3·2·∑k = 1 Δtk

(13)

with n being the current step and tkMC being the elapsed time excluding the 100 initial steps. vn, Δtn, and xn,i are the mean velocity, change in time, and the ith-component of the particle position of each step n of the simulation, respectively. Dsim and vs̅ im are calculated as the average over all simulated samples. As a model system, the self-diffusion of 40Ar is studied at a temperature of T = 500 K and a pressure of p = 1 bar. The diffusion coefficient Dref(40Ar) = 3.388 × 10−5 ± 4 × 10−8 m2/s was calculated according to the textbook of Kennard41 (see also the Supporting Information). The approximations in common textbooks deviate by 44% (Atkins and de Paula42) or by 17% (Wedler43) (see the Supporting Information). vr̅ ef(40Ar) = 514.69 m/s follows from the ideal gas law. Figures 3 and 4 show the decreasing deviation of Dsim and vs̅ im from the reference values with increasing Nint and Nvp,

(7)

Here, the collision probability is rewritten in terms of the dimensionless parameter vs:

vs =

∑n vn·Δtn

3

Dsample =

3/2

− (vs − 1)3/2

(11)

(8)

The outcome of the scattering process depends on vb and the angle θ. Because of its complicated functional form, vb(ρvb(vp)) was precomputed before the kMC simulation. For an equally spaced grid of Nvp values, vb(ρvb) was interpolated between ρvb = 0 and 1 using cubic splines. For a continuous spline, a minimal distance dmin between the values of ρvb was introduced. The results are independent of dmin between dmin = 10−15 and 10−6 (see the Supporting Information). During the kMC calculation, for a given vp and a random number ρ1a, vb was obtained by quadratic interpolation of vb(ρvb(vp) = ρ1a) using the three closest vp values of the precomputed grid. As shown in eq 7, θ only depends on vs and the probability ρvs. Equation 7 can therefore be solved for θ: θ = arccos( [−ρv (θ )((vs + 1)3/2 − (vs − 1)3/2 ) + (vs − 1)3/2 ]2/3 s

+ vs)

(9)

The collision process with a background molecule with a certain velocity and orientation is then chosen according to v b = v b(ρv = ρ1a ) p

Figure 3. Deviation of Dsim and vs̅ im from Dref and vr̅ ef with respect to Nint and constant Nvp = 2000 (NkMC = 3 × 103, Ntraj = 4 × 106, and number of calculations per point: 10; the error bars show the standard deviation of the simulated data with respect to the mean value).

(10a) C


Article


Figure 4. Deviation of Dsim and vs̅ im from Dref and vr̅ ef with respect to Nvp and constant Nint = 2 × 105 (NkMC = 3 × 103, Ntraj = 4 × 106, and number of calculations per point: 10; the error bars show the standard deviation of the simulated data with respect to the mean value).

Figure 6. Comparison of Dsim (black circles) and Dref (black line) and vs̅ im (blue squares) and vr̅ ef (blue line) for different temperatures at p = 1 bar (Nvp = 2000, Nint = 2 × 105, NkMC = 3 × 103, Ntraj = 105, and number of calculations per point: 20).

respectively. Taking the value of the last data point in Figure 4, the simulated diffusion coefficient is Dsim = 3.3853 × 10−5 m2/s, which deviates by 0.09% from Dref and is within the error range given by Kennard. The simulated velocity is vs̅ im = 514.70 m/s, which deviates by 0.002% from vr̅ ef. The The accuracy is expected to be sufficient for sintering calculations and could be further reduced by increasing Nint and Nvp. Figure 5 shows the differential scattering cross section of 40Ar for a constant collision energy of 0.062 eV at T = 500 K and

steady-state solution based on Fick’s laws can be obtained, thus providing a continuum model for comparison. 3.1. Parallel Surfaces. The first example is the diffusion of single atoms between two parallel surfaces with a surface area A and a distance Δz along the z-axis. The system is infinitely large in two dimensions and is described with a quadratic unit cell with the surface area A. Ultimately, A does not affect the outcome of the simulation (per a given surface area). However, A controls the rate constant for emission and adsorption for a given surface area. Choosing a very small surface area of A = 10−16 m2 results in a very small gas-phase volume in the unit cell. This, in turn, means that during the simulation, generally only one particle is present in the gas phase, which greatly simplifies the analysis of the flux, as opposed to a situation with a larger volume and a larger number of particles in the gas phase. From both surfaces, atoms can be emitted with rate constants K1,2 according to their chemical potentials μ1,2. In the case of surface collisions, atoms are adsorbed with a probability that is given by the sticking coefficient S. We will study the diffusion of Pt(g) in between clean Pt surfaces and in a pure 40Ar gas phase. In this model, S = 1 is chosen as is appropriate for Pt(g). More complex surface−gas interactions could be incorporated through an angle- and velocity-dependent sticking coefficient. The rate equation for atom emission is defined33 as

Figure 5. Differential scattering cross section for a collision energy of 0.062 eV (Nvp = 2000, Nint = 2 × 105, NkMC = 3 × 103, and Ntraj = 106).

p = 1 bar, which qualitatively agrees with Phelps et al.44 The probability to have either small or high scattering angles is increased, whereas the probability to find scattering angles between 20° and 160° is orders of magnitude smaller. Oscillations are due to poor sampling in the range where collisions are less likely. Figure 6 shows Dsim and Dref and the corresponding vs̅ im and vr̅ ef for different temperatures at a constant pressure of p = 1 bar. Both graphs are in excellent agreement over the whole temperature range. We have thus shown that ideal gas properties can be simulated accurately using the developed kMC model with the parameters Nvp = 2000 and Nint = 2 × 105, reproducing the mean velocity and the diffusion coefficient over a wide range of temperatures.

K1,2 =

μ1,2 − ΔG° A·S p° exp kBT 2π ·m ·kBT

(14)

with p° being the standard pressure and ΔG° being the standard Gibbs free energy in the gas phase. The atoms are emitted with a velocity distribution for vz corresponding to that for collision of the ideal gas with a surface with a normal vector in the z-direction so that emission and collision are in equilibrium at the interface. Importantly, this velocity distribution is not the 1D MB-distribution, ρ(vz), but proportional to vz·ρ(vz) (Figure 7).45−47 vx and vy are chosen according to the 1D MB distribution. For a system in equilibrium, the potential difference of the walls is Δμ = 0 eV. In consequence, all gas properties should be constant throughout the system. In Figure 8, the simulated mean squared velocity vrms,z of the z-component is plotted and is in good agreement with the theoretical value vref rms,z = 291.40 m/s at T = 2000 K. The calculated particle density (PD)

3. MODEL SYSTEMS To test the kMC model for diffusion in systems on a nanoscale, two model systems were chosen. In these systems, an analytical D


Article


reflected, depending on S. The positions of all other atoms in the gas phase are then updated using Δt. Formally, this kMC approach can be considered as a first-reaction method,49 where possible processes are collision with the surfaces, collision in the gas phase, and emission from the surfaces. Here, the wall collision processes are deterministic. The probability for any other process to occur follows from the BKL algorithm described above, using the total rate constants of these processes and eq 3. The total flux of platinum atoms Jnum total is calculated as num Jtotal =

Figure 7. Probability distribution vz·ρ(vz) at T = 2000 K for Pt atoms.

(NE1 − NC1) − (NE2 − NC2) 2·A ·tkMC

(16)

with NE(1,2) and NC(1,2) being the number of emitted and adsorbed particles on the surfaces, respectively. T = 2000 K was chosen as a temperature at which the formation of Pt(g) is relevant. The total flux for different pressures and distances at T = 2000 K is shown in Figure 9.

Figure 8. Bars: simulated PD for Δμ = 0. Red line: PD in thermodynamic equilibrium. Blue line with symbols: deviation of the simulated mean squared velocity from vref rms,z = 291.40 m/s (kMC parameters: Δz = 10−4 m, Nvp = 2000, Nint = 105, NkMC = 2 × 107, Ntraj = 25, T = 2000 K, and p = 1 bar).

remains constant over the whole range, in agreement with the one given by the thermodynamic equilibrium. Next, the chemical potential μ1 = −5.82 eV48 of the left surface was kept constant at the value of solid Pt. μ2 was chosen to be variable with the condition μ2 > μ1, causing a flux from the right surface to the left. The solution of Fick’s laws gives the following expression for the total flux Jtotal = Jout − Jin, which generally consists of the out-coming and in-coming particle flux: Jtotal = Δn·α =

−DΔneq αΔz + 2D

·α

Figure 9. Total flux between two parallel surfaces for different pressures and surface distances at T = 2000 K and Δμ = 100 meV. Solid lines with symbols: simulated flux. Solid lines: flux according to eq 15 (kMC parameters: Nvp = 2000, Nint = 105, NkMC = 109, and Ntraj = 10).

The solid lines represent the solution of eq 15, and the lines with point symbols represent the simulated values. The general behavior of both curves is similar and so is the effect of the background pressure. For distances Δz > 10−4 m for p = 10 bar or increased Δz for smaller pressures, both models give the same result and converge to zero for large distances because of a negligible concentration gradient. For distances Δz < 10−4 m, the kMC model shows a higher flux than the continuum model. In our case, the simulated flux is twice the flux given by the continuum model for Δz < 5 × 10−9 m. On this length scale, there is an increasing probability for particles to be emitted and adsorbed on the other wall without scattering processes. Carrillo et al.50 investigated experimentally how flowing air reduces the amount of Pt in a Pt/SiO2 sample, by studying the Pt particle size distribution. They concluded that the flowing air is sweeping volatile PtO2(g) off the surface. In our kMC model, the flowing air can be simulated by applying a drift velocity vdrift to the background gas. We simulate how volatile Pt(g) is swept off the surface. In order to achieve significant pressures of Pt(g), we employ again T = 2000 K, meaning a higher temperature than necessary to cause emission of PtO2(g), and apply a drift velocity in the x-direction vdrift,x. To find out for which distance range from a surface the Pt atoms adapt to the

(15)

−1/2

with α = kBT(2π·mkBT) and Δneq = (p2,eq − p1,eq)/kBT, with pi,eq being the equilibrium pressure at which surface i would be in equilibrium with the gas phase. We will now discuss how the kMC algorithm for the simulation of a pure gas phase is extended to include surface processes. In addition to the gas-phase scattering, the emission of an atom and the adsorption/reflection on the surfaces are considered as kMC processes. In the absence of gas-phase collisions, the atoms move linearly with a constant velocity until they collide with a surface and are either adsorbed or reflected. In the presence of an external potential (not considered in this work), the trajectory could be modified accordingly. The time to collision with a background-gas particle, Δt, is computed first according to eq 3 and is used to check if any atom in the gas phase collides with one of the surfaces before gas-phase collision. In the case of a collision within Δt, Δt is changed to the actual collision time Δtcol and the atom is adsorbed or E


Article


is chosen. The high Δμ gives an increased ratio of the probabilities of emissions from the inner sphere and outer sphere. This reduces the computational cost as the fraction of executed processes that are relevant for diffusion increases. The radius r1 of the inner sphere is 1 nm. In general, the Gibbs− Thomson equation51 could be used to estimate a radiusdependent chemical potential. The emission rate is calculated according to eq 14. As before, the velocity distribution is chosen to be in equilibrium with the gas phase. In contrast to the parallel surfaces, the surface areas of the spheres are not equal but depend on the radius. In analogy to a system involving real nanoparticles, we study the flux from the inner sphere J(r1), mimicking a sintering nanoparticle. From Fick’s laws, it follows

drift velocity, we show in Figure 10 the ratio of vx to vdrift,x for vdrift,x = 10, 100 and 1000 m/s. It can be seen that vx converges

J(r1) = α(neq r1 − nr1)

(17)

with Figure 10. Convergence of the ratio of the calculated velocity vx to the applied drift velocity vdrift,x at T = 2000 K, p = 1 bar, Δz = 10−4 m and Δμ = 0 eV, and vdrift,x = 10, 100, 1000 m/s (kMC parameters: Nvp = 2000, Nint = 105, NkMC = 2 × 107, and Ntraj = 50).

r1

n =

neq r1αr12r2 − neq r1αr1r2 2 − Dneq r1r12 − Dneq r2r2 2 αr12r2 − αr1r2 2 − Dr12 − Dr2 2

(18)

Equation 18 shows that J(r1) depends on not only r1 and μ1 but also r2 and μ2. With increasing r2, the surface area of the outer sphere increases quadratically. This generally leads to a faster adsorption of particles emitted from the inner sphere. Eventually, at large r2, the system resembles an open system where the background pressure of Pt(g) is determined by the outer sphere, where it follows from the chemical potential of solid platinum. The limit given by the open system is plotted as a dashed black line in Figure 12.

within 5 × 10−6 m to the applied drift velocity, independent of its magnitude. The convergence towards the drift velocity is important to estimate the rate with which Pt atoms are swept off the surface. In Figure 11, the simulated PD is shown as bars. It is in good agreement with the PD according to Fick’s law.

Figure 11. Bars: simulated PD for different distances from the emitting surface. Red line: PD according to Fick’s law at T = 2000 K, p = 1 bar, Δz = 10−4 m, and Δμ = 4 eV (kMC parameters: Nvp = 2000, Nint = 105, NkMC = 2 × 107, and Ntraj = 50).

Figure 12. Total flux between two centered spheres for different pressures and surface distances at T = 2000 K, r1 = 1 nm, and μ1 = −1.82 eV (kMC parameters: Nvp = 2000, Nint = 105, NkMC = 4 × 108, and Ntraj = 20).

3.2. Centered Spheres. As a second example, we study the diffusion of single particles between two centered spheres. Diffusion occurs through the gas phase between the inner sphere with radius r1 and the outer sphere with radius r2. The distance between the spheres is defined as Δz = r2 − r1. When the volume between the spheres becomes small, the gas phase consists in reality only of a small discrete number of molecules. The implicit background model does not provide an accurate description of this situation. However, because the number of collisions is small in this case anyway, we do not think that this is a significant problem, given that the situation of very small closed systems is not very common. Again, Pt is the surface material and 40Ar constitutes the gas phase. The chemical potential of the outer sphere is kept constant at μ2 = −5.82 eV, and for the inner one, μ1 = −1.82 eV

A competing effect is that, with increasing r2, diffusion limitations, which depend on the pressure, become more important. Figure 12 shows the flux from a sphere with r1 = 1 nm at T = 2000 K, different pressures, and surface distances. The solid lines show the flux according to eq 17. For 1 bar, the flux increases with increasing Δz and converges to the flux of an open system, which is also the case for smaller pressures. For 50 and 100 bar, the flux according to Fick’s laws is increasing and converges to a constant value for Δz > 10−7 m, which is smaller than the flux of the open system. From the kMC simulations, the flux of the inner sphere is calculated as F


Article

The Journal of Physical Chemistry C J num (r1) =

(NE1 − NC1) A ·tkMC

(3) Schaefer, C.; Jansen, A. P. J. Coupling of kinetic Monte Carlo simulations of surface reactions to transport in a fluid for heterogeneous catalytic reactor modeling. J. Chem. Phys. 2013, 138, 054102. (4) Matera, S.; Blomberg, S.; Hoffmann, M. J.; Zetterberg, J.; Gustafson, J.; Lundgren, E.; Reuter, K. Evidence for the Active Phase of Heterogeneous Catalysts through In Situ Reaction Product Imaging and Multiscale Modeling. ACS Catal. 2015, 5, 4514−4518. (5) Majumder, D.; Broadbelt, L. J. A multiscale scheme for modeling catalytic flow reactors. AIChE J. 2006, 52, 4214−4228. (6) Knudsen, M. Die Gesetze der Molekularströmung und der inneren Reibungsströmung der Gase durch Röhren. Ann. Phys. 1909, 333, 75−130. (7) Bird, G. Molecular Gas Dynamics and the Direct Simulation of Gas Flows; Oxford Science Publications: New York, 1994. (8) Nie, L.; Mei, D.; Xiong, H.; Peng, B.; Ren, Z.; Hernandez, X. I. P.; DeLaRiva, A.; Wang, M.; Engelhard, M. H.; Kovarik, L.; Datye, A. K.; Wang, Y. Activation of surface lattice oxygen in single-atom Pt/CeO2 for low-temperature CO oxidation. Science 2017, 358, 1419−1423. (9) Li, L.; Plessow, P. N.; Rieger, M.; Sauer, S.; Sánchez-Carrera, R. S.; Schaefer, A.; Abild-Pedersen, F. Modeling the Migration of Platinum Nanoparticles on Surfaces Using a Kinetic Monte Carlo Approach. J. Phys. Chem. C 2017, 121, 4261−4269. (10) Zhdanov, V. P. Kinetics of Migration and coalescence of supported nm-sized metal particles. Surf. Rev. Lett. 2008, 15, 217−220. (11) Sault, A. G.; Tikare, V. A new Monte Carlo model for supported-catalyst sintering. J. Catal. 2002, 211, 19−32. (12) Zhdanov, V. P.; Larsson, E. M.; Langhammer, C. Novel aspects of Ostwald ripening of supported metal nanoparticles. Chem. Phys. Lett. 2012, 533, 65−69. (13) Larsson, E. M.; Millet, J.; Gustafsson, S.; Skoglundh, M.; Zhdanov, V. P.; Langhammer, C. Real Time Indirect Nanoplasmonic in Situ Spectroscopy of Catalyst Nanoparticle Sintering. ACS Catal. 2012, 2, 238−245. (14) Behafarid, F.; Cuenya, B. R. Towards the Understanding of Sintering Phenomena at the Nanoscale: Geometric and Environmental Effects. Top. Catal. 2013, 56, 1542−1559. (15) Challa, S. R.; Delariva, A. T.; Hansen, T. W.; Helveg, S.; Sehested, J.; Hansen, P. L.; Garzon, F.; Datye, A. K. Relating rates of catalyst sintering to the disappearance of individual nanoparticles during Ostwald ripening. J. Am. Chem. Soc. 2011, 133, 20672−20675. (16) Zhdanov, V. P.; Schweinberger, F. F.; Heiz, U.; Langhammer, C. Ostwald ripening of supported Pt nanoclusters with initial size-selected distributions. Chem. Phys. Lett. 2015, 631−632, 21−25. (17) Ouyang, R.; Liu, J.-X.; Li, W.-X. Atomistic theory of ostwald ripening and disintegration of supported metal particles under reaction conditions. J. Am. Chem. Soc. 2013, 135, 1760−1771. (18) Simonsen, S. B.; Chorkendorff, I.; Dahl, S.; Skoglundh, M.; Sehested, J.; Helveg, S. Ostwald ripening in a Pt/SiO2 model catalyst studied by in situ TEM. J. Catal. 2011, 281, 147−155. (19) Lo, A.; Skodje, R. T. Kinetic and Monte Carlo models of thin film coarsening: Cross over from diffusion-coalescence to Ostwald growth modes. J. Chem. Phys. 2000, 112, 1966. (20) Wynblatt, P.; Gjostein, N. A. Particle growth in model supported metal catalystsI. Theory. Acta Metall. 1976, 24, 1165−1174. (21) Wynblatt, P. Particle Growth in Model Supported Comparison of Experiment. Acta Metall. 1976, 24, 1175−1182. (22) Datye, A. K.; Xu, Q.; Kharas, K. C.; McCarty, J. M. Particle size distributions in heterogeneous catalysts: What do they tell us about the sintering mechanism? Catal. Today 2006, 111, 59−67. (23) Porsgaard, S.; Merte, L. R.; Ono, L. K.; Behafarid, F.; Matos, J.; Helveg, S.; Salmeron, M.; Cuenya, B. R.; Besenbacher, F. Stability of platinum nanoparticles supported on SiO2/Si(111): a high-pressure X-ray photoelectron spectroscopy study. ACS Nano 2012, 6, 10743− 10749. (24) Adibi, P. T. Z.; Zhdanov, V. P.; Langhammer, C.; Grönbeck, H. Transient bimodal particle size distributions during Pt sintering on alumina and silica. J. Phys. Chem. C 2015, 119, 989−996.

(19)

The simulated flux is shown in Figure 12 as a solid line with point symbols. For 1 bar, the simulated flux is constant with increasing surface distance and agrees with the value of an open system. For 50 and 100 bar, the flux is smaller than the one given by an open system but larger than the flux according to eq 17. Our calculations show only a small dependence on the pressure. Similar to the parallel surfaces, we expect that the calculated flux will converge to the flux of Fick’s laws for high Δz.

4. SUMMARY AND CONCLUSIONS We have shown that ideal gas properties, in particular, the diffusion constant, can be simulated accurately using the developed kMC model. The model describes the collisions of an explicitly simulated molecule in a background gas. This approach can further be used to simulate mass transport through the gas phase between surfaces. For model systems, we have found that solutions based on Fick’s laws underestimate mass transport if the systems are on the nanometer scale. This can be explained by the long mean free path in the gas phase that leads to a sizable number of trajectories where particles are not scattered by the background gas. Such an effect can not be predicted by standard continuum models, and we expect that this is important in systems where gas-phase diffusion occurs on the nanoscale. Potential applications of the kMC approach are, in particular, simulations involving nanoparticles and processes such as vapor-mediated ripening.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b01816. Details regarding the kMC simulation and derivation of the flux from Fick’s laws (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Elisabeth M. Dietze: 0000-0001-9619-9242 Philipp N. Plessow: 0000-0001-9913-4049 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors acknowledge support by the state of Baden− Württemberg through bwHPC (bwunicluster and JUSTUS, RV bw16G001 and bw17D011). Financial support from the Helmholtz Association is also gratefully acknowledged.

■

REFERENCES

(1) Matera, S.; Maestri, M.; Cuoci, A.; Reuter, K. Predictive-Quality Surface Reaction Chemistry in Real Reactor Models: Integrating FirstPrinciples Kinetic Monte Carlo Simulations into Computational Fluid Dynamics. ACS Catal. 2014, 4, 4081−4092. (2) Matera, S.; Reuter, K. When atomic-scale resolution is not enough: Spatial effects on in situ model catalyst studies. J. Catal. 2012, 295, 261−268. G


Article


(47) Harvie, D. J. E.; Fletcher, D. F. J. Heat Transfer 2001, 123, 486− 491. (48) Kittel, C. Introduction to Solid State Physics, 7th ed.; Wiley: Hoboken, NJ, USA, 1996. (49) Gillespie, D. T. Exact Stochastic Simulation of couple chemical reactions. J. Phys. Chem. 1977, 81, 2340−2361. (50) Carrillo, C.; DeLaRiva, A.; Xiong, H.; Peterson, E. J.; Spilde, M. N.; Kunwar, D.; Goeke, R. S.; Wiebenga, M.; Oh, S. H.; Qi, G.; Challa, S. R.; Datye, A. K. Regenerative Trapping: How Pd Improves the Durability of Pt Diesel Oxidation Catalysts. Appl. Catal., B 2017, 218, 581−590. (51) Johnson, C. A. Generalization of the Gibbs-Thomson equation. Surf. Sci. 1965, 3, 429−444.

(25) Garetto, T. F.; Borgna, A.; Monzón, A. Modelling of sintering kinetics of naphtha-reforming Pt/Al2O3-Cl catalysts. J. Chem. Soc., Faraday Trans. 1996, 92, 2637−2640. (26) Wang, X.; van Bokhoven, J. A.; Palagin, D. Phys. Chem. Chem. Phys. 2017, 19, 30512−30519. (27) Carrillo, C.; Johns, T. R.; Xiong, H.; DeLaRiva, A.; Challa, S. R.; Goeke, R. S.; Artyushkova, K.; Li, W.; Kim, C. H.; Datye, A. K. Trapping of mobile Pt species by PdO nanoparticles under oxidizing conditions. J. Phys. Chem. Lett. 2014, 5, 2089−2093. (28) Hansen, T. W.; DeLaRiva, A. T.; Challa, S. R.; Datye, A. K. Sintering of Catalytic Nanoparticles : Particle Migration or Ostwald Ripening? Acc. Chem. Res. 2013, 46, 1720−1730. (29) Johns, T. R.; Gaudet, J. R.; Peterson, E. J.; Miller, J. T.; Stach, E. A.; Kim, C. H.; Balogh, M. P.; Datye, A. K. Microstructure of bimetallic Pt-Pd catalysts under oxidizing conditions. ChemCatChem 2013, 5, 2636−2645. (30) Johns, T. R.; Goeke, R. S.; Ashbacher, V.; Thüne, P. C.; Niemantsverdriet, J. W.; Kiefer, B.; Kim, C. H.; Balogh, M. P.; Datye, A. K. Relating adatom emission to improved durability of Pt-Pd diesel oxidation catalysts. J. Catal. 2015, 328, 151−164. (31) Jones, J.; Xiong, H.; DeLaRiva, A. T.; Peterson, E. J.; Pham, H.; Challa, S. R.; Qi, G.; Oh, S.; Wiebenga, M. H.; Hernandez, X. I. P.; Wang, Y.; Datye, A. K. Thermally stable single-atom platinum-on-ceria catalysts via atom trapping. Science 2016, 353, 150−154. (32) Xiong, H.; Peterson, E.; Qi, G.; Datye, A. K. Trapping mobile Pt species by PdO in diesel oxidation catalysts: Smaller is better. Catal. Today 2016, 272, 80−86. (33) Plessow, P. N.; Abild-Pedersen, F. Sintering of Pt Nanoparticles via Volatile PtO2: Simulation and Comparison with Experiments. ACS Catal. 2016, 6, 7098−7108. (34) Dadyburjor, D. B.; Marsh, S. P.; Glicksman, M. E. The role of multiparticle-adatom interactions on the sintering of supported metal catalysts. J. Catal. 1986, 99, 358−374. (35) Parker, S. C.; Campbell, C. T. Kinetic model for sintering of supported metal particles with improved size-dependent energetics and applications to Au on TiO2(110). Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 035430. (36) Yang, F.; Chen, M. S.; Goodman, D. W. Sintering of Au Particles Supported on TiO2(110) during CO Oxidation. J. Phys. Chem. C 2009, 113, 254−260. (37) Morgenstern, K.; Rosenfeld, G.; Comsa, G. Local correlation during Ostwald ripening of two-dimensional islands on Ag(111). Surf. Sci. 1999, 441, 289−300. (38) Theis, W.; Bartelt, N. C.; Tromp, R. M. Chemical Potential Maps and Spatial Correlations in 2D-Island Ripening on Si(001). Phys. Rev. Lett. 1995, 75, 3328−3331. (39) Geiser, J.; Blankenburg, S. Monte Carlo Simulations of Elastic Scattering with Applications to DC and High Power Pulsed Magnetron Sputtering for Ti3SiC2. Commun. Comput. Phys. 2012, 11, 1618−1642. (40) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. A new algorithm for Monte Carlo simulation of Ising spin systems. J. Comput. Phys. 1975, 17, 10−18. (41) Kennard, H. Kinetic Theory of Gases; McGraw-Hill Book Company, Inc., 1938. (42) Atkins, P. W.; de Paula, J. Physikalische Chemie, 4th ed.; WILEYVCH Verlag GmbH & Co. KGaA: Weinheim, 2006; p 839. (43) Wedler, G. Lehrbuch der Physikalischen Chemie, 5th ed.; WILEYVCH Verlag GmbH & Co. KGaA: Weinheim, 2007; pp 812−817. (44) Phelps, A. V.; Greene, C. H.; Burke, J. P., Jr Collision cross sections for argon atoms with argon atoms for energies from 0.01 eV to 10keV. J. Phys. B: At., Mol. Opt. Phys. 2000, 33, 2965−2981. (45) Somasundaram, T.; Lynden-Bell, R. M. The velocity distribution of desorbing molecules: a simulation study. Mol. Phys. 1999, 97, 1029−1034. (46) Knox, C. J. H.; Phillips, L. F. Maxwell versus Non-Maxwell Velocity Distributions for Molecules Emitted from a Liquid Surface. J. Phys. Chem. B 1998, 102, 10515−10520. H


Kinetic Monte Carlo Model for Gas Phase Diffusion in Nanoscopic

Recommend Documents