Micro-Kinetic Modelling of CO-TPD from Fe(100)—Incorporating Lateral Interactions

Abstract

The experimentally determined temperature programmed desorption profile of CO from Fe(100) is characterized by four maxima, i.e., α₁-CO, α₂-CO, α₃-CO, and β-CO (see e.g., Moon et al., Surf. Sci. 1985, 163, 215). The CO-TPD profile is modeled using mean-field techniques and kinetic Monte Carlo to show the importance of lateral interactions in the appearance of the CO-TPD-profile. The inclusion of lateral interactions results in the appearance of a new maximum in the simulated CO-TPD profile if modeled using the mean-field, quasi-chemical approach or kinetic Monte Carlo. It is argued that α₂-CO may thus originate from lateral interactions rather than a differently bound CO on Fe(100). A detailed sensitivity analysis of the effect of the strength of the lateral interactions between the species involved (CO, C, and O), and the choice of the transition state, which affects the activation energy for CO dissociation, and the energy barrier for diffusion on the CO-TPD profile is presented.

1. Introduction

Heterogeneous catalysts are the workhorses for the chemical process industry. They provide an alternative reaction pathway for the transformation of the reactants into desired products. Hence, catalysis is a kinetic phenomenon, and the search of improved catalysts and catalytic processes requires an in-depth understanding of the rate of the underlying processes on the catalyst surface. The overall rate of catalytic processes links the observed rate of reaction to the concentration of the species in the gas phase considering geometric factors of the catalyst to include possible mass transfer resistances. In its simplest form, the rate law is in the form of a power-rate law. This will only lead to a limited understanding of factors determining the catalyst performance.

Most catalytic reactions are modeled most commonly with a Langmuir–Hinshelwood mechanism or an Eley–Rideal type of mechanism [1], in which reactants (partially) in the adsorbed state interact with each other in a series of elementary reaction steps. The overall rate of reaction is then obtained by applying the quasi-steady-state approximation (QSSA) and the assumptions associated with the Langmuir isotherm.

The heterogeneous catalysis community has obtained a basic understanding of events occurring on surfaces using the Langmuir isotherm (not unlike the understanding the 19^th-century physicists of gases by using the ideal gas law). However, it is well known that adsorbates are affected by the presence of other adsorbates due to lateral interactions [2,3,4,5,6,7,8,9,10], which originate from so-called through-space and through-surface interactions [11]. The classical through-space interactions can be long range interactions, e.g., via induced dipole-dipole interactions [12], whereas the through-surface interactions are mainly influenced by nearby neighbor interactions. Hence, lateral interactions are important at high surface coverage, i.e., at conditions often obtained under industrially relevant conditions, and will need to be incorporated in an appropriate manner in the kinetic models [13,14,15,16]. Monte Carlo kinetic modeling seems to be the preferred option [14,15,16], although the widely different time scales of the various processes taking place may put some constraint on its use and simplified models may result in more transparent results.

Temperature programmed desorption (TPD) of CO from Fe(100) has been well described [17,18,19,20]: after adsorption of CO at 150 K the desorption profile shows four different states for adsorbed CO, viz. three different states for associatively adsorbed CO (with a maximum desorption rate at 220 K for α₁-CO, at 306 K for α₂-CO and at 420 K for α₃-CO) and dissociated CO (β-CO) with a maximum desorption rate at 820 K when using a heating rate of 10 K·s⁻¹ [18]. Molecular CO bonded in its most stable state, i.e., in the 4-fold hollow site, has been linked to α₃-CO [6,7,19,20,21,22]. CO located in the atop position has been linked to α₁-CO [19,23]. The desorption of carbon monoxide from Fe(100) associated with α₂-CO has been tentatively assigned to CO bonded to a bridge site [19], although this site is unstable at low CO-coverage and is a local minimum at high coverage, which is less stable than adsorption on the atop site [6].

In this study, we focus on incorporating lateral interactions in various kinetic models to describe temperature programmed desorption of CO from Fe(100). This system is of interest due to its relevance to the industrial Fischer–Tropsch process and may thus yield insight into the importance of lateral interactions under those conditions.

2. Results and Discussion

2.1. CO-Adsorption on p(4x4)-Fe(100)

The adsorption of CO on Fe(100) has been studied extensively using DFT (see e.g., [6,8]) and the adsorption energy is known to decrease with increasing coverage. Here we report on the adsorption of CO on p(4x4)-Fe(100). At a coverage of CO of 0.0625 monolayer (ML), the C-O bond length was found to be 1.304 Å with a tilt angle of 46.8° from the normal and a stretching frequency of 1169 cm⁻¹ was observed (in reasonable agreement with the experimentally reported tilt angle of 45° ± 10° and a stretching frequency of 1210 cm⁻¹ [1,20,21,22]). The Fe-C bond distance was determined to be 2.16 Å. The adsorption energy of a single CO on this unit cell (Θ_CO = 0.0625 ML) was found to be −1.91 eV/CO. Adsorption of CO leads to some deformation of the unit cell (primarily in the top layer) and deformation energy of 0.27 eV for the p(4x4)-Fe(100) unit cell was determined at this coverage.

The p(4x4)-Fe(100) cell was used to probe the multiplicity of possible configurations of adsorbed CO at a constant coverage of Θ_CO = 0.25 ML (see Figure 1a). The configurations differ in the number of shared surface iron atoms and the orientation of adsorbed CO-species.

The “diagonal” configuration corresponds to a c(2x2) surface structure for CO on Fe(100). In this configuration, CO with a bond length of 1.304 Å was found to have a Fe-C bond length of 2.11 Å and a tilt angle of 47.0° away from the normal. The adsorption energy was found to be slightly lower than for CO at Θ_CO = 0.0625 ML (−1.93 eV/CO vs. 1.91 eV/CO). A similar stabilization in this configuration has been observed for the co-adsorption of C-C, C-O, C-H, C-CH, C-CH₂ and C-CH₄ interactions [24]. This small stabilization of adsorbed CO may be attributed to reduced deformation energy of the surface structure per adsorbed CO (0.098 eV/CO), which offsets the repulsive interaction between adsorbed CO. The surface layer is slightly deformed upon adsorption of CO and the position of the surface iron atoms directly interacting with adsorbed CO shifts towards the adsorbed species. The presence of more CO on the surface will reduce the shift of the iron surface atoms.

In the “corner share” configuration, adsorbed CO takes in two unique geometries, each with a C-O bond length of 1.301 Å. These geometries differ in the Fe-C bond length (2.10 vs. 2.12 Å) and the tilt angle from the normal (49.2° vs. 44.9°).

The “0 share” configuration corresponds to the structure typically calculated on p(2x2)-Fe(100) and the obtained adsorption structure. Both the calculated energy and geometry compare well with the previously reported one using the same functional [7]. The adsorption energy and the deformation energy for the “diagonal”, “corner share” and “0 share” configurations do not differ much and may be considered the same within the error of calculation.

