5206
J. Phys. Chem. C 2009, 113, 5206–5211
Synchronization between Kinetic Oscillations Occurring on Neighboring Nanoscaled Supported Platinum Clusters via Thermal Conduction Ching-Cher Sanders Yan,† Yu-Wei Huang,‡ and Shyi-Long Lee*,‡ Institute of Chemistry, Academia Sinica, 128 Academia Road Sec. 2, Nankang Taipei 115 Taiwan, and Department of Chemistry and Biochemistry, National Chung Cheng UniVersity, 168 UniVersity Road, Min-Hsiung Chia-Yi 621, Taiwan ReceiVed: October 20, 2008; ReVised Manuscript ReceiVed: January 22, 2009
Monte Carlo technique is applied to simulate the synchronization between neighboring kinetic oscillators of nanoscaled supported platinum particles. Non-isothermal kinetics caused by local overheating provides an extra thermal conduction channel for synchronization between oscillators. Kinetic oscillation is produced when fast carbon monoxide oxidation is accompanied with slow subsurface oxygen formation and removal. Different effects are investigated on synchronization, including thermal conduction efficiency, heat dissipation efficiency, and distance between the two oscillators. Decreasing thermal conduction efficiency isolates each kinetically oscillating catalyst and results in complicated oscillation patterns. Heat dissipation has the opposite influence on synchronization comparing to thermal conduction. Low heat dissipation from the substrate to the environment maintains the temperature of the system and thus promotes synchronization. In addition, synchronization is destroyed when the distance between two oscillating particles becomes several-fold longer than the platinum particle size. 1. Introduction Kinetic oscillations of heterogeneous catalytic systems have been investigated for over 30 years.1 Major parts of oscillatory kinetics are found in supported catalysts systems with noble metals scattered over SiO2 or Al2O3. The main feature of these systems is that the oscillations periods are irregular and exhibit rich dynamic behaviors, such as quasi-periodicity, mixed mode oscillations, and chaos.2,3 The fundamental reason for such chaotic oscillatory periods is lack of synchronization. The main communicating channels for synchronization between kinetic oscillatory catalysts are mostly gas phase or reactant transportation;1 however, non-isothermal effects caused by local overheating found in 3D array supported catalyst systems also can provide thermal conduction as an extra synchronization channel. In this report, thermal governing factors influencing synchronization efficiency are identified, and their effect on the oscillatory catalysis is demonstrated by employing Monte Carlo simulations. For single crystalline metallic catalysts, the oscillation pattern is expected to be monotonic, based on the prediction of the model assuming a spatially uniform distributed surface. In reality, spatial heterogeneity exists and introduces different reactivity. Thus, different oscillation frequencies can be found in different regions on the same crystalline surface.4 Monotonic oscillation only occurs when efficient synchronization between different regions is achieved. On the other hand, for supported catalysts, assuming that all reactive regimes or spots are contained in a reactor, turnover rates of different catalytic particles oscillating independently would just lead to chaotic kinetic behavior. However, some processes, such as gas-phase diffusion or reactant transportation, provide channels for coupling among different oscillators. Different coupling strengths * Corresponding author.
[email protected]. † Academia Sinica. ‡ National Chung Cheng University.
would cause different kinds of oscillatory behaviors. Strong coupling produces main peaks in which different reaction spots oscillate harmoniously, even though different independently oscillating regimes are found.5 For weak synchronization, Liauw6 et al. proposed a reaction-diffusion model in the meanfield limit, in which gas-phase concentration provides the coupling channel and produces self-similar and mixed-mode oscillations. Therefore, complex catalytic kinetic behaviors result from different coupling strengths among different oscillators. Because of the tendency to smaller scales for catalyst fabrication, the finite size effect of nanoscaled supported catalysts makes the employment of reaction-diffusion equations investigating synchronization problems inappropriate. The main reason is that mean-field approaches cannot describe the interplay or heterogeneity between different faces on a single nanometer particle or on two nearby catalysts on a support. Phase separation of reactants is the special feature often found in the tiny reactors. Monte Carlo simulations conducted by Zhdanov et al.7,8 and Gracia et al.9 have the advantage of taking into account the distributions of reactants more closely to real systems. In addition, they revealed the correlation between reactant surface diffusion and synchronization. Phase separation of reactants indeed profoundly influences the kinetic behavior. However, the understanding of the governing factors influencing synchronization is still far from complete. In particular, it is not fully clear whether oscillatory kinetics on adjacent nanoscaled particles are synchronized or not, and if so, how do they communicate with each other? Besides gas phase and reactant transportation as the common communicating channel, the thermal conduction channel via substrate is often omitted, since overheating of nanoparticles is less possible because of the high surface/volume ratio.10 However, heat accumulates massively during the reaction process11,12 especially for the massive distributed three-dimensional reacting matrix. Because in a three-dimensional array of catalytic particles, rapid and consecutive reactions occur in a local region, the temperature
10.1021/jp809253k CCC: $40.75 2009 American Chemical Society Published on Web 03/12/2009
Kinetic Oscillations on Nanoscaled Pt Clusters
J. Phys. Chem. C, Vol. 113, No. 13, 2009 5207
difference among particles is not obvious; thus, the local region can be overheated. Therefore, a detailed investigation of the heat conduction channel for kinetic oscillation synchronization is necessary. In general, supported catalytic particles are dispersed in a three-dimensional matrix; kinetic oscillation is synchronized among a group of different particles. However, detailed investigation on a pair of kinetic oscillatory particles can help to demonstrate more clearly the different governing factors influencing synchronization. Increases in the distances between particles would decrease the heat exchange efficiency more rapidly among reaction spots in real systems than the paired system considered here. However, the parameters for heat exchange efficiency can be easily adjusted to mimic more realistic situations and do not change the main kinetic behavior of the system. The kinetic oscillations are produced by slow subsurface oxygen formation from adsorbed oxygen and removed by carbon monoxide oxidation accompanied by rapid carbon monoxide oxidation with adsorbed oxygen. The temperature variation accompanied by kinetic oscillation originates from the overheating of reacting particles. Heat redistribution among catalysts and substrate forms a channel for synchronization between different oscillators via thermal conduction. The report is organized as follows. Section 2 gives details of the simulation, including reactions, diffusion, and heat redistribution. In section 3, simulation results and discussion are given. The conclusion is given in the last section. 2. Methodology 2.1. Mechanism. The subsurface oxygen formation and removal model was first proposed by Sales13 et al. to explain the oscillatory kinetics of carbon monoxide oxidation over platinum particles. The mechanism contains several steps as follows:
COgas T COads
(1)
O2,gas f 2Oads
(2)
COads + Oads f CO2,gas
(3)
Oads f Osub
(4)
COads + Osub f CO2,gas
(5)
where Osub stands for subsurface oxygen, which forms a stronger bond with metal surfaces and becomes less reactive with COads. COads reacts with Oads via the Langmuir–Hinshelwood (LH) reaction of eq 3, which is a rapid step compared to others. Thus, the oxidation rate decreases when the Osub coverage increases. Osub can only be moderately removed by COads during the low rate state and provides new vacant adsorption sites for newly adsorbed oxygen. The reaction rate will rise again when the Oads coverage is high. Therefore, oscillatory behavior is observed when the system is transferring back and forth between the low and high rate states. 2.2. Monte Carlo Simulation of Reactions and Diffusion. In this study, reaction kinetics of the subsurface oxygen model was simulated over a squared lattice with open boundary conditions to mimic the exposed surface of a Pt pellet deposited on the support. Note that the periodic boundary condition is
not adopted here, because the open boundary condition is more appropriate to imitate the supported Pt(100) nanoscaled catalysts.14 The size of the lattice is selected as 30 by 30, which is comparable to the typical platinum pellet size. In addition, the selected lattice size is big enough to ensure that a single catalyst produces regular oscillations, not affected by a particle’s finite size effect15 producing irregular oscillation. When chaotic oscillations are found from a couple of catalysts, it is solely caused by the non-synchronization between them, but not the finite size effect. To introduce probabilities for each step for MC simulation, the sum of the rate constants of the LH step and CO desorption are selected and regarded as a normalization constant, because (a) these two steps are the next movement of an adsorbed CO in the simulation and (b) the LH step is fastest among all the elementary reaction steps. Probabilities for other steps are then taken from the ratio in proportion to the sum of the two rate constants.
pi ≡
ki kLH + kdes
(6)
In addition, one Monte Carlo step (MCS) tMC is related to a fixed real time interval ∆t as eq 7
tMC ) (kLH + kdes)∆t
(7)
Fast diffusion of COads step is running separately from the slow reaction steps by introducing a dimensionless parameter Ndiff to characterize the trial ratio between diffusion and reaction steps. Thus, the simulation time step is related only to the reaction steps but not to diffusion. Ndiff is estimated to be ∼400 for the diffusion coefficient of CO over the Pt surface (Df ) 1.2 × 10-7 cm2/s). Different Ndiff values, except zero, are found to cause only a slight difference in the kinetic behavior, and below Ndiff were taken to be 100. The CO diffusion trial is conducted by randomly selecting a site from the lattice. When a randomly selected neighboring site of the selected COads is empty, the diffusion process is completed by moving the COads to the vacant site. COads diffusing from the catalyst border to the substrate is not allowed because the binding energy of CO to Pt (38 kcal/mol) is higher than that ( Frea ≡ 1/(1 + Ndiff), diffusion is conducted. Otherwise, a reaction trial is executed. A new random number F2 is generated. A different step is conducted depending on the occupation of (xs, ys): (a) When (xs, ys) is empty, adsorption is possible. CO is adsorbed when F2 < pCO. For the case of pCO < F2 < pCO + pO2, O2 dissociative adsorption is performed when randomly selected next nearest neighboring site (xnnn, ynnn) is also empty for the oxygen repulsive interaction effect.14 (b) When (xs, ys) is occupied by Oads, it would convert to Osub for F2 < pox. (c) When (xs, ys) is occupied by CO, desorption would proceed for F2 < pdes. Otherwise, rapid LH reaction is accomplished and produces two empty sites when randomly selected (xnn, ynn) is occupied by Oads. If it is Osub at (xnn, ynn), a third random number F3 is generated to check whether removing Osub is accomplishable or not. It would become a successful trial for F3 < pred. By
5208 J. Phys. Chem. C, Vol. 113, No. 13, 2009
Yan et al.
TABLE 1: Thermal Parametersa density heat capacity Avogadro’s number closest Pt-Pt lattice separation heat of CO adsorption heat of O2 dissociative adsorption heat of CO desorption heat of CO oxidation
F Cp NA lc
21.45 g/cm3 3.107 × 10-5 kcal/(g · K) 6.0221367 × 1023 mol-1 3.9242 × 10-8 cm
∆HCO,ads ∆HO2,ads
-38 kcal/mol -54 kcal/mol
∆HCO,des ∆HLH
38 kcal/mol 24.37 kcal/mol
a From Cisternas, et al.17 Note: there is no heat released or absorbed accompanying the formation and deformation of subsurface oxygen.
TABLE 2: Kinetic Parametersa CO adsorption rate constant CO partial pressure O2 dissociative adsorption rate constant O2 partial pressure desorption of CO CO oxidation subsurface oxygen formationc subsurface oxygen back conversionc a
kCO pCO kO2
3.135 × 105 1/s · mbar 3.414 × 10-5 mbarb 2.343 × 105 1/s · mbar
pO2 0 kdes Edes 0 kLH ELH k0sub Esub 0 krev
4.567 × 10-5 mbarb 2 × 1016 1/s 38 kcal/mol 1.704 × 107 1/sb 10 kcal/mol 8.971 × 103 1/sb 8.2 kcal/mol 2.984 × 103 1/sb
Erev
8.2 kcal/mol
17
From Cisternas et al. unless otherwise indicated. b The partial pressure and prefactors is deduced from the probability sets from Zhdanov’s work.14 c From Oertzen et al.18
dTc ) dt
∑ (∆Ti · rate of step i) - Kcd(Tc - T0)
where reaction rate and heat dissipation constant Kcd are for the whole metal particle because of good heat conduction for a small metal pellet. T0 is the temperature of the environment. In the simulation, heat accumulates when each reaction step is executed successfully. Temperature decrease due to heat dissipation is calculated in each MCS. For simplicity, including all possible effective factors for heat dissipation, constant Kcd (1/s) is taken as a governing parameter. According to our experience, characteristic time scale (1/Kcd) has to been 10 s or longer to observe temperature variation, which is an extreme case, comparing to the typical value of about 10-12 s. However, a local overheating situation is still possible mainly due to the finite temperature difference between neighboring envioronments.10 Varying temperature influences the kinetics of the catalysis exponentially through the Arrhenius form. Probabilities for MC simulation are recalculated according to eq 6 using the values of prefactor and activation energy given in Table 2. Conventionally, reaction trials of the same number as site number on the surface are executed in one MCS when the simulation is performed over a 30 × 30 lattice at constant temperature. Under non-isothermal conditions, rate constants change with temperature variation. Thus, the number of reaction trials for each MCS changes according to the ratio of the sum of two rate constants at temperature Tc and at Tref.
Trial number in a MCs at Tc ) L × L ·
Figure 1. Construction of simulation lattices and definitions of symbols.
the way, pred is obtained by taking kred normalized to kLH only. The LH reaction is hindered due to less reactive Osub formation. 2.3. Temperature Variation. All the MC simulations were started with a clean surface. A non-isothermal effect is caused by self-heating, so that temperature variation ∆Ti accompanied by the ith step was estimated according to eq 8 with all the thermal parameters listed in Table 1.
∆Ti ) -
∆Hi F · Cp · NA · (lc · L)3
(8)
where the length of squared lattice L ) 30 in this report. The Pt particle size is assumed within a reasonable range without quantum effect. The geometry of the particle is assumed to be a cube. Heat dissipation to the environment would cause the temperature of catalyst Tc decrease. The governing equation used in our simulation was
(9)
i
kLH(Tc) + kdes(Tc) kLH(Tref) + kdes(Tref) (10)
where Tref can be selected arbitrarily. It was selected to be 517 K in our simulation. By doing this, the different kinetic pace at different temperature can be discriminated in the MC simulation. 2.4. Synchronization via Thermal Conduction. The simulation of reaction and diffusion processes is started independently on two clean squared lattices. At the same time, to describe thermal conduction between them, another one-dimensional lattice Ts (K) is adopted to record the temperature profile of the substrate between two catalysts (see Figure 1). It is divided into d unit, by which each unit length is equal to the linear size of one typical Pt pellet, representing the distance between the two oscillators. The temperature is initially set the same as the temperature of the environment T0 (K). Adsorptions, diffusion, or reactions are not allowed over the substrate but only on the catalysts, so that thermal conduction is the only channel coupling the two oscillatory catalysts. After diffusions and reaction trials are finished within one Monte Carlo step (MCS) for each catalyst, the temperature changes in one MCS of the two catalysts Tc (K) are calculated by eq 11 including the temperature raised by the heat of reaction Tr (K), heat dissipation to the environment, and thermal conduction to the substrate.
Kcd (MCS-1) and Kcs (MCS-1) stand for the coefficient of heat dissipation from the catalyst and thermal conduction to
Kinetic Oscillations on Nanoscaled Pt Clusters
J. Phys. Chem. C, Vol. 113, No. 13, 2009 5209
Figure 2. Temperature, reaction rate (CO2 molecules per site per MCS), and average reaction rate of two reactive 30 × 30 lattices calculated for Kcd ) 0.05, Ksd ) 0.03, Kcs ) 5, Kss ) 50, and d ) 5.
Figure 3. Temperature, reaction rate (CO2 molecules per site per MCS), and average reaction rate of two reactive 30 × 30 lattices calculated for Kcd ) 0.05, Ksd ) 0.03, Kcs ) 0.05, Kss ) 0.5, and d ) 5.
the substrate, respectively. Temperature changes in one MCS of the substrate are calculated by eq 12.
Ksd (MCS-1) is the coefficient of heat dissipation from the substrate. Note that, considering local overheating effect here,
5210 J. Phys. Chem. C, Vol. 113, No. 13, 2009
Yan et al.
adopted Kcd and Ksd are much smaller than in reality. After the heat redistribution, the diffusion and reaction processes are conducted again independently over two catalyst lattices with recalculated probabilities, since the temperature for each catalyst then is different from the last MCS. In order to characterize the extent of deviation from synchronization, the discrepancy is estimated by δ, which is defined
(∑ n
δ)
)
1 t |t - tBi |/ × 100% n i)1 Ai 2
(13)
tAi and tBi are the simulation times finding the ith peak of catalysts A and B. jt is the averaged oscillation period of the two catalysts. Only when the simulation temperature is appropriate and obvious oscillation is found with high and low rates are the periods of two oscillators very similar under similar conditions. When two oscillators have similar periods, δ defines
Figure 6. Deviation from in-phase synchronization, δ, as a function of heat dissipation efficiency of catalyst Kcd and substrate Ksd with Kcd/ Ksd ) 5:3, Kcs ) 1, Kss ) 10, and d ) 5.
the phase difference. δ ) 0 refers to good synchronization between two oscillators, namely, in-phase. Statistically, δ > 30% is regarded as poor coupling and without a good communicating channel. However, oscillation periods are not commonly regular, because distribution heterogeneity plays an important role for the nanosized system. It makes the mean-field approximation not appropriate for this system. Thus, the Monte Carlo algorithm is employed to take into account the random heat rising from two different oscillatory sources on the thermal conduction synchronization. The simulation time is long enough for 106 oscillation periods to make the statistics of deviation from synchronization meaningful. 3. Results and Discussion
Figure 4. Deviation from in-phase synchronization, δ, as a function of thermal conduction efficiency between catalyst and substrate, Kcs, and between substrate and substrate, Kss, with Kcs/Kss ) 1:10, Kcd ) 0.05, Ksd ) 0.03, and d ) 5.
Figure 5. Deviation from in-phase synchronization, δ, as a function of distance between two catalysts with Kcd ) 0.05, Ksd ) 0.03, Kcs ) 1, Kss ) 10.
For all the simulation work, the temperature for simulation T0 is selected as 517 K, as the typical temperature for supported catalyst finding kinetic oscillation.16 The catalyst heat dissipation coefficient Kcd is much smaller than in reality under a local overheating situation; thus, the temperature of catalysts rises a few degrees corresponding to the high rate as usual experimental observations, as shown in Figure 2. The simulation was first performed to examine the effect of different thermal conduction. Different conduction efficiency represents different materials used as supportive substrates. The ratio of Kcs and Kss was fixed to 1:10 to represent the difference of thermal conductivities of the interface between catalyst and substrate and the internal zone of the substrate. When the thermal conduction efficiency is good, the average reaction rate shown in Figure 2 is quite similar to the rate calculated from each catalyst, respectively. The system’s temperature also oscillates accordingly with reaction rate, shown in Figure 2. Only a few peaks with a small split are found. The estimated deviation δ from in-phase synchronization is smaller than 20%. The simulation in Figure 2 is regarded as a synchronized couple. When the thermal conduction efficiency is decreased to 0.01, as shown in Figure 3, the averaged high rate peaks become less synchronized. Larger peak splits are found due to the lack of coupling between two oscillatory reactors. Averaged high rate peaks are wider and shorter, because a rising rate only comes from one catalyst, but the value taken was the average over two catalysts. In Figure 4, the dependence of δ, the deviation from synchronization, on thermal conduction is shown. With
Kinetic Oscillations on Nanoscaled Pt Clusters increasing thermal conduction efficiency, synchronization becomes better and δ decreases from 50% to 25%. Thus, good thermal conduction provides the channel coupling different oscillators and results in more periodic and higher overall rate. In addition, when comparing the number of oscillation cycles from individual catalysts, both reactors A and B in Figure 3 completed one more cycle than those in Figure 2. Because of high thermal conduction, the substrate distributes heat more easily and makes the temperature of the catalysts lower. Lower catalyst temperature corresponds to longer oscillation periods. This simulation result is consistent with experimental observations. Further simulations of the effect of distance on synchronization are given in Figure 5. Increasing the distance between two oscillators diminishes the synchronization between them. A deviation up to 50% is regarded as two independent oscillators when distance d > 15. The changes of δ when d > 15 only represent the period fluctuation of two individual oscillators in statistics. Different environmental conditions of the substrate would change its own heat dissipation and further influence the coupling strength of two oscillating reactors. Figure 6 represents the simulation results of the extent of synchronization with increasing heat dissipation efficiency of the substrate. Kcd and Ksd are increased together at the ratio 5:3 with Kcs and Kss fixed in order to simulate the difference in heat dissipation efficiency between Pt and the support surface, and the distance is equal to 5. The reaction heat is maintained in the system and increases the potential of each catalyst to influence the other. Thus, synchronization is promoted. With increasing heat dissipation, the potential no longer exists, which causes increasing deviation from in-phase synchronization. 4. Conclusion The extent of synchronization is investigated between two oscillatory nanometer-sized supported catalysts, on which CO oxidation take place with temperature variation, by conducting Monte Carlo simulation. Thermal conduction via the support
J. Phys. Chem. C, Vol. 113, No. 13, 2009 5211 plays an important role in the coupling of two nearby oscillators. With heat transfer, neighboring oscillators build up communication channels, synchronization of oscillation is promoted, and simple oscillation cycles are produced, without split high rate peaks. An index, deviation from synchronization, δ, is designed to shed some insight into the strength of the coupling channel among oscillating reactors. Increasing distance or heat dissipation weakens the coupling strength and destroys the synchronization. Increasing thermal conduction efficiency opens the coupling channel, and thus synchronization is obtained. Acknowledgment. Financial assistance from National Science Council, Taiwan, is gratefully acknowledged. References and Notes (1) Schu¨th, F.; Henry, B. E.; Schmidt, L. D. AdV. Catal. 1993, 39, 51. (2) Lindstrom, T. H.; Tsotsis, T. T. Surf. Sci. 1986, 171, 349. (3) Hudson, J. L.; Hart, M.; Marinko, D. J. Chem. Phys. 1979, 71, 1601. (4) Schmitz, R. A.; D’Netto, G. A.; Razon, L. F.; Brown, J. R. In Chemical Instabilities; Nicolas G., Baras, F., Eds.; Springer: Boston, 1984. (5) Schu¨th, F.; Wiche, E. Ber. Bunsen-Ges. Phys. Chem. 1989, 93, 491. (6) Liauw, M. A.; Plath, P. J.; Jaeger, N. I. J. Chem. Phys. 1996, 104, 6375. (7) Zhdanov, V. P.; Kasemo, B. Chaos 2001, 11, 335. (8) Zhdanov, V. P.; Kasemo, B. Surf. Sci. 2002, 511, 23. (9) Gracia, F. J.; Wolf, E. E. Chem. Eng. Sci. 2004, 59, 4723. (10) Zhdanov, V. P.; Kasemo, B. Catal. Lett. 2001, 75, 61. (11) Yan, C. C. S.; Chuang, W. T.; Chaudhari, A.; Lee, S. L. Chem. Phys. Lett. 2004, 400, 245. (12) Yan, C. C. S.; Chuang, W. T.; Chaudhari, A.; Lee, S. L. Appl. Surf. Sci. 2005, 252, 784. (13) Sales, B. C.; Turner, J. E.; Maple, M. B. Surf. Sci. 1982, 114, 381. (14) Zhdanov, V. P. Catal. Lett. 2000, 69, 21. (15) Zhdanov, V. P. Catal. Lett. 2004, 93, 135. (16) Lauterbach, J.; Bonilla, G.; Pletcher, T. D. Chem. Eng. Sci. 1999, 54, 4501. (17) Cisternas, J.; Holmes, P.; Kevrekidis, I. G.; Li, X. J. Chem. Phys. 2003, 118, 3312. (18) Oertzen, A.; Mikhailov, A. S.; Rotermund, H. H.; Ertl, G. J. Phys. Chem. 1998, 102, 4496.
JP809253K
