Towards a first-principles chemical engineering:
Transport limitations and bistability in in situ CO oxidation at RuO(110)
We present a first-principles based multiscale modeling approach to heterogeneous catalysis that integrates first-principles kinetic Monte Carlo simulations of the surface reaction chemistry into a fluid dynamical treatment of the macro-scale flow structures in the reactor. The approach is applied to a stagnation flow field in front of a single-crystal model catalyst, using the CO oxidation at RuO(110) as representative example. Our simulations show how heat and mass transfer effects can readily mask the intrinsic reactivity at gas-phase conditions typical for modern in situ experiments. For a range of gas-phase conditions we furthermore obtain multiple steady-states that arise solely from the coupling of gas-phase transport and surface kinetics. This additional complexity needs to be accounted for when aiming to use dedicated in situ experiments to establish an atomic-scale understanding of the function of heterogeneous catalysts at technologically relevant gas-phase conditions.
First-principles kinetic Monte Carlo (1p-kMC) simulations have evolved into an important tool in the modeling of heterogeneous catalytic processes Reuter (2010). The success of the approach relies on the accurate treatment of two central aspects for the reactive surface chemistry: A first-principles description of the involved elementary processes and an evaluation of their statistical interplay that fully accounts for the correlations, fluctuations and spatial distributions of the chemicals at the catalyst surface. Particularly if suitably combined with sensitivity analyses Meskine et al. (2009), 1p-kMC simulations thus offer the prospect of an error-controlled and quantitative microkinetic modeling of the surface catalytic function. In particular for technologically relevant environments, i.e. at near-ambient reactant pressures and elevated temperatures with concomitant higher product formation rates, a third aspect comes into play that needs to be accounted for to properly describe the observable catalytic conversions. This is the flow of heat and mass in the given reactor geometry, for instance the transport of formed products away from the active surface and how efficiently the large amount of heat generated by the exothermic surface reactions can dissipate into the system.
Corresponding macro-scale flow structures are suitably described at the continuum-level by the transient Navier-Stokes equations together with energy and species governing equations. The methodological objective of the present work is then to couple 1p-kMC into such a fluid dynamical (FD) framework, thereby augmenting the accurate treatment of the reactive surface chemistry provided by the prior technique with the capability to account for heat and mass transport effects. With a brief account of the main results already given in ref. Matera and Reuter, 2009, we focus here in particular on a detailed description of this methodology. While the presented approach can readily be coupled with any computational fluid dynamics (CFD) software enabling the treatment of arbitrary reactor geometries, we develop it in the following specifically for a stagnation flow field in front of a flat-faced model catalyst. We argue that this is a suitable, though admittedly simplified reactor geometry to qualify transport effects in modern in situ studies aiming at an atomic-scale understanding of the catalytic function of single crystals in technologically relevant environments Stierle and Molenbroek (2007).
The focus of such studies lies on possible differences in the surface chemistry compared to operation in ultra-high vacuum (UHV), where the function of model catalysts has been extensively studied in the past. In order to discern corresponding so-called ”pressure gap” effects, it is important to assess if heat and mass transfer effects noticeably mask the true intrinsic reactivity in the in situ conditions. Particularly for the often studied CO oxidation reaction at late transition metal catalysts there are good reasons to suspect that transport limitations might not be negligible. Typically reported activities indicate a high rate of mass conversion at the surface concomitant with a large heat release. Using the established 1p-kMC model for the CO oxidation at RuO(110) Reuter et al. (2004); Reuter and Scheffler (2006) as a representative example we illustrate with the coupled 1p-kMC–FD approach that the peculiarities of the single-crystal reactor geometry lead indeed readily to heat dissipation and mass transport limitations that severely affect the observable catalytic function. Key factors are the degree of heat conduction at the backside of the thin single-crystal, and the propensity to build-up a product boundary layer above the flat-faced surface. Obivously, such reactor-dependent effects need to be disentangled, understood and controlled when aiming to compare data obtained by different experimental setups, and when aiming to conclude on the actual surface chemistry at technologically relevant gas-phase conditions.
Ii Macro-Scale Flow Structures: Continuum Fluid Dynamics
We develop our approach for a simple reactor geometry suitable to discuss heat and mass transfer effects at a flat-faced model catalyst. In the so-called stagnation flow geometry Kee et al. (2003) shown in Fig. 1 the gas mixture enters through an inlet at a macroscopic distance away from the surface. At this inlet the flow is directed towards the catalyst surface, with no variation of gas composition, velocity, and temperature in the direction perpendicular to the flow. As schematically indicated in Fig. 1 the advantage of such a geometry is that it results in an axisymmetric flow profile. Neglecting edge effects, i.e. the finite lateral extension of the model catalyst, this flow profile can effectively be described by a one-dimensional boundary-value problem. As further detailed below, this eases the analysis of the influence of the reactor setup significantly and allows to extract the relevant physics without being riddled by algebra and numerics.
However, it is not only this symmetry imposed simplification which makes the stagnation flow setup appealing or better its realization desirable. In the spirit of the Surface Science approach to heterogeneous catalysis, the fundamental objective of in situ studies of single crystals is to obtain insight into their intrinsic reactivity in near-ambient environments, for example as function of temperature and reactant partial pressures . For this to be well-defined, a central prerequisite is that all, or at least a dominant fraction of the active surface sees the same gas phase. If one considers e.g. a flow geometry where the stream of reactants would approach the surface from the side – thus in some sense an opposite scenario to the here discussed perpendicular stagnation flow – then this is clearly not the case. Due to the on-going conversion of reactants into products the gas-phase composition close to the surface would gradually change across the lateral extension of the single crystal. If there are non-negligible heat transfer limitations, this goes hand in hand with a non-uniform temperature profile parallel to the surface. Under such conditions, making a defined assignment of observed turnovers to specific pressures and temperatures becomes essentially intractable. In contrast, in the stagnation flow geometry at least the entire center part of the active surface sees the same gas phase, thereby facilitating the analysis.
In the stagnation flow setup shown in Fig. 1 the most relevant spatial coordinate is thus the direction perpendicular to the catalyst surface. In the axisymmetric problem we denote the other, radial coordinate with . Using the inlet height as zero reference for , the catalyst surface is then at and with a thickness of the single crystal its backside is located at . In a continuum mechanical description the system is correspondingly characterized by two spatial regions, the flow chamber () and the sample (), as well as three important interfaces, the inlet (), the surface () and the catalyst backside ().
In the following subsections we will successively describe the modeling carried out for each of these regions and interfaces. In this manuscript this modeling aims at a description of steady-state operation, in which chemicals get converted at the active surface at a stable, time-independent rate. This rate is the turnover frequency (TOF) in units of molecules per surface area and time. We note that this time-independent formulation is a convenience, not a necessity. In fact, in particular the coupling scheme integrating the surface chemistry into the FD environment is also applicable to transient problems and we will briefly mention routes in this direction throughout the manuscript. In the same spirit it is clear that the approach is, of course, not restricted to the simple CO oxidation reaction, on which we will concentrate from now on for clarity.
ii.1 Interface I: Inlet
At the inlet the gas flow is fully controlled. For the FD description this defines boundary values for the temperature , the partial pressures (with O, CO, CO), the total pressure , and the axial velocity . We set the radial velocity to as is the case in stagnation flow reactors with a sieve-like showerhead as inlet Kee et al. (2003), and we restrict our attention to flow situations with no circumferential motion. Instead of the partial pressures it is more common to use the total density and mass fractions as independent fields. For the present catalytic context the conversion is readily achieved through the ideal gas law
where is the mass of species and is the Boltzmann constant. Similarly, we will see below that it is more useful to consider the scaled radial velocity in the stagnation flow equations.
ii.2 Region I: Flow chamber
ii.2.1 Flow equations
The continuum mechanical description of the heat and mass transport in the flow chamber is based on balance equations for total mass, species mass, momentum and internal energy. For the present purposes, the general formulation is much simplified by treating the gas phase as a Newtonian fluid with vanishing bulk viscosity, by the absence of relevant gravitational forces and by the absence of significant gas-phase chemical conversions in the context of low temperature CO oxidation (i.e. the reactions are confined to the catalyst surface). Another significant simplification arises for the laminar flows of interest here in that we can work in the low Mach number approximation: For the corresponding small flow velocities the variations of the total pressure over the whole flow chamber domain will be much smaller than its absolute magnitude. Gradients of this small variation, denoted as dynamic pressure , can then be neglected in all equations except the momentum balance.
Working in the low Mach number approximation the specific equations for the steady-state stagnation flow problem are e.g. discussed in detail by Kee, Coltrin and GlarborgKee et al. (2003). The key assumption in the derivation is that variations of partial pressure, temperature and axial velocity in radial direction are much smaller than in axial direction, at least near the symmetry axis in the center of the catalyst. A lowest order expansion in the radial dependence leads then to a set of differential equations in which all fields depend only on the axial coordinate .
Here, is the so-called radial pressure curvature, the shear viscosity, the thermal conductivity, the specific heat capacity, the specific heat capacity of species , and the axial component of the diffusive mass flux of species .
Together with the ideal gas law, eq. (1), and in the low-Mach-number-limit the condition
this set of equations allows to solve for all dependent fields, , , , , , and . In fact, the axial momentum balance, eq. (5), is decoupled from the other equations. All fields except the dynamic pressure can be determined without it, and it can therefore be used afterwards to fix from the other calculated fields. The problem gets fully determined by boundary conditions at the inlet and at the surface. Specifically, the second order equations demand independent information about , and at both ends of the domain, and the first-order continuity equation requires information about on one boundary. As this information about is provided at both boundaries in the here discussed finite-gap stagnation flow, i.e. at inlet and surface, the resulting overdetermination is resolved by treating the unknown as an eigenvalue of the problem, i.e. its magnitude is adjusted to satisfy the additional boundary condition.
ii.2.2 Thermophysical and transport parameters
|(meV)||196||269||291, 167, 83 (2x)||[Linstrom and Mallard, 2009]|
|(eV)||0||-1.179||-4.074||[Stull and Prophet, 1971]|
What remains for the numerical solution are values for the transport coefficients and , and a diffusion theory relating the diffusive mass flux to composition gradients. Further, we need expressions for the isobaric specific heat capacities and , as well as for the specific enthalpy (see section IIC below). For the gas phase considered here the mixture specific heat is simply the mass-weighted sum of the species specific heats
For the themselves we consider translational, rotational and vibrational degrees of freedom, treating the prior two in the classical limit and the latter in the harmonic approximationBird et al. (1960)
with , and and the number of vibrational and rotational degrees of freedom, respectively. For the three linear molecules O, CO, and CO , and Table 1 lists their vibrational frequencies as taken from ref. Linstrom and Mallard, 2009.
Computing the transport coefficients for a multi-component gas phase in complete generality is a very complex task of its own. For the accuracy level and gas-phase conditions of interest for our study effective semi-empirical molecular transport models provide fortunately a fully sufficient description. From the manifold of such existing models we adopt a strategy described by Cloutman that is also used in the reactive flow CFD program COYOTECloutman (2000). In a first step, this strategy relies on standard mixture-averaged approaches that relate and of the multi-component gas mixture to the pure species values. In the case of the viscosity, this is the Wilke formula Cloutman (2000); Bird et al. (1960)
where are the pure species conductivities. Note that is exactly the same as in Wilkes formula, i.e. it depends on the viscosities and not on the thermal conductivities. From the kinetic theory of dilute gases the expression Cloutman (2000); Bird et al. (1960)
provides the species viscosities in terms of two empirical parameters, a characteristic diameter and a characteristic energy . For a gas with Lennard-Jones interaction between the molecules and are the two parameters defining the interaction potential Kee et al. (2003). For the general case referencing to a Lennard-Jones model just provides a convenient way of representing the temperature dependence of the transport coefficients and the two parameters need not to have a microscopic meaning. Values for these parameters for a wide range of species are found in data bases and we list in Table 1 those for the species O, CO, CO needed hereCloutman (2000). From the thus determined viscosities the thermal conductivities are derived via the Eucken correctionCloutman (2000); Bird et al. (1960)
In the present application the remaining diffusive mass fluxes are predominantly driven by the concentration gradients and can then be implicitly obtained from corresponding Stefan-Maxwell equations Kee et al. (2003); Bird et al. (1960) at each point in the flow field
with , and
which also needs only the empirical characteristic diameter and energy of each species, cf. Table 1.
ii.3 Interface II: Surface
As already mentioned, we need to specify a number of boundary conditions at the solid surface to fully determine the set of stagnation flow equations. For the radial fluid velocity this is the standard no-slip condition, i.e. . The normal component of the velocity at the surface ( in our case) is commonly termed ”Stefan velocity” Kee et al. (2003); Deutschmann (2008). It is calculated by considering the mass balance at the surfaceMüller (1985), and a finite Stefan velocity accounts for transient storage or release of species due to a changing average surface composition. For the here discussed steady-state operation this is not the case and we have .
As the CO oxidation reactions at the surface are the only processes that consume reactants and yield products in stationary operation, the boundary condition for the diffusive mass flux is
where the chemical source term for each species is simply given by the overall rate of reaction events (per area and time) times the stoichiometry coefficient and mass of the species in the reaction. For the simple CO oxidation reaction, , , and .
The heat release connected to the exothermic conversions must also enter the heat balance at the surface. With the previously introduced boundary conditions the surface energy balance reduces to the mere requirement that the normal heat flux is continuous at the surface
Here, the two terms on the left hand side account for the heat transported away by the gas phase and the heat released by the surface chemical reactions, respectively, which must balance the heat flux into the solid . The temperature-dependent contribution to the specific enthalpies is determined completely equivalently to the specific heat capacities, i.e. by considering translational, rotational and vibrational degrees of freedom as described above. The additionally required zero Kelvin component including the vibrational zero point energies can be drawn from thermochemical tables Stull and Prophet (1971), and for completeness they are included in Table 1 for the species O, CO, and CO.
At first glance this use of experimental thermochemical data might seem inconsistent with the use of first-principles energetics in the 1p-kMC simulations. Our choice is motivated by the consideration that all other transport parameters are equally derived from experiment. The empirical transport models provide a fully sufficient and controlled description of the macro-scale flow structures at the accuracy level of interest for this study. As argued in more detail below such a description is reached for the reactive surface chemistry through the use of first-principles based 1p-kMC modeling. We are thus faced with the problem of matching these two descriptions. In our approach this match occurs uniquely through the TOFs. Whatever energetics is required for their determination comes exclusively and consistently from first-principles, here density-functional theory with a semi-local exchange-correlation functional. Vice versa, the entire transport description is exclusively and consistently based on experimental numbers, such that the balance leading to its effective accuracy is not disturbed by occasional parameters as e.g. the coming from approximate first-principles theory. While we feel that such considerations are necessary for the envisioned error-controlled multiscale modeling, we note that in practice the semi-local functional employed in the 1p-kMC overestimates the zero Kelvin enthalpy change for the CO oxidation reaction in Table 1 by only about 10%, such that none of the conclusions reported below would be touched when using this energetics instead of the experimental one.
Since we are considering heat transport in both the solid and the gas phase, we finally need two conditions describing the change of the temperature field when crossing the surface. Equation (21) can only serve as boundary condition for the gas-phase heat transport. Another one is needed for the heat transport in the solid. Here, we want to assume that the temperature is continuous across the gas-solid interface. Within our interest in modeling the reactivity at near-ambient conditions, this should be rather well ensured by the frequent gas-surface collisions.
ii.4 Region II and Interface III: Sample and sample backside
Neglecting possible heat-induced deformations or phase transformations we only account for heat conduction through the single crystal. For the stagnation flow problem there is no radial variation of the gas-phase temperature. We maintain this description within the sample and correspondingly also model the heat transport as a one-dimensional problem described by Fourier’s law
where is the heat conductivity of the sample, which for simplicity we assume to be temperature-independent. Solution of this equation requires fixing a boundary value at the backside of the single crystal , which thus describes the degree of heat dissipation that is possible e.g. through radiative loss or the contact of the single crystal with the sample holder. Specifying this for real experimental setups is a demanding task and we suspect that in present in situ experiments this value will vary largely. Addressing these specificities in a quantitative way is clearly outside the capabilities of the present idealized reactor model. With real experimental setups lying anywhere in between we therefore analyze the relevance of this factor for thin single crystals by considering two opposite extremes: First a fixed temperature at the sample backside to mimic a highly efficient heat coupling of the crystal to the system, and second a zero heat flux at the sample backside to represent a well insulated sample. For either boundary condition eq. (22) can be solved analytically and provides in turn the missing boundary condition for the surface heat balance, eq. (21).
|Isolated sample (adiabatic limit)|
|Connected sample (isothermal limit)|
where is the temperature at the back of the crystal. While in principle this temperature could have any value (e.g. through controlled heating) we will assume here that it is identical to the one of the outside system and thus to the temperature of the incoming gas at the inlet, i.e. . The remaining material parameter to specify in the resulting boundary condition is the bulk heat conductivity. An experimental value for bulk RuO is = 0.50 W cm K.Ferizovic et al. (2009) However, almost all in situ work on RuO(110) has in fact been performed on ultra-thin films grown on Ru crystals, which would suggest that the value = 1.17 W cm K is more appropriate Lide (1995). Fortunately, for the results reported below it makes no difference which value for this quantity is used. At the rather high thermal conductivity of either metallic Ru or the metallic oxide RuO this dissipation channel is so dominant, that – as soon as it is enabled in the surface heat balance, eq. (21) – it simply ensures that the surface temperature remains at the nominal value . Regardless of the specific value of the second boundary condition is therefore simply equivalent to modeling an isothermal reactor limit, while the first boundary condition would correspond to the adiabatic limit.
ii.5 Numerical solution
where are the so-called differential components and the algebraic component. In our case the differential components are density , mass fractions , temperature , velocity components and , diffusive mass fluxes , intrinsic heat flux , pressure curvature , and . The algebraic component is the gradient of the density, which is determined by the requirement that the total pressure is constant between inlet and surface.
We solve the DAE boundary value problem numerically using the COLDAE package Ascher and Spiteri (1994). This package uses piecewise orthogonal collocation at Gaussian points and has an adaptive mesh strategy allowing for an error-controlled solution. We use the default option controlling the number of intermediate points and an initial equidistant mesh with 10 spacings. In all simulations we employ a tolerance of for each differential component. In order to improve the stability and have full error control, all variables, dependent or independent, are presented in appropriate units. The employed units are the inlet-surface distance for length, and for the velocity, density, and mass fractions their values at the inlet. The employed temperature scale is Kelvins. We use the representation for the temperature, so that this renormalized temperature is always zero at the inlet. The mass fluxes, heat flux, and density gradient are expressed in multiples of , , and , respectively, where is the mixture averaged diffusion coefficientKee et al. (2003). Finally, the radial pressure curvature and are scaled with and , respectively.
The software uses a (damped) Newton strategy to find a solution starting from an initial guess. The central features of the initial guess for the first simulation were constant , and , as well as as a third order polynomial for that fulfills the boundary conditions. All remaining unknowns were approximated according to these assumptions. For subsequent simulations we used the results of previous simulations where appropriate. In those cases the adapted mesh from the previous simulation was coarsened to contain only half as many grid points as initially.
Iii Surface reaction chemistry: First-principles kinetic Monte Carlo simulations
The actual surface catalytic activity enters the FD simulations through the TOFs in the boundary value eq. (20) for the partial mass fluxes at the surface. The corresponding calculation of the rate of product formation per surface area and time from information about the elementary processes in the catalytic cycle is the realm of microkinetic theories. The most prominent such approach relies on phenomenological rate equations which only consider the mean-field averaged concentrations (coverages) of the reaction intermediates at the surfaceDumesic et al. (1993). This level of modeling of the surface chemistry is the prevalent standard in reactor engineeringKee et al. (2003); Deutschmann (2008). There, the kinetic quantities entering the rate expressions are in fact often treated as adjustable parameters. In a top-down fashion the idea is thus to use macroscopic reactor data to derive effective insight into the on-going surface catalytic activity. In more bottom-up oriented work the kinetic quantities are alternatively drawn from independent detailed experiments, or in modern hybrid approaches increasingly from first-principles calculationsGokhale et al. (2004); Maestri et al. (2009). The idea of such integrated approaches is correspondingly to model how both the intrinsic surface chemistry and transport effects contribute together to the macroscopically observable activity in a given reactor setup.
The latter is also the central objective of the present study. However, for the here aspired quantitative modeling present-day hybrid approaches are not sufficient. Use of scattered experimental and first-principles kinetic data from different sources and in the latter case potentially from different levels of approximate theory incurs a rather uncontrollable error. Even if the kinetic parameters of all involved elementary processes were reliable, there is still the error from the approximate mean-field treatment underlying the rate equation approach. In fact, for the here studied CO oxidation reaction at RuO(110) this error has recently been shown to be qualitative with deviations of the mean-field TOFs spanning up to several orders of magnitudeTemel et al. (2007). Aiming at an error-controlled multiscale modeling of predictive quality we therefore employ for the description of the surface kinetics 1p-kMC as most accurate approach with explicit account of the correlations, fluctuations and detailed spatial distributions of the chemicals at the surface.
iii.1 1p-kMC model of CO oxidation at RuO(110)
The molecular-level basis of 1p-kMC is a microscopically correct first-principles description of the elementary processes involved in the catalytic cycle. In the established model of CO oxidation at RuO(110) Reuter et al. (2004); Reuter and Scheffler (2006) this is specifically the set of 26 elementary processes defined by all non-correlated site and element specific adsorption, desorption, diffusion and reaction events that can occur on a lattice spanned by two different active sites offered by the surface, so-called bridge (br) and coordinatively unsaturated (cus) sites. For all these processes density-functional theory in conjunction with harmonic transition-state theory is used to compute the kinetic parameters. The resulting 26 first-principles rate constants form then the essential input for the actual 1p-kMC simulations which evaluate the statistical interplay among the surface chemical processes by following the long-term time evolution of the open catalytic system through numerical solution of a Markovian master equationReuter (2010).
Using exactly the computational setup as detailed beforeReuter et al. (2004); Reuter and Scheffler (2006); Temel et al. (2007), these simulations are carried out for the present purposes for a given local temperature and reactant partial pressures and directly at the surface, which in particular fixes the impingement and therewith the rate constants of the adsorption processes. Under such conditions, the system eventually reaches a unique steady-state, in which the detailed surface composition and occurrence of the individual elementary processes still exhibit the correct microscopic fluctuations, yet when averaged over simulation cells exceeding the characteristic correlation lengths at the surface have well-defined and constant values. These values thus comprise the average rate of reactant to product conversion under the given gas-phase impingement and local temperature, i.e. exactly the TOFs that enter the partial mass flux boundary condition of eq. (20). While it is only this averaged quantity that matters for the macroscopically described flow field, it is important to note that it is still properly derived from microscopic simulations that fully account for the site heterogeneity and distributions at the surface. This is thus distinctly different to the mentioned mean-field based phenomenological descriptions that are commonly integrated in the CFD modeling of macro-scale flow structures.
iii.2 Integration of 1p-kMC into the fluid dynamical environment
The 1p-kMC and FD simulations are intricately coupled. On the one hand, the TOFs required in eq. (20) to close the stagnation flow equations are provided by the 1p-kMC simulations. On the other hand, fixing the surface impingement in the 1p-kMC simulations requires the local temperature and gas-phase partial pressures directly at the surface, which are determined by the heat and mass transport modeled at the continuum level. A straightforward approach to this interdependence is a simultaneous solution until self-consistency between flow and 1p-kMC boundary condition is achievedVlachos (1997). For the here discussed stagnation flow equations this approach is in principle feasibleKissel-Osterrieder et al. (1998), albeit potentially numerically unstableChatterjee et al. (2004). However, for more complex reactor geometries it would quickly become intractable, as usually several independent 1p-kMC simulations would be required for every spatially resolved cell at the surface. Precisely due to the necessity to continuously rerun 1p-kMC simulations the approach would also be very inefficient for the here intended simulation of catalytic activity at a large variety of flow conditions.
We instead achieve a computationally much more efficient formulation by decoupling the interdependence through an instantaneous steady-state approximation Deutschmann (2008). The kMC simulations are first carried out to determine the steady-state TOFs for a wide range of surface impingement and local temperature conditions. The resulting grid data is then interpolated to a continuous representation, which in turn provides the entire necessary boundary condition for the stagnation flow problem. For the steady-state operation targeted in this manuscript this divide-and-conquer type approach is exact and may easily be applied to more complex reactor geometries. It could even be extended to transient phenomena under the assumption that on the time scale of relevant flux variations the surface chemistry adapts quasi instantaneously to the new steady-state, hence the name.
For the CO oxidation at RuO(110) we thus first computed 1p-kMC steady-state TOFs for the entire relevant range of temperatures and gas-phase composition. With a negligible CO readsorption probability these TOFs are independent of the CO partial pressure. Specifically we then used a dense grid with 25 K spacing for the temperature range 400 K 850 K and with logscale spacing to cover the range atm atm with 30 and the range atm atm with 42 spacings. Through modified quadratic Shepard interpolationRenka (1988a, b) this is converted into a reliable continuous representation TOF() that is finally presented as boundary condition to the stagnation flow solver.
The intrinsic activity resulting from the 1p-kMC model for the CO oxidation at RuO(110) has been analyzed for a wide range of temperature and pressure conditions before Reuter et al. (2004); Reuter and Scheffler (2006). Not surprisingly, high catalytic activity is only observed for a rather narrow range of gas-phase conditions, which stabilize O and CO simultaneously at the surface in appreciable amounts. For more O-rich feed the surface is poisoned by oxygen, for more CO-rich feed the surface is poisoned by CO, and little CO is formed in either case. For gas-phase conditions in the UHV regime, which allow direct measurements of the intrinsic activity, the model reproduces existing experimental TOF data quantitatively Reuter et al. (2004); Reuter and Scheffler (2006). It must be stressed though that for the description of the CO-poisoned regime the model has clear limitations. In corresponding CO-rich feeds the oxide surface would in reality eventually be reduced. By construction, this and the catalytic activity connected to such a phase transformation can not be grasped by the present 1p-kMC model assuming an intact underlying RuO(110) lattice.
Notwithstanding, this restriction concerns the modeling of the surface chemistry. The focus of the present work is instead to quantitatively integrate a given microkinetic description of this surface chemistry into a FD framework to assess heat and mass transport effects in a reactor geometry representative for in situ studies of model catalysts. For this endeavor the existing account of the surface chemistry is – despite its noted limitation – very suitable, in particular as it exhibits a number of features that we consider rather generic for a high-TOF reaction like the CO oxidation: (i) The intrinsic catalytic activity is narrowly peaked in a small range of gas-phase conditions. (ii) This activity is not sufficiently described by standard rate equation formulations. For predictive quality the first-principles based microkinetic modeling must therefore be based on an approach like 1p-kMC that explicitly accounts for the detailed spatial distribution of the chemicals at the surface. In turn, it is the latter type of approach to which the FD environment must be coupled, e.g. through the instantaneous steady-state approximation employed here. (iii) The peak activity at optimum partial pressures increases rapidly in the temperature range of interest. Towards the upper end at around 500-600 K and together with the high exothermicity of the CO oxidation reaction, this leads to a degree of mass conversion and heat release at the surface that is prone to transport limitations in the reactor.
As we will see the amount of heat dissipation possible at the back of the thin single-crystal is a crucial factor for such limitations that can mask the true intrinsic activity in in situ model catalyst studies. Not aiming (nor being able to) give a detailed account for one specific experimental setup we will analyze this in the following in more generic terms for the two already described opposite limits. In the adiabatic limit there is no heat flux at all through the sample backside, mimicking to some degree the situation that could e.g. arise from an insulating sample holder. In contrast, in the isothermal limit we assume the sample to be sufficiently well connected to the outside system that a constant nominal temperature is maintained throughout. Here, this is chosen to be the same temperature as that of the gases at the inlet.
Considering the peaked structure of the intrinsic catalytic activity in -space we can conveniently study heat and mass transfer effects in these two limits focusing on prototypical sets of gas-phase conditions: For defined temperature, essentially zero CO concentration ( atm throughout) and near-ambient oxygen partial pressure at the inlet these sets comprise a range of inlet CO partial pressures. They cover the O-poisoned regime at the lowest , the CO-poisoned regime at the highest , and span over the conditions of highest intrinsic activity.
iv.1 Adiabatic limit
iv.1.1 Surface heating
The heat flux through the back of the sample is completely suppressed in the adiabatic limit. The only dissipation channel left for the heat released by the exothermic surface reactions is then into the surrounding gas phase. Compared to the heat conduction through a metallic sample this channel is rather inefficient. One can therefore suspect that it may not be efficient enough to maintain the nominal surface temperature, once the TOF and therewith the generated heat rate exceeds a certain critical value. Instead, the surface will heat up and give rise to gas-phase temperature gradients from inlet to active surface. This is indeed what we find in the coupled 1p-kMC–FD simulations. When using representative parameters for the inlet distance cm and axial inlet flow velocity cm/s significant deviations from the nominal surface temperature set in for near-ambient environments at TOFs exceeding site s. Such peak TOFs are reached for optimum partial pressure conditions at temperatures above about 500 K.
The effect of the ensuing temperature increase at the active surface on the observable steady-state conversion rates is quite dramatic already at this threshold temperature and is illustrated in Fig. 2. While the intrinsic activity is peaked close to stoichiometric feeds, the heat transport limitations lead to the stabilization of a high-activity operation mode that extends to significantly more CO-rich conditions. In this branch of the observable TOF-profile shown in Fig. 2 the surface temperature is up to 150 K higher than the nominal temperature at the inlet. Quantitatively this and the concomitant TOFs in the newly established high-activity mode depend on the details of the reactor setup. This is exemplified in Fig. 2 with the TOF-profiles that result when increasing or decreasing the inlet distance by one order of magnitude. For a larger cm the extension of the high-activity branch is reduced at overall lower conversion rates, while for a smaller mm the branch extends to much higher CO partial pressures and the observable conversion rate exceeds the peak intrinsic activity by more than one order of magnitude. Similar, but quantitatively smaller variations are obtained when changing the inlet velocity by one order of magnitude up or down (not shown). For a lower mm/s the high-activity branch exhibits slightly smaller TOFs and extends up to slightly smaller than for the cm/s displayed in Fig. 2. An increase to cm/s, on the other hand, increases the extension of the branch to higher CO partial pressures at also higher TOFs than those shown in Fig. 2.
Regardless of these quantitative variations, the net effect of the surface heating resulting from the transport limitations is in all cases a substantially changed observable TOF-profile compared to the true underlying intrinsic reactivity. The absolute TOFs in the relevant high-activity regime are significantly different and the inlet gas-phase conditions for which highest conversions are obtained are shifted to much more CO-rich feeds. Obviously, if these observable TOFs were mistaken for the intrinsic reactivity, wrong mechanistic conclusions about the on-going surface chemistry in such in situ environments would be derived. Furthermore, the observable TOF-profile in Fig. 2 exhibits another feature that is completely absent in the intrinsic reactivity: For a range of CO partial pressures we obtain two stationary solutions: The high-activity branch and in addition a low-activity branch that coincides with the intrinsic activity. As the underlying 1p-kMC model of CO oxidation at RuO(110) has no multiple steady-statesTemel et al. (2007), this bistability arises solely from the coupling of macroscopic transport and surface chemistry.
iv.1.2 Formation of a boundary layer
An analysis of the observed range of transport effects (high-activity branch with concomitant bistability and variations with reactor setup) starts from the anticipated formation of a finite boundary layer above the flat-faced model catalyst. As schematically drawn in Fig. 3 non-vanishing gradients of temperature and partial pressures are in general restricted to this boundary layer. This is due to the convective nature of the transport. The gas streams towards the surface and the chemical reactions at the surface have no influence on the fields far away. Heat and species move with the flow, which dominates over any non-convective transport. Only when the axial velocity approaches zero near the surface, diffusion and heat conduction kick in, as they are now of similar size as the convective transport. For the here discussed adiabatic limit the continuous catalytic conversions at the surface lead to a continuous heat release into the gas phase. For sufficiently large TOFs above the critical value this heat rate is too high to be efficiently transported away by the flow. The result is a steady-state temperature profile above the surface as sketched in Fig. 3, where the temperature rises within the boundary layer from to the values for shown in Fig. 2. On the other hand, surface mass conversions of the order of magnitude as those of Fig. 2 are still too low to significantly affect the gas composition in the boundary layer. The formed products are transported away sufficiently quickly so that the partial pressures remain essentially unchanged, i.e. we find for all conditions discussed in Fig. 2.
The latter feature allows us to suitably discuss the observed transport effects in terms of the intrinsic reactivity summarized in Fig. 4. Shown is the TOF profile as obtained with 1p-kMC for the same atm as in Fig. 2 and as a function of surface temperature and CO partial pressure. The region of highest-activity is narrowly peaked and the peak TOFs at this rim increase steadily with temperature. A crucial feature for the following is that the rim does not go along constant partial pressure ratio for higher temperatures (i.e. vertically up in Fig. 4), but shifts continuously towards more CO-rich mixtures (i.e. diagonally up in Fig. 4). This is because the most active state is characterized by the coexistence of O and CO at the surface, which follows more a constant ratio of chemical potentials than partial pressures. Imagine now that the system starts at the nominal inlet and surface temperature of 500 K, corresponding to the bottom horizontal line in Fig. 4. As long as the intrinsic TOF does not exceed the critical value, the system is able to maintain this temperature and the observable TOF is identical to the intrinsic one, cf. Fig. 2. For the range 0.4 atm atm the activity is, however, above the critical value. The system cannot transport the generated heat away sufficiently quickly, and surface and boundary layer start to heat up. As directly apparent by focusing on e.g. atm in this critical pressure range the intrinsic TOFs increase with increasing surface temperature (vertically up in Fig. 4). This gives rise to a runaway process. Higher TOFs generate more heat, which increases the surface temperature and therewith leads to even higher TOFs. This will only stop once the system is over the highest activity rim in Fig. 4. Further increase in surface temperature (i.e. moving even more vertically up along the line of constant in Fig. 4) leads then to a decrease in intrinsic reactivity. This allows the system to find a new steady-state, viz. the high-activity branch corresponding to the increased surface temperature in Fig. 2.
Where on this downward TOF slope the system precisely settles down depends on how efficiently the generated heat can be transported away, i.e. which surface temperature results for the heat rate connected with a given intrinsic TOF. This is crucially controlled by the thickness of the boundary layer and corresponding thickness variations rationalize the entire observed dependence on the reactor setup. A smaller inlet distance compresses the boundary layer and therewith enables a better heat convection. Correspondingly, at smaller the system ends up closer to the most active rim in Fig. 4 and the high-activity operation mode exhibits higher observable TOFs as seen in Fig. 2. A higher inlet velocity has the same effect on the boundary layer, and therewith also leads to higher observable TOFs as described above.
The intrinsic TOF profile shown in Fig. 4 also helps to understand why the high-activity branches eventually break down and how their extension varies with the reactor setup. Apart from the highest-activity rim dominated by the oxidation reactions between CO and O adsorbed at the most active cus sites one can discern at the upper edge of Fig. 4 a weak new rise of the intrinsic TOFs. At these highest temperatures shown this is due to the ”entropic” widening of the kinetic ”phase” transition region where the appropriate chemical potential ratio in the gas-phase stabilizes the coexistence of both reactants at the surface Reuter and Scheffler (2006, 2003). When over the rim the hitherto described decrease of the intrinsic TOF with increasing temperature will therefore eventually change into a new increase. If the system reaches this change of slope, a new runaway cycle of increasing temperature and TOF will start and no stationary operation mode can be stabilized (at least not in the temperature range of interest for the present study). Correspondingly, in Fig. 4 the high-activity branches always break down at positions of such a gradient change in the intrinsic TOF profile. As apparent from Fig. 4 the further away from the rim the high-activity branch is situated, the earlier it hits this slope change, i.e. its extension reaches only up to smaller . In the present adiabatic limit, modifications of the reactor setup that compress the boundary layer (either by higher axial inlet velocity or smaller inlet distance) lead therefore to higher absolute observable TOFs in a high-activity branch that extends up to higher CO partial pressures.
Finally, some remarks about the observed bistability are at place. The structure of the TOF profile in Fig. 4 rationalizes why for some pressure conditions two steady-state solutions are obtained. One, in which the system exhibits a low activity that coincides with the intrinsic one, and one, in which significant surface heating has brought the system above the highest-activity rim. While intuitive, the rationalization in terms of thermal runaway is at present clearly an interpretation. A verification would require the extension of the present steady-state approach to transient operation. Only corresponding time-resolved simulations will then give access to the wealth of phenomena that are now only suggested by the observed bistability. Notably this is the possibility of oscillations between the two modes. In contrast to e.g. purely surface reaction–surface diffusion driven oscillations on single-crystals in UHV Imbihl (2009) the mechanistic details behind corresponding reactor–reaction oscillations in the ambient pressure regime are only poorly understood Schüth et al. (1993). Obviously, extending the present model in this direction offers the prospect of a detailed analysis, on which we will concentrate in future work.
iv.2 Isothermal limit
In the opposite isothermal limit the high thermal conductivity of the metallic sample allows for such an efficient removal of the generated heat that even at higher temperatures around 600 K, where the intrinsic peak TOFs at optimum partial pressures exceed sites, no significant surface heating results. The temperature remains at the nominal inlet value throughout the entire system. As shown in Fig. 5 the intrinsic activity is nevertheless noticeably masked, this time by mass transfer limitations in the boundary layer. At cm/s and cm the maximum observable TOFs are lower than the peak intrinsic ones, and a high-activity branch extends now to much more CO-poor feeds. As in the adiabatic limit there is a range of CO partial pressures for which we observe a bistability, and the results depend again quantitatively on the reactor setup. This is illustrated in Fig. 5 by showing the data also for an increased ( cm) and decreased ( mm) inlet distance. For smaller inlet distances higher observable TOFs result. Highly comparable variations are obtained when changing the axial inlet velocity by one order of magnitude up or down, with smaller yielding larger observable TOFs (not shown). The varying reactor conditions also affect the extension of the high-activity branch, which at maximum reaches down to atm, i.e with atm to stoichiometric feed. In addition, there is now in principle also a dependence of the results on the employed sample thickness , which enters through the boundary condition, eq. (23). However, due to the high thermal conductivity this dependence is in practice negligible. Compared to the results in Fig. 5, which were obtained using mm, changing the thickness by one order of magnitude up or down has virtually no effect on the observable catalytic activity.
iv.2.1 Mass transfer limitations
The entire range of transport effects observed in the isothermal limit is this time due to mass transfer limitations in the boundary layer. At the high intrinsic TOFs around peak performance the mass conversion at the active surface is so large that these limitations lead to noticeable changes of the partial pressure profiles from inlet to surface. As schematically shown in Fig. 3 there is essentially a build-up of a significant product concentration that is no longer sufficiently quickly removed. This goes hand in hand, cf. eq. (8), with a decrease in reactant partial pressures, i.e. O and CO are hindered in their access to the active surface. Due to the similar transport parameters and mass of both diatomic reactants, cf. Table 1, this limitation affects both species similarly. As long as the nominal inlet composition is different from stoichiometric feed, a corresponding roughly equal reduction of both reactant partial pressures close to the surface will then effectively change the ratio. As also apparent from Fig. 4 at K the range of peak intrinsic activity corresponds to quite CO-rich feeds. In this range the mass transfer limitations therefore lead to a noticeable increase of the ratio compared to the nominal inlet composition as shown in Fig. 5. At a nominal inlet composition that would correspond to optimum intrinsic activity, atm in Fig. 5, the surface then sees a comparatively more CO-rich feed and the observable TOFs are lowered compared to the intrinsic ones. On the other hand, at a nominal inlet composition only slightly more CO-rich than stoichiometric feed, where the intrinsic activity would already have collapsed in Fig. 5, the significantly more CO-rich feed effectively seen by the surface corresponds in fact to conditions close to optimum intrinsic activity. The observable TOF is much increased and the high-activity branch of Fig. 5 results. Exactly at stoichiometric feed this effective CO enrichment close to the surface ends, and for even more CO-poor mixtures possible mass transfer limitations would rather suppress the CO minority species. However, for such partial pressure ratios the intrinsic activity is low anyway and no mass transfer limitations arise. At the latest the high-activity branch therefore breaks down at stoichiometric feed. With this understanding the observed variations with the reactor setup are also easy to rationalize. A smaller boundary layer as resulting from increased axial inlet velocity or reduced inlet distance reduces the mass transfer limitations. The partial pressure ratio at the surface gets closer to the nominal one. In turn, the observable TOFs approach the intrinsic ones as illustrated in Fig. 5 for the varying inlet distances.
In the presence of such mass transfer limitations a natural question is to what degree they mask the intrinsic surface reactivity. Is the observable TOF profile the result of a complex mixture of the on-going surface chemistry and the gas-phase transport, or are the flow conditions in the reactor such that the measured activity conveys little information about the actual catalyst anymore. To qualify this it is useful to assess the upper TOF limit that results if mass flux is the only limitation. Such an estimate can be obtained by realizing that the steady-state mass conversion by the catalyst can never be higher than the one that completely depletes the minority species at the surface. Rather than using the catalyst specific boundary condition eq. (20) that depends on the actual intrinsic TOFs, we then simply employ for O-rich feeds
and for CO-rich feeds
For the respective majority species the boundary condition is still eq. (20), but the TOF entering this equation is now determined by the mass flux for the minority species, i.e. the conversion is completely dictated by the amount of impinging minority species. Figure 6 shows the upper TOF limit that results from this estimate for the afore discussed gas-phase conditions of Fig. 5. Apparently, the observable TOFs come very close to this upper limit for most of the active region. This means that in this regime the measurable profile has very little to do with the actual RuO(110) catalyst, it rather reflects only the flow conditions in the employed reactor. Obviously, and similar to the adiabatic limit discussed before, if corresponding effects are not appropriately accounted for in in situ studies, wrong conclusions about the surface chemistry at technologically relevant gas-phase conditions will be derived.
V Summary and conclusions
In summary we have presented an efficient approach to couple first-principles kinetic Monte Carlo and fluid dynamical simulations in the context of heterogeneous catalysis. This augments the accurate description of the surface chemistry achieved by 1p-kMC with a continuum account of the heat and mass transfer in a given reactor geometry, or vice versa it integrates the accurate 1p-kMC microkinetics into reactor-level modeling. In prevalent chemical engineering approaches the latter modeling instead incorporates phenomenological microkinetic treatments based on mean-field rate equations. In contrast to such a description, the presented 1p-kMC based multiscale modeling approach derives the average flux quantities required for the macroscopically described flow field properly from microscopic simulations that fully account for the site heterogeneity and distributions at the catalyst surface. As such it has the potential to carry the predictive power of the underlying electronic structure calculations for the elementary processes all the way to the reactor level.
On the way to such a first-principles chemical engineering we have applied the approach to the problem of in situ studies of model catalysts using the ambient pressure CO oxidation at RuO(110) as a representative example. As a suitable, though idealized reactor model to discuss the transport at the flat-faced surface we have chosen a stagnation flow geometry. The observed catalytic function depends sensitively on the employed reactor geometry (mimicked in this study by varying the inlet-surface distance) and the applied throughput rate (i.e. the streaming velocity at the inlet). For the thin, well heat conducting single crystal the degree of heat dissipation possible at the back of the sample (e.g. through radiative loss or contact to the sample holder) is a further crucial factor. Not aiming, nor being able to address specific experimental realizations this was considered through two opposite extremes: The isothermal limit mimicking a highly efficient heat coupling of the crystal to the system, and the adiabatic limit to represent a well insulated sample. In both limits transport effects were found to readily mask the intrinsic catalytic function at the high conversion rates reached at near-ambient gas-phase conditions. In the adiabatic limit this is due to a significant surface heating, in the isothermal limit due to mass transfer limitations that lead to the build-up of a significant concentration of products in the boundary layer above the active surface. With the single-crystal in real experimental setups neither perfectly heat-coupled nor isolated, these two effects discussed here separately will obviously be intricately intermingled and need to be disentangled by dedicated measurements and setups. Furthermore, we obtained in both limits a range of gas-phase conditions where the system exhibits two stationary operation modes, a low-activity branch corresponding to the intrinsic reactivity and a high-activity branch which arises from the coupling of the surface chemistry to the surrounding flow field. A corresponding bistability obtained here in the steady-state limit clearly suggests that the system could oscillate between the two modes, possibly even inhomogeneously in form of reaction fronts over the single-crystal surface. In case of heat transfer limitations, an intuitive propagation mechanism would hereby be via the formation of local hot spots, while in the mass transfer case it would be via gas-phase coupling, with the presented approach establishing the intriguing possibility to quantify these model conceptions with first-principles based simulations.
The main objective of in situ studies of model catalysts is a detailed, atomic-scale analysis of the catalytic function at technologically relevant gas-phase conditions, thereby bridging the pressure gap to the at present often much better characterized function in UHV. The range of transport effects discussed in the present study qualifies the additional complexity that needs to be accounted for in corresponding work to prevent wrong mechanistic conclusions about the surface chemistry at high pressure. That this complexity has potentially not yet been sufficiently appreciated may very well be the reason for the many existing controversies in the field. Also because of the limitation of the employed 1p-kMC RuO(110) model with respect to a reduction of the oxide catalyst we have refrained from comparing our simulations to already published experimental data. Nevertheless, we note that the gas-phase conditions and TOFs discussed here are of the order of magnitude presented in a manifold of in situ studies of CO oxidation at late transition metal catalysts. In this respect it is important to recognize that the CO oxidation reaction – that has been a fruitfly reaction in UHV Surface Science due to its alleged ”simplicity” and model character – requires particular attention. The high turnovers that can be reached precisely because of this ”simplicity” make this reaction much more prone to transport effects than other more complex, selective ones.
Acknowledgements.Funding within the MPG Innovation Initiative Multiscale Materials Modeling of Condensed Matter and the DFG Cluster of Excellence Unifying Concepts in Catalysis is gratefully acknowledged.
- Reuter (2010) K. Reuter, in Modeling Heterogeneous Catalytic Reactions: From the Molecular Process to the Technical System, edited by O. Deutschmann (Wiley-VCH, Weinheim, 2010).
- Meskine et al. (2009) H. Meskine, S. Matera, M. Scheffler, K. Reuter, and H. Metiu, Surf. Sci. 603, 1724 (2009).
- Matera and Reuter (2009) S. Matera and K. Reuter, Catal. Lett. 133, 156 (2009).
- Stierle and Molenbroek (2007) A. Stierle and A. Molenbroek, MRS Bull. 32, 1001 (2007).
- Reuter et al. (2004) K. Reuter, D. Frenkel, and M. Scheffler, Phys. Rev. Lett. 93, 116105 (2004).
- Reuter and Scheffler (2006) K. Reuter and M. Scheffler, Phys. Rev. B 73, 045433 (2006).
- Kee et al. (2003) R. Kee, M. Coltrin, and P. Glarborg, Chemically Reacting Flows (Wiley, Hoboken NJ, 2003).
- Cloutman (2000) L. Cloutman, Tech. Rep. UCRL-ID-139893, Lawrence Livermore National Laboratory (2000).
- Linstrom and Mallard (2009) P. Linstrom and W. Mallard, eds., NIST Chemistry WebBook, NIST Standard Reference Database Number 69 (National Institute of Standards and Technology, 2009), URL http://webbook.nist.gov.
- Stull and Prophet (1971) D. Stull and H. Prophet, JANAF Thermochemical Tables (U.S. Bureau of Standards, Washington D.C., 1971), 2nd ed.
- Bird et al. (1960) R. Bird, W. Stewart, and E. Lightfood, Transport Phenomena (Wiley, Hoboken NJ, 1960), 1st ed.
- Deutschmann (2008) O. Deutschmann, in Handbook of Heterogeneous Catalysis, edited by G. Ertl, H. Knözinger, F. Schüth, and J. Weinkamp (Wiley-VCH, Weinheim, 2008), 2nd ed.
- Müller (1985) I. Müller, Thermodynamics (Pitman Publishing, 1985), 1st ed.
- Ferizovic et al. (2009) D. Ferizovic, L. K. Hussey, Y.-S. Huang, and M. Munoz, Appl. Phys. Lett. 94, 131913 (2009).
- Lide (1995) D. R. Lide, ed., Handbook of Chemistry and Physics (CRC, 1995), 76th ed.
- Ascher and Spiteri (1994) U. Ascher and R. Spiteri, Siam J. Sci. Comp. 15, 938 (1994).
- Dumesic et al. (1993) J. A. Dumesic, D. F. Rudd, and L. M. Aparicio, The Microkinetics of Heterogeneous Catalysis (American Chemical Society, Washington, DC, 1993).
- Gokhale et al. (2004) A. A. Gokhale, S. Kandoi, J. P. Greeley, M. Mavrikakis, and J. A. Dumesic, Chem. Eng. Sci. 59, 4679 (2004).
- Maestri et al. (2009) M. Maestri, D. G. Vlachos, A. Beretta, G. Groppi, and E. Tronconi, AICHE J. 55, 993 (2009).
- Temel et al. (2007) B. Temel, H. Meskine, K. Reuter, H. Metiu, and M. Scheffler, J. Chem. Phys. 126, 204711 (2007).
- Vlachos (1997) D. Vlachos, AIChE J. 43, 3031 (1997).
- Kissel-Osterrieder et al. (1998) R. Kissel-Osterrieder, F. Behrendt, and J. Warnatz, Proc. Comb. Inst. 27, 2267 (1998).
- Chatterjee et al. (2004) A. Chatterjee, M. Snyder, and D. Vlachos, Chem. Eng. Sci. 59, 5559 (2004).
- Renka (1988a) R. Renka, ACM Trans. Math. Software 14, 139 (1988a).
- Renka (1988b) R. Renka, ACM Trans. Math. Software 14, 151 (1988b).
- Reuter and Scheffler (2003) K. Reuter and M. Scheffler, Phys. Rev. B 68, 045407 (2003).
- Imbihl (2009) R. Imbihl, Surf. Sci. 603, 1671 (2009).
- Schüth et al. (1993) F. Schüth, E. Henry, and L. Schmidt, Adv. Catal. 39, 51 (1993).