A significant weakening of the adsorption of CO was observed in the “2 share” configuration with an adsorption energy of −1.79 eV/CO. The increase in the deformation energy contributes only ca. 0.01 eV to the weakening of the adsorption of CO. In this configuration, adjacent CO adsorbates are orientated perpendicular to each other’s C-O bond. The adsorbate-adsorbate interaction may have resulted in a re-orientation of the adsorbate structure [25]. Two unique geometries for adsorbed CO were obtained with Fe-C bond length of 2.10 Å and 2.15 Å and with a tilt angle from the normal of 47.7° and 44.4° (the C-O bond lengths were equal at 1.300 Å).

The “4 share” configuration corresponds to a full coverage on 25% of the p(4x4)-Fe(100) unit cell and zero coverage at 75% of the unit cell. Adsorbed CO has now adjacent CO with the C-O bond perpendicular to the C-O bond, both C-O bonds parallel to each other and C-O bonds directing in the opposite direction. In this configuration, the Fe-C bond is elongated (2.21 Å) and the C-O bond distance shortened to 1.285 Å in comparison to the other configurations. The tilt angle of C-O bond from the surface normal is now 50.2°.

The adsorption at Θ_CO = 0.25 ML is weakest in the “cluster” configuration (−1.69 eV/CO) with all adsorbed CO in a row and the C-O bonds in parallel. There are still two unique geometries for adsorbed CO (both with a C-O bond length of 1.291 Å), viz. one with a Fe-C bond length of 2.11 Å and a tilt angle from the normal of 45.9° and one with a Fe-C bond length of 2.15 Å and a tilt angle from the normal of 42.9°.

A few empirical observations can be made on the adsorption energy for the adsorption of CO in different configurations. Configurations with CO sharing none or a single iron surface atom, i.e., in the ‘diagonal’, ‘corner share’ or ‘0 share’ configuration, the adsorption energy does not seem to be affected significantly (<0.02 eV/CO) by the presence of nearby adsorbed CO (and seems similar to the adsorption energy for the adsorption of a single CO on p(4x4)-Fe(100)). The configuration where CO adsorbates are in adjacent hollow sites and lie perpendicular to each other causes a decrease in the adsorption energy by 0.18 eV per interaction. The configuration where CO adsorbates are in adjacent hollow sites and lie parallel to each other causes a decrease in the adsorption energy by 0.24 eV per interaction. The configuration where CO adsorbates are in adjacent hollow sites and tilted away from each other causes a decrease in the adsorption energy by 0.10 eV per interaction. This empirical observation was tested with four additional configurations (see Figure 1b) showing a reasonable description of the adsorption energy (in particular seeing the large difference in the computational time) allowing a quick estimation of the adsorption energy for all 805 possible unique configurations on p(4x4)-Fe(100) at all possible coverages.

Lateral interactions are in principle temperature dependent. However, the difference in the enthalpy between 0 K and 500 K between the ‘diagonal’ configuration, and the ‘cluster’ configuration for CO adsorbed on Fe(100) at Θ_CO = 0.25 ML is only 3 meV/CO and the difference in the entropy between 0 K and 500 K is only 1.35x10⁻² meV/K/CO. Hence, the influence of the temperature variation on the lateral interaction may be neglected.

2.2. Co-Adsorption of CO-C, CO-O, C-C, C-O, and O-O on Fe(100)

Interactions between nearest neighbor and next-nearest neighbor co-adsorbed species both on their most stable sites were considered (adsorbed CO, C, and O all prefer the four-fold hollow site on Fe(100) [6,8,24,26]). Three different types of configurations are considered to investigate these “A-B” interactions, viz. two configurations on a p(2x2) cell at a total coverage of 0.5 ML, in which the adsorbates are either diagonally aligned (diagonal configuration in Figure 1) or adjacent to each other (analogue to the cluster configuration in Figure 1), and one configuration on a p(3x2) cell at a total coverage of 0.33 ML. Each configuration has equal amounts of species A and B. In a configuration on the p(3x2) cell at a coverage of 0.33 ML, each adsorbate only has one nearest neighbor while for the 0.5 ML adjacent configuration on a p(2x2) cell, each adsorbate has two nearest neighbors. Each adsorbate has two next-nearest neighbors in the diagonal configuration at a coverage of 0.5 ML on a p(2x2) cell.

The adsorption energy for atomic carbon on a p(2x2)-Fe(100) unit cell at a coverage of 0.25 ML was determined to be −8.00 eV/C relative to atomic carbon in the gas phase in good agreement with previously reported values [24,27,28,29,30]. The Fe-C distance was found to be 1.961 Å, which decreases with increasing coverage (see

Table 2. Nearest and next-nearest neighbor interactions in the system containing CO, C, and O on Fe(100).

Configuration		0.33 ML 1	0.5 ML 2 Adjacent	0.5 ML2 Diagonal
C-C	dFe-C, Å	1.955	1.931	1.898
	Edeformation 3, meV	116	150	103
	Eexcess 4, meV	−50	131	−68
	ΔHexcess (1000K) 5, meV	3	2	24
	TΔSexcess (1000 K) 5, meV	−24	−16	52
O-O	dFe-O, Å	2.100	2.122	2.079
	Edeformation 3, meV	224	49	134
	Eexcess 4, meV	−46	380	0
	ΔHexcess (1000 K) 5, meV	1	0	6
	TΔSexcess (1000 K) 5, meV	−28	−2	−3
CO-C	dFe-C(O), Å	2.109	2.151	2.015
	dC- O, Å	1.317	1.287	1.317
	dFe-C, Å	1.931	1.962	1.937
	Edeformation 3, meV	448	172	220
	Eexcess 4, meV	38	287	−200
	ΔHexcess (1000 K) 5, meV	9	22	35
	TΔSexcess (1000 K) 5, meV	−141	−88	−20
CO-O	dFe-C(O), Å	2.127	2.134	2.031
	dC- O, Å	1.299	1.286	1.319
	dFe-O, Å	2.125	2.107	2.097
	Edeformation 3, meV	85	61	116
	Eexcess 4, meV	24	387	−124
	ΔHexcess (1000 K) 5, meV	3	−13	−4
	TΔSexcess (1000 K) 5, meV	13	−77	−42
C-O	dFe-C, Å	1.931	1.917	1.931
	dFe-O, Å	2.160	2.129	2.049
	Edeformation 3, meV	448	171	220
	Eexcess 4, meV	185	640	−190
	ΔHexcess (1000 K) 5, meV	4	2	17
	TΔSexcess (1000 K) 5, meV	−56	0	2

¹ adsorption on a p(3x2)-Fe(100) unit cell, ² adsorption on a p(2x2)-Fe(100) unit cell, ³ deformation energy of a unit cell: ${\mathrm{E}}_{\mathrm{deformation}}={\mathrm{E}}_{\left(\mathrm{Fe}\left(100\right)+{\displaystyle \sum }_i{\mathrm{adsorbate}}_i\right)-{\displaystyle \sum }_i{\mathrm{adsorbate}}_i}-{\mathrm{E}}_{\mathrm{Fe}\left(100\right)}$; ⁴ excess energy relative to individually adsorbed species at Θ = 0.25 ML: ${\mathrm{E}}_{\mathrm{excess}}=\left({\mathrm{E}}_{\mathrm{Fe}\left(100\right)+{\displaystyle \sum }_i{\mathrm{adsorbate}}_i}-{\mathrm{E}}_{\mathrm{Fe}\left(100\right)}\right)-{\displaystyle \sum }_i\left({\mathrm{E}}_{\mathrm{p}\left(2x2\right)-Fe\left(100\right)+{\mathrm{adsorbate}}_i\left(\Theta =0.25ML\right)}-{\mathrm{E}}_{\mathrm{p}\left(2x2\right)-Fe\left(100\right)}\right)$; ⁵ change in the enthalpy and entropy contribution at 1000 K to the energy of the co-adsorbed state and the state in which the adsorbates are adsorbed individually at a coverage of 0.25 ML on p(2x2)-Fe(100).

). The adsorption of atomic carbon results in some deformation in particular of the surface layer, however, the deformation energy does not scale with the number of atomic carbon species in a unit cell. The lateral interaction for the nearest neighbor interaction at low coverage (Θ_C = 0.33 ML) is slightly attractive. The excess energy (defined as the difference of the energy of the co-adsorbed system and the energy obtained at a coverage of 0.25 ML) of −50 meV at this coverage implies an attractive lateral interaction of −50 meV per interaction. The lateral interaction for the nearest neighbor interaction becomes repulsive at higher coverage of Θ_C = 0.5 ML. A lateral interaction of 66 meV per C-C interaction can be derived from excess energy of 131 meV in the p(2x2) cell. The repulsive interaction increases to 85 meV per C-C interaction upon increasing the coverage to Θ_C = 1 ML. The next-nearest neighbor interaction is attractive by −34 meV per C-C interaction as reported previously [24,27,29] resulting in the dominance of this configuration in agreement with experimental LEED-studies [31,32]. The lateral interaction varies with temperature favoring the next-nearest neighbor interaction even more over the nearest neighbor interaction.

The adsorption of atomic oxygen on Fe(100) has been studied well using DFT [27,30,33]. Here, we report a resulting Fe-O bond length at this coverage of 2.138 Å. The obtained adsorption energy of atomic oxygen at a coverage of 0.25 ML (−3.77 eV/O relative to gas phase O₂) is in reasonable accordance with the reported adsorption energies (in particular noting the difficulty in determining the energy of oxygen in the gas phase). The adsorption of atomic oxygen on Fe(100) at Θ_O = 0.25 ML is associated with some deformation of the p(2x2) unit cell, although the deformation energy is lower than upon the adsorption of atomic carbon at the same coverage. Increasing the coverage to 0.33 ML on a p(3x2) cell increases the deformation energy, but the overall lateral interaction remains slightly attractive. At a coverage of Θ_O = 0.5 ML, the nearest neighbor interaction is repulsive, whereas the next-nearest neighbor interaction is negligibly small. The temperature dependency of lateral interaction is negligible at this coverage.

Co-adsorption of CO and atomic carbon was probed on a p(3x2) and a p(2x2) unit cell. It does not seem to affect the Fe-C bond length if C and CO are adsorbed adjacent to each other, but this bond length is shortened if they are adsorbed diagonally from each other. The C-O bond length is shortened for the adsorption of CO and C adjacent to each other in a p(2x2) unit cell. Co-adsorption of CO and C is associated with quite a large deformation in particular of the surface layer resulting in relatively high deformation energy. The nearest neighbor interaction between C and CO is repulsive, whereas as reported before [6,24,27,30] the next-nearest neighbor interaction is attractive.

Similar trends were obtained in the co-adsorption of CO and atomic oxygen, albeit with less deformation of the surface layer resulting in lower deformation energy (more details in Table 2).

Co-adsorption of atomic carbon and atomic oxygen results in large deformation energy with a significantly shortened Fe-C bond length. As noted before, the nearest neighbor interaction between C and O is repulsive, whereas the next-nearest neighbor interaction is attractive. The temperature dependency of the lateral interaction between C and O is rather small and may be neglected.

2.3. Micro-Kinetic Model of CO-TPD

2.3.1. Setting-Up of Micro-Kinetic Models

The temperature programmed desorption of CO from Fe(100) was modeled using two mean-field micro-kinetic models, viz. Bragg–Williams approximation (BWA) and the quasi-chemical approach (QCA), and a kinetic Monte Carlo (kMC) method [34]. The mean-field approach presumes a desorption process as a result of a single averaged effect. In the Bragg–Williams approximation (BWA) the likelihood for a certain species at an adjacent site is given by the fractional coverage of that species. This kind of model may work well for systems with low or negligible mobility of the adsorbed species (either through diffusion or via desorption and re-adsorption) resulting in a ‘frozen’ surface structure since surface species cannot diffuse to an energetically more favorable site. In the quasi-chemical approach (QCA), the likelihood for a certain species at an adjacent site is dependent on the interaction energy between the species at the central site and the species at the neighboring sites. This type of models may work sufficiently well in systems with high mobility of surface species. The more complex kinetic Monte Carlo (kMC) model can capture subtleties in the catalytic process better, e.g., induced effects by finite rates of diffusion.

Essentially only three reaction steps need to be considered when describing temperature programmed desorption of CO taking place in a vacuum, i.e.,

$$\mathrm{CO-desorption}\hspace{1em}{\mathrm{CO}}^{*}\to \mathrm{CO}\left(\mathrm{g}\right)+ \ast \hspace{1em}{\mathrm{r}}_1={\mathrm{k}}_{\mathrm{des},\mathrm{CO}}·{\mathsf{\theta }}_{\mathrm{CO}}$$

(1)

$$\mathrm{CO-dissociation}\hspace{1em}{\mathrm{CO}}^{*}+ \ast \to { \mathrm{C}}^{*}+{ \mathrm{O}}^{*}\hspace{1em}{\mathrm{r}}_2={\mathrm{k}}_2·{\mathsf{\theta }}_{\mathrm{CO}}·{\mathsf{\theta }}_{*}$$

(2)

$$\mathrm{C}/\mathrm{O-recombination}\hspace{1em}{\mathrm{C}}^{*}+{ \mathrm{O}}^{*} \to { \mathrm{CO}}^{*}+ \ast \hspace{1em}{\mathrm{r}}_3={\mathrm{k}}_3·{\mathsf{\theta }}_{\mathrm{C}}·{\mathsf{\theta }}_{\mathrm{O}}$$

(3)

In addition to these reaction steps, diffusion of the various surface species is considered in the kMC-model:

$$\mathrm{Diffusion} \mathrm{of} \mathrm{CO}\hspace{1em}{\mathrm{CO}}^{*}+ \ast \to \ast +{ \mathrm{CO}}^{*}\hspace{1em}{\mathrm{r}}_4={\mathrm{k}}_{\mathrm{dif},\mathrm{CO}}·{\mathsf{\theta }}_{\mathrm{CO}}·{\mathsf{\theta }}_{*}$$

(4)

$$\mathrm{Diffusion} \mathrm{of} \mathrm{C}\hspace{1em}{\mathrm{C}}^{*}+ \ast \to \ast +{ \mathrm{C}}^{*}\hspace{1em}{\mathrm{r}}_5={\mathrm{k}}_{\mathrm{dif},\mathrm{C}}·{\mathsf{\theta }}_{\mathrm{C}}·{\mathsf{\theta }}_{*}$$

(5)

$$\mathrm{Diffusion} \mathrm{of} \mathrm{O}\hspace{1em}{\mathrm{O}}^{*}+ \ast \to \ast +{ \mathrm{O}}^{*}\hspace{1em}{\mathrm{r}}_6={\mathrm{k}}_{\mathrm{dif},\mathrm{O}}·{\mathsf{\theta }}_{\mathrm{O}}·{\mathsf{\theta }}_{*}$$

(6)

The rate constants for these steps are given by the transition state theory. A thermodynamic consistent expression for the rate constant for CO desorption is obtained from the equilibrium constant for adsorption of CO and the rate constant for the adsorption of CO. The latter is given by the impingement and sticking of CO on to a surface site with the area, A_site. Assuming a sticking probability of one for the adsorption of CO, the rate constant for desorption of CO is then given by:

$${\mathrm{k}}_{\mathrm{des},\mathrm{CO}}=\frac{{\mathrm{k}}_{\mathrm{ads},\mathrm{CO}}}{{\mathrm{K}}_{\mathrm{ads},\mathrm{CO}}}=\frac{{\mathrm{A}}_{\mathrm{site}}}{\sqrt{2\mathsf{\pi }·{\mathrm{m}}_{\mathrm{CO}}·{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}}·\frac{{\mathrm{Q}}_{\mathrm{CO}\left(\mathrm{g}\right)}}{{\mathrm{Q}}_{{\mathrm{CO}}^{*}}}·\exp \left(-\frac{{\mathrm{E}}_{\mathrm{des},\mathrm{CO}}\left(\mathsf{\theta }\right)}{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}\right)$$

(7)

where Q_CO(g) and Q_CO* are partition functions for CO in the gas phase and adsorbed state, respectively, taking into account the vibrational and for CO in the gas phase also the translational and rotational partition function and E_des,CO(θ) is the coverage dependent desorption energy for CO. The rate constants for the surface reactions and diffusion steps are given by

$${\mathrm{k}}_{\mathrm{i}}=\frac{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}{\mathrm{h}}\frac{{\mathrm{Q}}_{\mathrm{TS}}^{‡}}{{\mathrm{Q}}_{\mathrm{IS}}}\exp \left(-\frac{{\mathrm{E}}_{\mathrm{a}}\left(\mathsf{\theta }\right)}{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}\right)$$

(8)

with ${\mathrm{Q}}_{\mathrm{TS}}^{‡}$ and ${\mathrm{Q}}_{\mathrm{IS}}$ are partition functions for the transition state and initial state, respectively, only considering the vibrational contribution while h is Planck’s constant and E_a(θ) the coverage dependent activation energy. Hence, the pre-exponential factors in all our models are dependent on temperature. It should be noted that in principle lateral interactions also affect the pre-exponential factor due to a change in the vibrational partition function, but this can be neglected seeing the relatively small change in the partition function.

The coverage dependent reaction energies can be obtained from the energy without lateral interaction (here taken as the reaction energy obtained from DFT at a coverage of 0.25 ML—see

Table 3. Reaction energies, energy barriers and pre-exponential factors at 200 K in the absence of lateral interaction for the reaction steps in CO-TPD.

Reaction	ΔErxn, eV	Ea,f, eV	Ea,r, eV	Af(200 K), s−1	Af/Ar
$\mathrm{CO}\left(\mathrm{g}\right)\to {\mathrm{CO}}^{*}+\ast$	−1.88	0.0	1.88	2.89 × 108 1	2.50 × 10−10
${\mathrm{CO}}^{*}\to { \mathrm{C}}^{*}+{\mathrm{O}}^{*}$	−1.05 2	1.11 3	2.16	1.14 × 1013	1.20
${\mathrm{CO}}^{*}+ \ast \to \ast +{ \mathrm{CO}}^{*}$	0.00	0.75 4	0.75	6.41 × 1011	1.00
${\mathrm{C}}^{*}+ \ast \to \ast +{ \mathrm{C}}^{*}$	0.00	1.49 4	1.49	8.24 × 1012	1.00
${\mathrm{O}}^{*}+ \ast \to \ast +{ \mathrm{O}}^{*}$	0.00	0.57 4	0.57	7.32 × 1012	1.00

¹ Units bar⁻¹·s⁻¹; ² Reaction energy to create adsorbed C and adsorbed O at infinite distance; ³ as reported by [6]; ⁴ as determined from the energy for adsorption in the bridge position as reported by [6].

) and the interaction energy, ε_i-j (see

Table 4. Interaction energy, ε_i-j (meV) for nearest and next-nearest neighbor interactions.

Nearest neighbor interaction
	CO	C	O
CO	180	144	194
C	144	66	320
O	194	320	190
Next-nearest neighbor interaction
	CO	C	O
CO	−14	−100	−62
C	−100	−34	−95
O	−62	−95	0

). The interaction energy is dependent on the orientation of the species involved and their environment. This complex relationship may be simplified by considering only the interaction energy obtained at a coverage of 0.5 ML. This may be a suitable simplification for the description of the CO-desorption process if interactions are repulsive and diffusion of the surface species is sufficiently high resulting in a low probability of nearest neighbor interactions at low coverage. Under these assumptions, the interaction energy for the nearest neighbor interaction is obtained from the excess energy obtained in the co-adsorption of CO-CO, CO-C, CO-O, and C-O at 0.5 ML in the p(2x2) unit cell with the adsorbates in adjacent positions. The interaction energy for the next-nearest neighbor interaction is obtained from the excess energy in the co-adsorbed systems with the adsorbates diagonally orientated.

The interaction energies, ε_i-j, are in principle dependent on temperature. The temperature dependency of CO-CO and the C-O interaction may be neglected (vide supra). During CO-TPD, the system is expected to have an initially a high coverage with CO and relatively low coverage at high temperature. Hence, the temperature dependency is not incorporated.

The coverage dependent activation barriers can be obtained by invoking the Brønsted−Evans−Polanyi correlation and Hammond’s postulate [35], i.e., implying that the curvature of the potential energy surface does not change significantly due to lateral interactions [36]:

$${\mathrm{E}}_{\mathrm{a},\mathrm{f}}\left(\mathsf{\theta }\right)={\mathrm{E}}_{\mathrm{a},\mathrm{f}}^0+\mathsf{\alpha }·\left[{\mathsf{\Delta }\mathrm{E}}_{\mathrm{rxn}}^{\mathrm{X}}\left(\mathsf{\theta }\right)-{\mathsf{\Delta }\mathrm{E}}_{\mathrm{rxn}}^0\right]$$

(9)

with ${\mathrm{E}}_{\mathrm{a}}^0$ and ${\mathsf{\Delta }\mathrm{E}}_{\mathrm{rxn}}^0$ the activation barrier and the reaction energy respectively in the absence of lateral interactions, ${\mathsf{\Delta }\mathrm{E}}_{\mathrm{rxn}}^{\mathrm{X}}\left(\mathsf{\theta }\right)$ the reaction energy in the presence of lateral interactions. Hence, it is assumed that the change in the activation energy due to lateral interactions is proportional to the change in the reaction energy due to lateral interactions.

The factor α is in the limit case 0 for reactions with an early transition state, such as CO dissociation, and 1 for reactions with a late transition state, such as C/O recombination. The factor α is taken to be 0.5 for the diffusion of CO, C, and O.

The coverage dependent change in energy upon reaction is given by the difference in the energy of the product and reacting species involved in the reaction taking into account their interaction with neighboring species, and thus the difference in the cluster energy of the product and the reacting species. According to the Bragg–Williams approximation (BWA) the coverage dependent cluster energy for a system containing m species is determined from the energy of the species without lateral interactions (${\mathrm{E}}_{\mathrm{i}}^0$) and the interaction of a central adsorbate i with the nearest neighbor j (${\mathsf{\epsilon }}_{\mathrm{i}-\mathrm{j}}^{\mathrm{NN}}$) and its interaction with the next-nearest neighbor k (${\mathsf{\epsilon }}_{\mathrm{i}-\mathrm{k}}^{\mathrm{NNN}}$):

$${\mathrm{E}}_{\mathrm{i}}^{\mathrm{BWA}}={\mathrm{E}}_{\mathrm{i}}^0+{\mathrm{n}}_{\mathrm{NN}}·{\displaystyle \sum }_{j=1}^m{\mathsf{\theta }}_{\mathrm{j}}·{\mathsf{\epsilon }}_{i-j}^{NN}+{\mathrm{n}}_{\mathrm{NNN}}·{\displaystyle \sum }_{k=1}^m{\mathsf{\theta }}_{\mathrm{k}}·{\mathsf{\epsilon }}_{i-k}^{NNN}$$

(10)

with n_NN the number of nearest neighboring sites (n_NN = 4 on Fe(100)) and n_NNN the number of next-nearest neighboring sites (n_NNN = 4 on Fe(100)).

In the quasi-chemical approach (QCA) the coverage dependent cluster energy for a system containing m species is given by [37]:

$$\begin{array}{ll}{\mathrm{E}}_{\mathrm{i}}^{\mathrm{QCA}}={\mathrm{E}}_{\mathrm{i}}^0+{\mathrm{k}}_{\mathrm{B}}& ·\mathrm{T}\\ & ·({\mathrm{n}}_{\mathrm{NN}}·\ln \left[{\displaystyle {\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{m}}}\mathrm{P}\left({\mathsf{\theta }}_{\mathrm{i}}{|\mathsf{\theta }}_{\mathrm{j}}^{\mathrm{NN}}\right)·\exp \left(\frac{{\mathsf{\epsilon }}_{\mathrm{i}-\mathrm{j}}^{\mathrm{NN}}}{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}\right)\right]+{\mathrm{n}}_{\mathrm{NNN}}\\ & ·\ln [{\displaystyle {\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{m}}}\mathrm{P}\left({\mathsf{\theta }}_{\mathrm{i}}{|\mathsf{\theta }}_{\mathrm{k}}^{\mathrm{NNN}}\right)·\exp \left(\frac{{\mathsf{\epsilon }}_{\mathrm{i}-\mathrm{k}}^{\mathrm{NNN}}}{{\mathrm{k}}_{\mathrm{B}}·\mathrm{T}}\right)])\end{array}$$

(11)

where ${\mathsf{\theta }}_{\mathrm{i}}$ and $\mathrm{m}$, are fractional coverage of species-$\mathrm{i}$ and the total number of surface species, respectively. $\mathrm{P}({\mathsf{\theta }}_{\mathrm{i}}{|\mathsf{\theta }}_{\mathrm{j}})$ is the probability that the neighboring site is occupied by species-$\mathrm{j}$ given that the central site is occupied by species $\mathrm{i}$. This probability is dependent on the probability of finding a pair of sites occupied by species $\mathrm{i}$ and species $\mathrm{j}$. A detailed description to calculate these probability functions is given by Hellman and Honkala [37] based on the work by Zhdanov [3].

The determination of the CO-TPD-profile using mean-field models involves solving a set of three ODEs. In contrast, the kinetic Monte Carlo (kMC) algorithms identify all possible sets of events which may occur within a given system state, <σ>. An event, reaction step, is chosen depending on its kinetic rate constant and its multiplicity—faster reactions are chosen more frequently. The algorithm then evolves from a given initial state, <σ₀>. In this study, we used the Zacros, version 2.0, software package [38,39,40,41] which is based on the graph-theoretical kinetic Monte Carlo framework [39]. Detailed descriptions of this package and its implementation can be found elsewhere [40,41]. Briefly, the total lattice energy Hamiltonian, H(σ), is calculated using the cluster expansion approach [40]:

$$\mathrm{H}\left(\mathsf{\sigma }\right)={\displaystyle \sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{c}}}\frac{{\mathrm{ECI}}_{\mathrm{k}}}{{\mathrm{GM}}_{\mathrm{k}}}·{\mathrm{NCE}}_{\mathrm{k}}\left(\mathsf{\sigma }\right)$$

(12)

where $\mathsf{\sigma }$ describes the occupancy of each lattice site, N_c is the number of clusters considered to account for lateral interaction effects on the total lattice energy. The summation is done over all $k$ clusters where ECI_k, GM_k and NCE_k are the effective cluster energy, cluster multiplicity and the number of occurrences of cluster k within a given lattice state, <σ> [42]. The cluster energy, ECI_k, is obtained from the energy of the central adsorbate (see Table 3) taking lateral interactions with the adjacent species into account (see Table 4). The cluster multiplicity, GM_k, the number of occurrences of cluster k are determined for a given lattice state using a detection algorithm [40].

2.3.2. Describing CO-TPD Using Mean Field Models

Figure 2 shows simulated spectra of temperature programmed desorption of CO from Fe(100) obtained from different starting coverages of CO using a mean-field approach. In the absence of lateral interactions, the TPD-profile starting with an initially saturated surface is characterized by two peaks with a peak maximum at 500 K and 810 K, respectively. These peaks can be attributed to the desorption of molecularly adsorbed CO and the desorption of CO after re-combination of surface carbon and surface oxygen, which have been formed upon dissociation of adsorbed CO. Decreasing the initial coverage of the surface to Θ_CO,initial = 0.75 ML decreases the magnitude of the first peak, but not its position. It can be further noted that the maximum in the high temperature region remains unchanged in both its magnitude and its position. The first peak almost completely disappears upon decreasing the initial coverage further to Θ_CO,initial = 0.50 ML and the second peak remains again unaltered. This can be rationalized since the rate of CO-dissociation at high coverage is limited by the availability of vacant sites. At sufficiently low initial coverage, CO will dissociate before it desorbs as the activation barrier for CO dissociation is lower than the activation barrier for CO desorption. Dissociation of CO for experiments started with an initial coverage of Θ_CO,initial > 0.50 ML will result in a surface almost saturated with atomic carbon and atomic oxygen. Decreasing the initial coverage to Θ_CO,initial = 0.25 ML results in a small shift of the high temperature peak to ca. 830 K. Starting with a low initial coverage of CO results in only a partial coverage of the surface with carbon and oxygen, thus limiting the rate of recombination and thus the desorption of CO from the surface.

The simulated TPD-spectra of CO from Fe(100) using the Bragg–Williams approximation also shows in principle two maxima. Starting from Θ_CO,initial = 1.00 ML the on-set of the first peak is estimated to be at ca. 300 K showing the highest rate of desorption at 360 K. Decreasing the initial CO-coverage to Θ_CO,initial = 0.75 ML and to Θ_CO,initial = 0.50 ML increases the onset temperature for the first peak to 330 K and 370 K, respectively. The maximum rate of desorption is then obtained at 403 K and 416 K, respectively. The occurrence of CO-desorption at significantly lower temperatures predicted by the Bragg–Williams approximation in comparison to the simulated CO-TPD spectra in the absence of lateral interactions is attributed to the repulsive nearest-neighbor interactions in this system, and in particular to the repulsive nature of the nearest neighbor CO-CO interaction. The increased broadness of the peak is attributed to the gradual change in the strength of adsorption due to the change in the surface coverage. It can be further noted that the tail of the first peak coincides for the experiments with an initial coverage of CO of 1.00 ML, 0.75 ML, and 0.50 ML. At these temperatures, the coverage of the surface is independent of the initial coverage with CO, and hence the desorption profiles become indistinguishable from this temperature onwards.

In contrast to the simulated CO-TPD profile with Θ_CO,initial = 0.25 ML in the absence of lateral interactions, the same system simulated using the Bragg–Williams approximation does show a low temperature peak with a peak maximum at ca. 450 K. CO-dissociation taking place at these temperatures result in a higher surface coverage (with atomic carbon and atomic oxygen) thereby increasing lateral interactions thus favoring desorption of some CO from the surface.

The high temperature peak in the simulated CO-TPD profile using the Bragg–Williams approximation is strongly shifted towards lower temperatures in comparison to the simulated CO-TPD profile in the absence of lateral interactions. Strong repulsive nearest-neighbor interactions between in particular C and O will favor the recombination of C and O, a reaction with a late transition state, resulting in the formation and desorption of CO at relatively low temperatures. Reduction of the lateral interactions, through desorption of CO, will thus reduce the rate of C-O recombination and slow down the rate of CO-desorption. The second peak maximum in the simulated CO-TPD profile for the experiment with Θ_CO,initial = 0.25 ML is again shifted slightly to higher temperatures. This can be attributed to the reduced likelihood of nearest neighbor interactions in this system due to the lower surface coverage. This was confirmed by reducing the initial coverage to 0.05 ML (not shown), which results in peak maximum at high temperature close to the one obtained ignoring lateral interactions.

The simulated spectra of CO-TPD from Fe(100) starting from an initially saturated surface using the quasi-chemical approach (QCA) shows three maxima at 327 K, 440 K and a broad maximum around 820 K. The appearance of two separate peaks at low temperature rather than a broad peak as observed in the simulated CO-TPD using the Bragg–Williams approximation, can be attributed to the instantaneous rearrangement of the surface upon generation of vacancies on the surface to the lowest energy state minimizing strongly repulsive nearest neighbor interactions. Decreasing the initial coverage to Θ_CO,initial = 0.75 ML results in a strong decrease of the first peak (as well as a slight shift in the peak maximum to 338 K) as the number of nearest neighbor interactions would have been lower from the start of the experiment. Decreasing the initial coverage further to Θ_CO,initial = 0.50 ML results in the disappearance of the low temperature peak. Only a rather small amount of CO desorbs at a low temperature (440 K) if the initial coverage is lowered to Θ_CO,initial = 0.25 ML.

The observed new low temperature peak in addition to the 2^nd low temperature peak in the simulated CO-TPD using the quasi-chemical method starting with an initially saturated surface compare well with the appearance of a new peak with a maximum at 306 K (α₂-CO) in addition to a peak at 440 K (α₃-CO) in the CO-TPD from Fe(100) with increasing dosing of CO [18]. The tentative assignment of α₂-CO in CO-TPD from Fe(100) to CO bonded to a bridge site [19] may rather be related to lateral interactions in the system.

The high temperature peak in the simulated CO-TPD using the quasi-chemical method is broad, and its maximum is shifted significantly towards higher temperatures in comparison to the simulated TPD using the Bragg–Williams approximation. Strong repulsive interactions can be minimized in systems with sufficient vacant sites, thus minimizing the rate of desorption from the surface.

The input parameters in the kinetic models are associated with the uncertainty inherent to the method of determining these factors. Figure 3 shows the influence of the variation of the strength of the nearest neighbor interaction on the simulated CO-TPD-profile. Increasing the strength of the lateral interactions between CO and CO by 10% shifts the broad, low temperature peak obtained using the Bragg–Williams approach to lower temperatures, but the essential feature in this profile, i.e., the broad single peak at a low temperature remains. It can be further noted that the tail of the first peak is hardly affected. A change in the CO-CO lateral interaction shifts only the low temperature peak in the CO-TPD profile simulated using the quasi-chemical approach, but not the peak with a maximum at ca. 440 K. This implies that at the adsorption of CO, which desorbs at this temperature is no longer affected by lateral interactions between adsorbed CO and CO on adjacent sites.

The low temperature region of the CO-TPD simulated using the Bragg–Williams approximation is only marginally affected by a change in the lateral interactions between CO and C and the lateral interactions between CO and O. Increasing the lateral interactions forces the desorption of some more CO at lower temperatures starting to result in the formation of a new maximum in the rate of desorption at ca. 440 K. The recombination of C and O is affected by the lateral interaction as indicated by a shift in the peak maximum at high temperature region towards higher temperatures with increasing strength of CO-X lateral interactions. This is not seen in the CO-TPD simulated using the quasi-chemical approach as these lateral interactions are minimized when applying this framework.

The CO-TPD profiles were simulated assuming an early transition state for CO dissociation. A shift in the transition state from a state resembling more the initial state to a state resembling more the product state will affect the energy barrier for CO-dissociation through a change in the difference in the reaction energy incorporating lateral interactions and the reaction energy in the absence of lateral interactions (see Equation (9)). This difference in the reaction energy for CO-dissociation becomes higher if the extent of CO dissociation increases whilst keeping the fraction of surface which is vacant constant, and upon decreasing the fraction of the surface which is vacant keeping the ratio of CO adsorbed to C/O adsorbed constant. Hence, the activation energy for CO dissociation will become higher for a late transition state. Similar reasoning can be put forward for the reverse reaction taking into account that the difference in the reaction energy incorporating lateral interactions is less than the reaction energy in the absence of lateral interactions. Hence, the activation energy for C/O recombination is lower with an early transition state than with a late transition state.

The tail of the first peak in the CO-TPD profile simulated using the Bragg–Williams approximation shifts towards higher temperatures upon changing from an early to a late transition state (see Figure 4). This can be attributed to an increase in the activation energy for CO-dissociation upon increasing the α value. The increase in the activation energy for CO-dissociation also results in a decrease of the surface coverage with C and O, e.g., at temperatures around 500 K. A higher coverage with atomic C and O on the surface will according to the Bragg–Williams approximation result in higher repulsive lateral interactions between the adsorbed species and thus favoring the C/O recombination. Hence, the high temperature peak (β) shifts towards lower temperatures upon decreasing α.

The observed variations in the CO-TPD profiles simulated using the quasi-chemical approach upon changing the transition state from an early transition state to a late transition state are much subtler than in the same variation in the CO-TPD profile simulated using the Bragg–Williams approximation. The first peak is not affected by changing the assumption on the transition state for CO-dissociation. Interestingly, the second peak with a maximum at ca. 440 K becomes broader and start showing features of two separate events assuming a late transition state. Increasing the value of α, whilst keeping the fraction of vacant sites approximately constant, will increase the reaction barrier for CO dissociation. As a consequence, the fractional coverage of the surface with CO at 400 K is slightly higher if CO dissociation has a late transition state, which thus yields a slightly higher rate of desorption of CO at this temperature. The delay in the build-up of appreciable amounts of carbon and oxygen on the surface further implies that the influence of those lateral interactions on the low temperature peak will be shifted towards slightly higher temperatures.

The high temperature peak in the CO-TPD profile becomes sharper upon changing from an early to a late transition state. Increasing the value of α will result in a shift from an early transition state to a later transition state increasing the energy barrier for C/O recombination.

2.3.3. Describing CO-TPD Using Kinetic Monte-Carlo (kMC)

The applied kinetic Monte Carlo explicitly considers diffusion of the species on the surface and as a consequence the local structure. The simulated TPD-profile of CO from Fe(100) using kinetic Monte Carlo (see Figure 5) shows similarity to the profile generated using a mean-field quasi-chemical approach (cf. Figure 3) with three main desorption maxima. However, there are some differences between the two CO-TPD profiles. The maximum of the first peak is shifted towards slightly higher temperatures at high coverage (from 320 K to 330 K). This might be attributed to the local local stabilization accounted for in the kinetic Monte Carlo simulation through a diffusional shift of adsorbed CO. The second peak maximum is shifted to a significantly lower desorption temperature (from 430 K to 390 K). At this temperature, CO starts to dissociate. The slow diffusion of C/O from adjacent CO will increase the rate of desorption of CO.

A particular feature of kinetic Monte Carlo is the inclusion of local effects rather than global averaging effect. It can be noted that upon dissociation (and given sufficient time), the system seems to prefer C-C and O-O neighboring couples rather than C-O as nearest neighbors (as could be deduced the data given in Table 4). Thus, CO desorption at temperatures higher than 700 K thus requires first diffusion of the reactants to adjacent positions. This slows down the desorption of dissociated CO and the maximum in the rate of desorption is seen to be shifted towards higher temperatures in the kinetic Monte Carlo simulation compared to the simulation using the mean-field, quasi-chemical approach.

Figure 6 shows the sensitivity of the simulated CO-TPD profile using the kinetic Monte Carlo method on the choice of the various parameters. The absence of lateral interactions results in a CO-TPD profile with only two maxima. This profile is similar to the one calculated using a mean-field approach, although the high temperature peak (β-CO) is shifted to a slightly lower temperature when using kMC. This can be attributed to a dilution effect. In the mean-field approach, the average fractional coverage of C and O is taken into account to estimate the rate for the C/O recombination, whereas kMC takes local coverages into account. The incorporation of nearest neighbor lateral interactions results in the low temperature region of the CO-TPD profile in a low temperature peak, as a result of CO desorption experiencing strong lateral interaction with co-adsorbed CO, and a peak with a maximum at ca. 405 K, due to the desorption of more strongly bound CO. The low temperature peak in the CO-TPD profile shifts towards higher temperatures upon incorporation of the slightly attractive, next-nearest neighbor interactions. The high temperature peak (β-CO) is also affected by the inclusion of lateral interactions. The inclusion of the repulsive, nearest-neighbor interactions results in the CO-desorption already at significantly lower temperatures: the repulsive interaction forces the C/O recombination. This is countered by the attractive next-nearest neighbor interactions.

Increasing the lateral interaction by 10% results in a shift in the low temperature peak in the CO-TPD with 19 K towards lower temperature. Similarly, a decrease in energy associated with CO-CO interaction causes the low temperature peak to shift to a higher temperature by ca. 19 K. A similar effect was observed in the sensitivity analysis using the quasi-chemical approach (QCA), but the position of the peak maximum was ca. 10 K higher in the kinetic Monte Carlo simulation in comparison to the simulation using the quasi-chemical approach (QCA). This is a result of the finite rate of diffusion of the species in the Monte Carlo simulation compared to the instantaneous equilibration of the lateral interactions on the surface assumed in the QCA-simulation.

The CO-TPD profile simulated using kinetic Monte Carlo was virtually unaffected by the strength of the lateral interactions between CO and C and CO and O (as also noted in the CO-TPD profile simulated using the quasi-chemical approach). The diffusion of the various species is fast enough at the temperatures where these interactions start to play a role to ensure that the repulsive interactions can be minimized.

The energy barrier for the diffusion over the surface of the various species does not affect the first peak in the TPD profile of CO from Fe(100). Also, the second peak in the CO-TPD profile is not affected by increasing the rate of diffusion (i.e., by decreasing the energy barrier for diffusion). However, slowing down the diffusion by increasing the energy barrier for the diffusion with a factor of 2 results in a significant broadening of the second peak with some evidence for peak splitting. CO on a site is then more likely to have a species on an adjacent site (e.g., formed upon dissociation of another CO), since it will diffuse away rather slowly. Thus, the increased coverage surrounding adsorbed CO may induce further CO-desorption. Furthermore, co-adsorbed C and O formed in adjacent sites upon CO-dissociation will not separate easily, which may result in recombination and subsequent desorption. The high temperature peak in the CO-TPD profile (β-CO) hardly affected by increasing the energy barrier for diffusion by a factor of two. Interestingly, the high temperature peak is affected by lowering the diffusion barrier. This may be attributed to the difference in the mobility of C and O, with surface carbon being much less mobile than surface oxygen (see Table 3). Decreasing the energy barrier for surface diffusion implies that both species will have similar, high mobility. They will try to avoid each other to minimize lateral interactions, thus resulting in a shift of the high temperature peak further upon increasing the rate of diffusion.

Shifting the assumption of an early transition state for CO dissociation to a late transition state affects mainly the second peak in the CO-TPD profile. A late transition state implies higher activation energy for CO dissociation, whilst at the same time, CO on adjacent sites will diffuse to energetically more favorable sites. This will lead to an increase in the desorption energy shifting the maximum for the second peak. It can be further noted that the second peak is significantly broadened since the stability of adsorbed CO depends on its immediate neighborhood (with the absence of nearest neighbors resulting in the strongest CO adsorption).

The experimentally observed CO-TPD profile from Fe(100) [18] seems to be modeled best using either the quasi-chemical approach (QCA) or the kinetic Monte Carlo method (kMC). The quasi-chemical approach is a mean-field model assuming a rapid equilibration between the interacting species. The kinetic Monte Carlo method accounts for both local effects and the finite rate of diffusion. Reducing the activation energy for the diffusion of CO did not change the low temperature region of the simulated CO-TPD profile using kMC significantly implying that the origin of the low temperature desorption maximum in the kMC model is due the local effect (i.e., the presence of CO still surrounded by co-adsorbed CO at the nearest neighbor position, which has relatively low activation energy for desorption). The local effect simulated in the kMC model is further demonstrated by slowing down the diffusion of in particular surface oxygen (vide supra). The similarity between the CO-TPD profile from Fe(100) simulated using the quasi-chemical approach (QCA), and the kinetic Monte Carlo method (kMC) may seem, therefore, fortuitous. However, the energy as modeled using the quasi-chemical approach (QCA—see Equation (11)) is in principle a sigmoid function rather than a linear function as predicted by the Bragg–Williams approximation. Hence, the increase in the strength of adsorption with decreasing coverage will be initially less in the QCA model compared to the BWA-model favoring desorption of less strongly adsorbed species. The strength of adsorption as predicted using increases rapidly upon decreasing the coverage further bringing the desorption process virtually to a halt. Thus, QCA starts mimicking the behavior observed with kMC for this system with repulsive lateral interactions and a low rate of surface diffusion. QCA appears to be the next best approach if the required parameter set for kMC is not available or if the kMC-model becomes computationally too cumbersome.

3. Methodology

Spin-polarized periodic density functional theory (DFT) calculations were conducted using the Cambridge Sequential Total Energy Package (CASTEP) [43] implemented in the Materials Studio^® software package. We used the generalized gradient approximation (GGA) exchange-correlation functional and the revised Perdew-Burke-Ernzerhof (rPBE) [44]. The ion-electron interactions were approximated using ultra-soft pseudopotentials (USSP) [45] with core corrections. The cutoff energy was set at 400 eV. For electronic structure minimization, Gaussian smearing was utilized with a smearing width, $k_B·T$, of 0.1 eV.

For all our DFT calculations we made use of five-layered Fe(100) slabs with a vacuum layer of at least 10 Å. The bottom three layers were fixed in their bulk optimized positions whilst all other atoms, including adsorbates, were allowed to relax until the residual forces on them were below 0.02 eV/Å. Adsorption was performed on only one side of the slab. This induced a dipole due to charge rearrangement on the surface which was corrected using the SCF method by Neugebauer and Scheffler [46]. We used various supercell sizes to estimate the interaction energy between two co-adsorbed reaction intermediates per interaction. The k-point grid was generated using the Monkhorst-Pack procedure [47] with a k-point spacing of <0.03 Å⁻¹. A partial Hessian analysis was conducted for the adsorbates [48,49]. This is a valid approximation as the Fe atoms are significantly heavier than those of C, H and O. The atoms were perturbed by 0.005 Å in Cartesian space.

The strength of adsorption of CO was investigated in different configurations on a p(4x4) cell and was determined from the difference in the electronic energy of the system with adsorbed CO and CO in the gas phase and a bare surface. The value was normalized with respect to the number of CO molecules adsorbed per unit cell

$${\mathrm{E}}_{\mathrm{ads}}=\left({\mathrm{E}}_{\mathrm{p}\left(4\mathrm{x}4\right)-\mathrm{Fe}\left(100\right)+\mathrm{nCO}}-{\mathrm{E}}_{\mathrm{p}\left(4\mathrm{x}4\right)-\mathrm{Fe}\left(100\right)}-\mathrm{n}·{\mathrm{E}}_{\mathrm{CO}\left(\mathrm{g}\right)}\right)/\mathrm{n}$$

(13)

Adsorption also results in a change in the position of the metal atoms in the slab. The deformation energy was determined by removing the adsorbates from the optimized geometry and determine the electronic energy of this system relative to the electronic energy of the bare Fe(100) surface:

$${\mathrm{E}}_{\mathrm{deformation}}={\mathrm{E}}_{(\mathrm{Fe}(100)+{\displaystyle \sum }_i{\mathrm{adsorbate}}_i)-{\displaystyle \sum }_i{\mathrm{adsorbate}}_i}-{\mathrm{E}}_{\mathrm{Fe}(100)}$$

(14)

The deformation energy is dependent on cell size. Lateral interaction can be described in terms of excess energy, which was evaluated relative to individually adsorbed species at Θ = 0.25 ML:

$${\mathrm{E}}_{\mathrm{excess}}=({\mathrm{E}}_{\mathrm{Fe}(100)+{\displaystyle \sum }_i{\mathrm{adsorbate}}_i}-{\mathrm{E}}_{\mathrm{Fe}(100)})-{{\displaystyle \sum }}_i\left({\mathrm{E}}_{\mathrm{p}(2x2)-Fe(100)+{\mathrm{adsorbate}}_i(\Theta =0.25ML)}-{\mathrm{E}}_{\mathrm{p}(2x2)-Fe(100)}\right)$$

(15)

4. Simulation Setup

We performed various simulations using the three microkinetic approaches to study CO-TPD over an Fe(100) surface. In a typical simulation, we ramped the temperature linearly from 200 K to over 1000 K at 10 K/s rate. Since the BWA and QCA models only involve the integration of a set of ODEs, a single run for a given condition is sufficient—we probed various starting CO coverages. We used Python’s SciPy ode package more specifically its LSODA ODE solver [50] which automatically detects stiffness and switches between the non-stiff Adams’ method and stiff backward differential formula (BDF) method [51]. The Jacobian is approximated using the finite difference method. In all QCA/BW simulations, the absolute and relative tolerances were set to 10⁻¹².

For the kMC model, we performed five repeat simulations each with a unique random number of generation seeds and reported average TPD profile. We used a 100 × 100 lattice cell, with 10,000 sites, seeded with a monolayer (except where stated otherwise) of CO* molecules—for a coverage less than 1 ML, the seeding was done randomly, and we used the same seeding for all simulations at a given initial coverage.

5. Conclusions

Lateral interactions show a variety and complexity depending on the considered coverage, which may be masked by countering effects, such as surface deformation. Kinetic models incorporating these lateral interactions are thus necessarily still simplified models – in this study, only average lateral interactions determined at a coverage of 0.5 ML were incorporated and the influence of the strength of lateral interaction on the simulated CO-TPD profile was probed using a sensitivity analysis.

The resulting models using either quasi-chemical approach (QCA) or kinetic Monte Carlo (kMC) result in the appearance of a separate low temperature desorption peak. The desorption peak at low temperature in the kinetic Monte Carlo method is caused by local effects, i.e., the presence of significant amounts of CO present on Fe(100) with multiple CO on next nearest neighboring sites, since surface diffusion does not appear to affect the appearance of this desorption peak. The appearance of the low temperature desorption peak in the simulated CO-TPD profile from Fe(100) is related to the sigmoid function describing the adsorption energy as a function of coverage.