A computational study of CO₂ hydrogenation on single atoms of Pt, Pd, Ni and Rh on In₂O₃(111)

Abstract

Metal promoted indium oxide (In₂O₃) catalysts are promising materials for CO₂ hydrogenation to products such as methanol and carbon monoxide. The influence of the dispersion of the promoting metal on the methanol selectivity of In₂O₃ catalysts is a matter of debate, which centers around the role of atomically dispersed single metal atoms vs. metal clusters as catalysts for methanol formation. In this study, we used density functional theory calculations to compare the role of single atoms (SAs) of Ni, Pd, Pt and Rh placed on the In₂O₃(111) surface to study CO₂ hydrogenation to CO and methanol. Direct and hydrogen-assisted CO₂ dissociation pathways leading to CO as well as methanol formation via either formate or CO intermediates are explicitly considered. Microkinetic simulations show that all SA models mainly catalyze CO formation via a redox pathway involving oxygen vacancies where adsorbed CO₂ dissociates followed by CO desorption and water formation. The higher barriers for hydrogenation of formate intermediates compared to the overall barrier for the rWGS reaction explain the negligible CH₃OH selectivity.

---

Introduction

The reduction of CO₂ emissions has become a crucial technological challenge because of the negative impact on the climate.^1–3 To address this challenge, the conversion of CO₂ captured from combustion processes and, in the future, directly from the air into liquid fuels and chemicals presents a viable solution, thereby reducing emissions and closing the carbon cycle. One of the most promising pathways for obtaining valuable products from CO₂ is via its hydrogenation, provided that the required hydrogen is obtained from renewable resources.^4–6 Not only does the large-scale hydrogenation of CO₂ into fuels offer an opportunity to decrease our reliance on fossil fuels for energy, this technology can also be leveraged to obtain chemical intermediates.^7–10 In this context, methanol is particularly attractive, because it can be used as a fuel or chemical building block.^11–13

Currently, large-scale methanol production is achieved by converting synthesis gas (CO/CO₂/H₂) over a Cu/ZnO/Al₂O₃ catalyst at temperatures of 473–573 K and pressures of 50–100 bar.¹⁴ Challenges arise however when CO₂ is hydrogenated instead of CO. The Cu-based catalyst shows significant activity for the reverse water-gas-shift reaction (rWGS), leading to the formation of unwanted CO.¹⁵ Moreover, the catalyst suffers from sintering, resulting in slow deactivation.^16,17 Many efforts have been made to identify new catalyst formulations with improved performance for the hydrogenation of CO₂ to CH₃OH.¹⁸ Recently, indium oxide (In₂O₃) has emerged as a promising catalyst for CO₂ hydrogenation to methanol.^19,20 Several density functional theory (DFT) studies emphasize the role of oxygen vacancies in the mechanism of CO₂ hydrogenation to methanol on In₂O₃.^19,21–25 Although such catalysts enable high selectivity to methanol by suppressing the competitive rWGS reaction, CO₂ conversion is limited by the intrinsically low activity in H₂ activation. To improve this, the use of metal promoters has been investigated. Recently, we reported by microkinetic modelling that small Ni clusters on In₂O₃ mainly catalyse methanol formation, whereas single atoms of Ni either doped in or adsorbed on In₂O₃ mainly catalyse CO.^26,27 In a recent study, Pinheiro Araújo et al. reported that, atomically dispersed metal atoms (Pd, Pt, Rh, Ni, Co, Au, Ir) on the In₂O₃ surface promote hydrogen activation and methanol production while hindering CO formation.²⁸ Frei et al. suggested by DFT calculations that single Ni atoms on In₂O₃ could catalyse the rWGS reaction but not methanol synthesis.²⁹ The same authors proposed that single atoms of Pd doped in In₂O₃ can stabilize Pd clusters on the In₂O₃ surface, which can enhance H₂ activation and, therefore, CH₃OH productivity.³⁰ Han et al. suggested that single atom (SA) of Pt doped in In₂O₃ would result in high methanol selectivity, whereas Pt nanoparticles on In₂O₃ would improve the rWGS reaction.³¹ Sun et al. showed by DFT calculations that the energetically most favourable pathway for methanol synthesis on a Pt₄/In₂O₃ model occurs via hydrogenation of a CO intermediate.³² Wang et al. reported that highly dispersed Rh phases on In₂O₃ can enhance dissociative H₂ adsorption and oxygen vacancy formation.³³ Dostagir et al. found that Rh doped into In₂O₃ catalyzes CO₂ hydrogenation to formate species.³⁴ Furthermore, earlier DFT works pointed out that SAs on the In₂O₃ surface do not catalyze the formation of CH₄.^{28–30,35,36} In nearly all these studies, the formation of atomically dispersed metal atoms was effective in shifting the selectivity towards the desired product. Although the role of supported clusters in metal-promoted In₂O₃ catalysts has been clarified for many cases,^32,37 generic understanding about the role of SAs is still lacking. Gaining a comprehensive understanding of the reaction pathways involved in the hydrogenation of CO₂ to synthesize methanol requires a thorough investigation of the role of SAs in In₂O₃.

In the present work, we study the reaction mechanism of CO₂ hydrogenation to CH₃OH and CO on SA-In₂O₃ model catalysts. Considering relevant literature, we have identified Ni, Pd, Pt, and Rh as potential promoters.^28–31,34 The reaction pathways explored in our DFT studies include a direct route for CO₂ hydrogenation to methanol (formate pathway), a pathway to methanol via CO hydrogenation and the competing rWGS reaction (both direct and hydrogen-assisted). We performed DFT calculations to determine the reaction energetics for all elementary reaction steps, which serve as input for microkinetic modelling. Microkinetic simulations show that the SA-In₂O₃ model surfaces produce CO as the main product with a very low methanol selectivity. The observed outcome can be primarily attributed to the presence of oxygen vacancies, which enable the adsorption of CO₂ and prompt activation towards facile direct C–O bond scission. In contrast, the formation of methanol is hampered by to the high activation energies associated with the hydrogenation of formate intermediates, resulting in relatively low rates for the preferred methanol product.

Computational methods

Density functional theory calculations

The calculations reported in this work were done in the same manner as reported before.^26,27 All DFT calculations were conducted using the projector augmented wave (PAW) method³⁸ and the Perdew–Burke–Ernzerhof (PBE)³⁹ exchange-correlation functional as implemented in the Vienna ab initio simulation package (VASP) code.^40,41 The valence 5s and 5p states of In were treated explicitly as valence states within the scalar-relativistic PAW approach. By using a plane-wave basis set with a cut-off energy of 400 eV solutions to the Kohn-Sham equations were calculated. All calculations were spin-polarized. The Brillouin zone was sampled using a 3 × 3 × 1 Monkhorst-Pack grid. By using Gaussian smearing with a smearing width (σ) of 0.1 eV electron smearing was employed. The stoichiometric In₂O₃(111) surface was modelled as a 2D slab with periodic boundary conditions. To avoid the spurious interaction with neighbouring super cells, a 15.0 Å vacuum region was introduced in the c-direction. It was verified that the electron density approached zero at the edges of the periodic super cell in the c-direction. The bottom two layers were frozen, while the top two layers were allowed to perturb in all calculations. The supercell has dimensions of 14.57 Å × 14.57 Å × 26.01 Å. The In₂O₃(111) slab consisted of 96 O atoms and 64 In atoms, distributed in four atomic layers on top of which the SA species were placed.

The stable states of the elementary reaction steps pertaining to CO₂ hydrogenation were calculated using the conjugate-gradient algorithm. Eqn (1) defines the adsorption energies of adsorbates (ΔE^ads_X):

ΔE^ads_X = E_X+In2O3(111) − E_In2O3(111) − E_X. (1)

where E_X+In2O3(111) is the electronic energy of the catalyst plus adsorbate system, E_In2O3(111) is the reference energy of the In₂O₃(111) slab and E_X is the DFT-calculated energy of the adsorbate in the gas phase.

The influence of oxygen vacancies on the reaction energetics was investigated by removing oxygen atoms from the In₂O₃(111) lattice. Since oxygen vacancies in chemical reactors are formed via surface reduction by H₂, the energy required to remove an oxygen from In₂O₃ (ΔE_Ov) was calculated using H₂O as reference, according to the following equation

ΔE_Ov = E_{defective slab} − E_{stoichiometric slab} + E_H2O − E_H2, (2)

where E_{defective slab} is the electronic energy of the catalyst containing one oxygen vacancy, E_{stoichiometric slab} is the reference energy of the stoichiometric slab. E_H2O and E_H2 are the DFT-calculated energies of gas-phase H₂O and H₂, respectively. Herein, we include the electronic energy, the zero-point energy correction and a finite temperature correction of translational and rotational energy of each gas-phase molecule.

Transition states were determined using the climbing-image nudged elastic band (CI-NEB) method.⁴² A frequency analysis was performed to all states to verify that that stable states have no imaginary frequencies and transition states have a single imaginary frequency aligned with the reaction coordinate.⁴³ The Hessian matrix in this frequency analysis was constructed using a finite difference approach with a step size of 0.015 Å for displacement of individual atoms along each Cartesian coordinate. The normal mode vibrations found with these calculations were used to evaluate the zero-point energy (ZPE) correction and the vibrational partition functions.

Partial density of state (pDOS) analysis are conducted to analyze the electronic structure of each SA-In₂O₃ model using the Lobster package.^44,45 The atomic charges were calculated using the Bader charge method.⁴⁶

Microkinetic simulations

Microkinetic simulations were carried out based on the DFT-calculated activation and reaction energies. The kinetic network was modelled using a set of ordinary differential equations involving rate constants, surface coverages and partial pressures of gas-phase species. To integrate the differential equations with respect to time the linear multistep backwards differential formula method was employed with a relative and absolute tolerance of 10⁻⁸.^47–49

The net rate of the adsorption process of a gas-phase species i was calculated as:

r_i = k_i,adsθ*P_i − k_i,desθ_i (3)

where θ* and θ_i are the fraction of free sites and the fraction of coverage species i, respectively. k_i,ads/des is the rate constant for the adsorption/desorption process and P_i is the partial pressure of species i.

To model adsorption processes, we assumed that the adsorbate loses one translational degree of freedom in the transition state with respect to the initial state. As a result of this, the rate of adsorption derived from transition state theory is:

where A_st and m_i are the effective area of an adsorption site and the molar mass of the gas species, respectively. P and T are the total pressure and temperature, respectively, and k_B is the Boltzmann constant. To calculate the gas-phase entropy of the adsorbates we employed the thermochemical Shomate equation as given bywhere S⁰ is the standard molar entropy.⁵⁰ The parameters A–G from eqn (5) were obtained from the NIST Chemistry Webbook.⁵¹ To model the desorption processes, we assumed that the desorbing species acquires two translational degrees of freedom and three rotational degrees of freedom in the transition state compared to the initial state. The resulting rate of desorption derived from transition state theory is:Herein, k_des is the rate constant for the desorption of the adsorbate, h is the Planck constant, σ is the symmetry number and is equal to 1, θ_rot the rotational temperature, and ΔE_ads the desorption energy. The value of A_st is equal to 9 × 10⁻¹⁹ m⁻².

The rate constant (k) of an elementary reaction step is given by the Eyring equation:

where Q^≠ and Q are the partition functions of the transition state and its corresponding initial state, respectively, and ΔE_act is the ZPE-corrected activation energy.

We employed the concepts of the degree of rate control (DRC) developed by Kozuch and Shaik^52,53 and popularized by Campbell⁵⁴ as well as the degree of selectivity control (DSC)^54–56 to identify the steps that control the CO₂ consumption rate and the product distribution.

Herein, the degree of rate control coefficient is defined as

A positive DRC coefficient indicates that the elementary reaction step is rate-controlling, whereas a negative coefficient corresponds to a rate-inhibiting elementary step.

Results and discussion

Structure of SA-In₂O₃ models

To study the catalytic properties of single atoms (SAs) of Ni, Pd, Pt and Rh, we constructed four models in which the SA was placed on the In₂O₃(111) surface, which is the most stable termination of In₂O₃.⁵⁷ We first determined the most stable location of the SA on the In₂O₃(111) surface. For this, we calculated the adsorption energies of the SA in different positions (Table S1†). The most stable adsorption sites are depicted in Fig. 1. Single Ni, Pd and Pt atoms preferentially adsorb in bridge sites between two oxygen atoms, in agreement with previous theoretical studies (Fig. 1a–c).^30,58 The Rh SA prefers to coordinate to three oxygens (Fig. 1d). A partial density of states (pDOS) analysis of the SA-In₂O₃ models shows that the bonding is mainly through interactions of the SA metal d orbitals and O 2p orbitals (Fig. S2†). The supported Ni, Pd, Pt and Rh SAs possess positive charges of +0.56, +0.35, +0.30 and +1.03, respectively. The energy needed to remove the SA atom to infinite distance from the surface is reported in Table S1.† Ni is most strongly bound to the surface (E_ads = −512 kJ mol⁻¹), followed by Pt (E_ads = −481 kJ mol⁻¹), Rh (E_ads = −425 kJ mol⁻¹) and Pd (E_ads = −339 kJ mol⁻¹). To understand the mobility of SAs on the In₂O₃(111) surface, we determined the barrier for migration between stable adsorption sites (Table S2†). Surface migration of SAs is associated with relatively high barriers for Ni, Ph and Rh (E_a > 100 kJ mol⁻¹), whereas much lower barriers are found for Pd (E_a < 50 kJ mol⁻¹). This is in line with the weaker binding energy of Pd with the In₂O₃(111) surface compared to Pt, Ni and Rh.

Fig. 1 Structural models used in DFT calculations composed of single atoms on the In₂O₃(111) surface. (a) Ni-In₂O₃, (b) Pd-In₂O₃, (c) Pt-In₂O₃ and (d) Rh-In₂O₃. (green: Ni, light blue: Pd, purple: Pt; brown: Rh; red: O; grey: In).

Oxygen vacancy formation

Oxygen vacancies (Ov) in the In₂O₃ surface have been identified as the active sites for CO₂ hydrogenation.^19,21,22 Here, we determine the effect of surface SAs on the formation of such vacancies. To this purpose, we calculate the energy required to remove an oxygen atom on each SA-In₂O₃ model and compare the results with data for bare In₂O₃ (Fig. 2, Table S3 and Fig. S3†). Since oxygen vacancies during CO₂ hydrogenation are formed via surface reduction by H₂, the computed energies are referenced to gaseous H₂O. For the bare In₂O₃(111) surface (Fig. 2a), the energy needed to remove an oxygen atom (E_Ov) ranges between −61 and 8 kJ mol⁻¹, in line with previous theoretical calculations.⁵⁹ The presence of a Ni atom on In₂O₃ (Fig. 2b) significantly increases the oxygen vacancy formation energy, resulting in E_Ov values between 139 and 190 kJ mol⁻¹. The effect of a Pd SA is much less pronounced, resulting in values for E_Ov comparable to those for the bare In₂O₃(111) surface, namely between −77 and 80 kJ mol⁻¹. For the Pt-In₂O₃ case (Fig. 2d), E_Ov assumes higher values than for Pd ranging between −35 and 217 kJ mol⁻¹. The Rh-In₂O₃ model (Fig. 2e) has intermediate values for E_Ov ranging between 39 and 170 kJ mol⁻¹. Overall, the removal of an oxygen atom on the bare In₂O₃(111) surface is less energetically demanding in comparison to SA-In₂O₃ models, except for Pd, where the oxygen formation energy is very similar.

Fig. 2 Oxygen vacancy formation energies (E_Ov, kJ mol⁻¹) for the 12 surface oxygens on (a) bare In₂O₃, (b) Ni-In₂O₃, (c) Pd-In₂O₃, (d) Pt-In₂O₃, and (e) Rh-In₂O₃. The coloring of the 12 surface oxygen atoms corresponds to their E_Ov indicated in the color bar. Vacancy formation energies are given with respect to H₂O formation. All other atoms are colored in grey. The SAs are highlighted inside dashed circles.

Further inspection of Fig. 2b–d shows that, for each of the SA-promoted surfaces, formation of vacancies is easier for lattice oxygen atoms directly bonding to the SA than for O atoms bonding only to In atoms. Comparing the H₂O-referenced E_Ov values for O binding to the SA, we computed as lowest lowest E_Ov values of 141, −71, −35 and 39 kJ mol⁻¹ for Ni-, Pd-, Pd- and Rh-In₂O₃, respectivley (Table S4†). On bare In₂O₃, the lowest E_Ov value is −61 kJ mol⁻¹. The relatively more endothermic E_Ov values for Ni–, Pt– and Rh–In₂O₃ are in line with the strong adsorption energy of the SA on the surface (Table S1†). Because Pd binds to the surface less strongly, E_Ov is relatively speaking more exothermic than for the other SAs. To understand whether oxygen vacancies will be formed on these models under reaction conditions and what role they play in CO₂ hydrogenation, we study the mechanism of oxygen vacancy formation via surface reduction by H₂ as well as the CO₂ hydrogenation reaction pathways by means of microkinetic simulations (vide infra).

Reaction mechanism

We performed DFT calculations to determine the reaction mechanism of CO₂ hydrogenation to methanol (CH₃OH), carbon monoxide (CO) and water (H₂O). Hydrogenation of CO₂ to methane (CH₄) was not included, because methane formation was not observed in experiments with low-loaded metal-promoted In₂O₃ catalysts.^28,32,34 Furthermore, earlier DFT works pointed out that SAs on the In₂O₃ surface do not catalyse the formation of CH₄.^28–30,35 The mechanism is depicted in Fig. 3 in the form of a reaction network diagram. Following previous computational studies, we investigated the pathways for (i) the formation of oxygen vacancies through H₂O formation, (ii) the formation of CH₃OH via formate, and (iii) the reverse water-gas shift (rWGS) pathway leading to CO₂.^{32,57,60–62} The latter rWGS reaction can take place via either direct C–O bond cleavage in adsorbed CO or an H-assisted pathway involving COOH intermediate. In addition to methanol formation via formate, we considered the direct hydrogenation pathway of CO intermediate into CH₃OH. We included a comparison with the DFT results of Frei et al. on In₂O₃(111).⁵⁷Fig. 4 depicts the potential energy diagrams (PED) for each of these routes.

Fig. 3 Complete reaction networks of CO₂ hydrogenation to CO and CH₃OH for (a) Ni-, (b) Pd-, (c) Pt- and (d) Rh-In₂O₃. Each elementary step is numbered. The forward and backward reaction energies, as well as the adsorption and desorption energies are given in kJ mol⁻¹.

Next, we discuss the elementary reaction steps on each SA-In₂O₃ surface and highlight differences in activation energies and transition-state structures. The activation barriers are given with respect to the most stable adsorbed state for each intermediate.

The geometries corresponding to initial, transition and final states can be found in the ESI† (Tables S6–S9). The computed forward and backward barriers for all elementary reaction steps are given in Tables S10–S13.†

Oxygen vacancy formation

We first discuss the formation of oxygen vacancies (steps 1–3). On all SA-In₂O₃ surfaces, the oxygen corresponding to the lowest E_Ov value is removed. For Ni-, Pd and Pt-In₂O₃, this corresponds to forming an oxygen vacancy in position Ov₃ (Table S2 and Fig. S3†). For Rh–In₂O₃ formation of an oxygen vacancy in position Ov₂ is preferred. On all SA-In₂O₃ models, H₂ is heterolytically dissociated into SA–H and O–H moieties (step 1). Hydrogen activation is facile on Ni-In₂O₃ (E_act,Ni = 19 kJ mol⁻¹) and more difficult on Pd-, Pt- and Rh-In₂O₃ (E_act,Pd = 62 kJ mol⁻¹, E_act,Pt = 48 kJ mol⁻¹ and E_act,Rh = 57 kJ mol⁻¹). Dissociative adsorption of H₂ is more exothermic on Pt- and Rh-In₂O₃ (ΔE_r,Pt = −94 kJ mol⁻¹, ΔE_r,Rh = −144 kJ mol⁻¹) than on Ni- and Pd-In₂O₃ (ΔE_r,Ni = −36 kJ mol⁻¹, ΔE_r,Pd = −10 kJ mol⁻¹). Oxygen vacancy formation proceeds via proton migration from a SA–H to an O–H moiety, resulting in adsorbed H₂O (step 2). This step has similar barriers for Ni-In₂O₃, Pt-In₂O₃ and Rh-In₂O₃ (E_act,Ni = 93 kJ mol⁻¹, E_act,Pt = 68 kJ mol⁻¹ and E_act,Rh = 76 kJ mol⁻¹). On these models, this step is endothermic by 27, 69, and 27 kJ mol⁻¹, respectively. The lowest barrier is obtained for Pd-In₂O₃ (E_act = 47 kJ mol⁻¹; ΔE_r = 80 kJ mol⁻¹). H₂O desorption completes the oxygen vacancy formation pathway (step 3). The highest desorption energy is computed on Ni-In₂O₃ (ΔE_des,Ni = 177 kJ mol⁻¹), while this step is considerably easier on the other SA models (ΔE_des,Pd = 12 kJ mol⁻¹, ΔE_des,Pt = 41 kJ mol⁻¹ and ΔE_des,Rh = 26 kJ mol⁻¹). The overall reaction energy for oxygen vacancy formation with respect to gas-phase H₂ amounts to +150 kJ mol⁻¹ (Ni), −72 kJ mol⁻¹ (Pd), −41 kJ mol⁻¹ (Pt) and +16 kJ mol⁻¹ (Rh). The overall barriers with respect to gas-phase H₂ are 150, 72, 48 and 57 kJ mol⁻¹ for Ni-, Pd-, Pt-, and Rh-In₂O₃, respectively. Notably, the higher overall barrier for Ni is in line with the stronger Ni–O bond as compared to the more noble metals considered here (Table S14†). With respect to the same reference, the unpromoted In₂O₃(111) surface has an overall barrier for oxygen vacancy formation of 67 kJ mol⁻¹.^59,63 The adsorption configurations of the SA on the In₂O₃(111) surface after formation of an oxygen vacancy are reported in Fig. S4.†

Methanol synthesis via formate

For CO₂ hydrogenation, we first discuss the most stable adsorption configuration for CO₂ on the SA-In₂O₃ models. We have explored different adsorption configurations of CO₂ on the SA-In₂O₃ models. The results are reported in Table S5 and Fig. S5,† respectively. On Ni-, Pd- and Pt-In₂O₃ CO₂ adsorbs near the SA with one of its O atoms occupying the oxygen vacancy and the C atom bonding to the SA. On Rh-In₂O₃, the C atom coordinates directly to a lattice O instead of the SA. The adsorption energies of CO₂ are −144, −31, −33 and −87 kJ mol⁻¹ for Ni-, Pd-, Pt- and Rh-In₂O₃, respectively. Thus, the interaction of CO₂ with Ni and Rh is relatively strong compared to Pd and Pt. This is in line with trends found in a previous computational study for CO₂ adsorption on low-index surfaces of transition metals.⁶⁴

Dissociative adsorption of H₂ (step 5) is heterolytic on all SA-In₂O₃ models and leads to O–H and SA–H moieties adjacent to adsorbed CO₂. Homolytic dissociation on top of the metal atom was also investigated but resulted in the formation of a SA–H and a O–H moiety. This state constitutes the initial configuration for CO₂ hydrogenation. Dissociative adsorption is facile on Ni-In₂O₃ (E_act,Ni = 19 kJ mol⁻¹), whereas higher barriers were computed for Pd-, Pt- and Rh-In₂O₃ (E_act,Pd = 104 kJ mol⁻¹, E_act,Pt = 48 kJ mol⁻¹, E_act,Rh = 51 kJ mol⁻¹). Next, methanol can be obtained via a mechanism involving a formate intermediate (steps 6–11 in Fig. 3). The PEDs of the formate pathway for the SA-In₂O₃ models are reported in Fig. 4a. The C atom in CO₂ is hydrogenated by a SA–H to HCO₂ (step 6). The forward activation energies for Ni-, Pd- Pt, and Rh-In₂O₃ are respectively 66, 84, 34 and 42 kJ mol⁻¹. This reaction is endothermic for Pd-, Pt- and Rh-In₂O₃ with reaction energies of 48, 32 and 54 kJ mol⁻¹ and slightly exothermic for Ni-In₂O₃ (ΔE_r = −2 kJ mol⁻¹). The C atom in HCO₂ is protonated again via migration of a O–H to a SA–H moiety, resulting in H₂CO₂ (step 7). Forming the H₂CO₂ intermediate from HCO₂ is endothermic on all SA-In₂O₃ models and involves relatively high activation energies of 114, 220, 146, and 155 kJ mol⁻¹ for Ni-, Pd- Pt- and Rh-In₂O₃, respectively. Notably, higher barriers are associated with migration from O–H to SA–H compared to direct hydrogenation of the C atom by a SA–A moiety. Furthermore, step 7 is the most difficult one in the formate pathway for all SAs. Direct migration of a O–H to hydrogenate the C atom in HCO₂ was not considered, as it will lead to even higher barriers. Next, scission of one of the C–O bonds in the H₂CO₂ intermediate results in CH₂O and O (step 8). This step leads to healing of the oxygen vacancy. This step is endothermic for Ni–, and Pd– and Pt-In₂O₃ surfaces with reaction energies of 68, 54 and 100 kJ mol⁻¹, respectively, whereas it is slightly exothermic for Rh-In₂O₃ (ΔE_r = −3 kJ mol⁻¹). Furthermore, lower barriers are found for Pd– (E_act = 62 kJ mol⁻¹) than for Ni– (E_act = 101 kJ mol⁻¹), Pt– (E_act = 156 kJ mol⁻¹) and Rh-In₂O₃ (E_act = 105 kJ mol⁻¹). Another pathway preceding C–O bond scission involves hydrogenation of HCO₂ to HCO₂H. We could not find a TS for this reaction because of the very weak adsorption of HCO₂H. A similar issue was encountered in investigating the hydrogenation of H₂CO₂ to H₂CO₂H. On all SA-In₂O₃ models, dissociative adsorption of another H₂ leads to a SA–H and O–H moieties adsorbed next to CH₂O (step 9). The CH₂O moiety resulting from H₂CO₂ dissociation is then hydrogenated to CH₃O by proton migration from a SA–H moiety (step 10). This step has an activation energy of 55, 62, 80 and 149 kJ mol⁻¹ for Ni-, Pd- Pt, and Rh-In₂O₃, respectively. Furthermore, it is endothermic for Pd-, Pt and Rh-In₂O₃ (E_act,Pd = 18 kJ mol⁻¹, E_act,Pt = 15 kJ mol⁻¹, E_act,Rh = 14 kJ mol⁻¹), exothermic for Ni-In₂O₃ (E_act,Ni = −77 kJ mol⁻¹). Finally, CH₃OH is obtained in a single concerted step from CH₃O via proton migration from the OH (step 11). This elementary reaction step has activation energies of 61, 106, 118 and 135 kJ mol⁻¹ for Ni-, Pd- Pt, and Rh-In₂O₃, respectively. Furthermore, it is moderately exothermic for Ni-, Pd- and Pt-In₂O₃ (ΔE_r,Ni = −7 kJ mol⁻¹, ΔE_r,Pd = −35 kJ mol⁻¹, ΔE_r,Pt = −13 kJ mol⁻¹), slightly endothermic for Rh-In₂O₃ (ΔE_r,Rh = 16 kJ mol⁻¹).

Fig. 4 Potential energy diagrams of (a) formate pathway to CH₃OH, (b) CO hydrogenation pathway to CH₃OH and (c) redox pathway to CO and (d) carboxyl pathway to CO on different models.

Frei et al. speculated based on potential energy diagrams that the most energetically favourable path to methanol involves three consecutive additions of hydrides and protons and features HCO₂H and CH₂(OH)₂ as intermediates.⁵⁷ However, as pointed out earlier, the very weak adsorption of HCO₂H precludes this pathway on the SA-In₂O₃ models we investigated.

Methanol synthesis via CO. The presence of metal atoms on the surface of In₂O₃ can allow for methanol formation via CO hydrogenation, which does not take place on unpromoted In₂O₃.²³ We thus included this pathway and show the corresponding PEDs in Fig. 4b. After formation of an oxygen vacancy (steps 1–3) and CO₂ adsorption (step 4), direct C–O bond scission can take place (step 12). This step leads to healing of the oxygen vacancy and adsorbed CO and has activation energies of 91, 65, 65 and 81 kJ mol⁻¹ for Ni–, Pd– Pt, and Rh-In₂O₃, respectively. Furthermore, it is exothermic for Pd-, Pt- and Rh-In₂O₃ with reaction energy of −26, −26 and −52, respectively, while it is slightly endothermic for Ni-In₂O₃ (ΔE_r,Ni = 9 kJ mol⁻¹). In the final state, the resulting CO moiety is linearly adsorbed on the SA (SA-C distances: 1.73 Å (Ni), 1.85 Å (Pd), 1.82 Å (Pt) and 1.82 Å (Rh), while the O atom from CO₂ is used to heal the oxygen vacancy in the In₂O₃ surface. We did not find a TS for the two subsequent hydrogenation steps from CO to HCO and H₂CO for Pt- and Rh-In₂O₃. Accordingly, we only report methanol formation via CO for Ni- and Pd-In₂O₃.

For Ni- and Pd-In₂O₃, H₂ is heterolytically dissociated into a SA–H and an O–H moiety (step 13). Such mode of dissociative adsorption of H₂ is more facile on Ni-In₂O₃ (E_act,Ni = 19 kJ mol⁻¹) than on Pd-In₂O₃ (E_act,Pd = 41 kJ mol⁻¹). Such state where CO is adsorbed next to a SA–H and a O–H moiety constitutes the initial state of the CO hydrogenation pathway. CO hydrogenation to HCO proceeds using the H adsorbed on the SA (step 14) with activation energies of 71 kJ mol⁻¹ (Ni-In₂O₃) and 148 kJ mol⁻¹ (Pd-In₂O₃). In line with this difference, the reaction is much more endothermic for Pd (ΔE_r,Pd = 93 kJ mol⁻¹) than for Ni (ΔE_r,Ni = 6 kJ mol⁻¹). Subsequently, HCO is hydrogenated to CH₂O (step 15) with barriers of 177 kJ mol⁻¹ (Ni-In₂O₃) and 128 kJ mol⁻¹ (Pd-In₂O₃). We found that the preferential pathway involves migration of H from the OH moiety to the SA, followed by hydrogenation of HCO to CH₂O. Direct migration of a OH to hydrogenate the C atom in HCO would face a higher barrier and was therefore not studied. This elementary step is endothermic for both SA-In₂O₃ models (ΔE_r,Ni = 142, ΔE_r,Pd = 114 kJ mol⁻¹, respectively). The resulting CH₂O intermediate also occurs in the formate pathway to methanol discussed above. Thus, further hydrogenation of CH₂O to methanol will occur via steps 9 and 10.

CO formation via rWGS

On SA-In₂O₃ models, the formation of CO can occur either via direct C–O bond cleavage of CO₂ (redox pathway, steps 4 and 20) or via a H-assisted pathway involving a COOH intermediate (carboxyl pathway, steps 16–19). In the redox pathway, after formation of H₂O (steps 1–3) and CO₂ adsorption (step 4) one of the C–O bonds in CO₂ is cleaved (step 12). These steps have been previously discussed. The desorption of CO (step 20) closes the rWGS catalytic cycle. The highest barriers for CO desorption are found on Pt– (ΔE_des,Pt = 158 kJ mol⁻¹) and Rh-In₂O₃ (ΔE_des,Rh = 215 kJ mol⁻¹) whereas lower values are found on Ni– (ΔE_des,Ni = 121 kJ mol⁻¹) and Pd-In₂O₃ (ΔE_des,Pd = 31 kJ mol⁻¹). Concerning the unpromoted In₂O₃, Ye et al. reported that the process replenishing an oxygen vacancy by direct cleavage of the C–O bond in CO₂ has a barrier of 134 kJ mol⁻¹ on an In₂O₃(110) model surface.⁶¹

The carboxyl pathway proceeds by protonation of adsorbed CO₂ to form COOH by a SA–H moiety (step 16). The highest activation energies are computed for Ni-In₂O₃ (E_act = 82 kJ mol⁻¹) and Pt-In₂O₃ (E_act = 147 kJ mol⁻¹) followed by Pd-In₂O₃ (E_act = 61 kJ mol⁻¹) and Rh-In₂O₃ (E_act = 59 kJ mol⁻¹). Furthermore, this elementary step is endothermic for Ni- and Pt-In₂O₃ (ΔE_r,Ni = 74 kJ mol⁻¹, ΔE_r,Pt = 55 kJ mol⁻¹), whereas it is exothermic for Pd- and Rh-In₂O₃ (ΔE_r,Pd = −20 kJ mol⁻¹, E_r,Rh = −78 kJ mol⁻¹). Formation a COOH intermediate is associated with higher barriers compared to formation of HCOO, in line with an earlier study on In₂O₃.⁶⁵ Next, COOH dissociates into CO and OH (step 17). This step features activation energies of 45, 42, 65, and 25 kJ mol⁻¹ for Ni-, Pd-, Pt- and Rh-In₂O₃, respectively. Furthermore, it is endothermic for all SA-In₂O₃ surfaces (ΔE_r,Ni = 17 kJ mol⁻¹, ΔE_r,Pd = 19 kJ mol⁻¹, ΔE_r,Pt = 33 kJ mol⁻¹, and ΔE_r,Rh = 20 kJ mol⁻¹). Water is formed via proton transfer to the OH moiety obtained from COOH dissociation (step 18). This step features activation energies of 140, 140, 133 and 186 kJ mol⁻¹ for Ni-, Pd-, Pt- and Rh-In₂O₃, respectively. Moreover, it is endothermic for Rh-, and Pd-In₂O₃ (ΔE_r,Rh = 107 kJ mol⁻¹, and ΔE_r,Pd = 58 kJ mol⁻¹), and exothermic for Ni– and Pt-In₂O₃ (ΔE_r,Ni = −58 kJ mol⁻¹, and ΔE_r,Pt = −16 kJ mol⁻¹). Notably, this step features the highest activation energy of the carboxyl pathway for all SA-In₂O₃. The elementary step of COOH dissociation on the bare In₂O₃(111) surface is associated with a relatively high barrier (146 kJ mol⁻¹), similar to the ones found on our SA-In₂O₃ models. Desorption of CO (step 19) features a barrier of 51, 25, 24 and 120 kJ mol⁻¹ for Ni-, Pd-, Pt- and Rh-In₂O₃, respectively, closing the carboxyl cycle.

Microkinetic simulations

Overall kinetics. To compare the catalytic activities of the various SA-In₂O₃ model surfaces, we calculate the CO₂ hydrogenation reaction rate using a microkinetic model based on DFT-computed reaction energetics. The active sites in our model consist of isolated single atoms, stabilized on the In₂O₃ support. We do not take migration of intermediates between different active sites into account. In other words, co-adsorbed species are modelled as distinct variations of a single active site. In this approximation, all elementary reaction steps are unimolecular, with exception of the adsorption and desorption steps. A detailed list of the elementary reaction steps is provided in the ESI† (Tables S10–13).

The CO₂ consumption rate and CH₃OH selectivity as a function of temperature for each SA-In₂O₃ model are displayed in Fig. 5. The turnover frequencies (TOF) towards CH₃OH and CO (TOF_CH3OH and TOF_CO, respectively) are reported in Fig. S6.† We also constructed a microkinetic model for unpromoted In₂O₃ based on the published DFT data of Frei et al.⁵⁷Fig. 5a shows that the CO₂ conversion rate decreases in the order Pd > Ni > Pt > Rh. Furthermore, none of the SA-In₂O₃ models show any appreciable methanol selectivity (Fig. 5b), which is due to the much lower TOF_CH3OH in comparison to TOF_CO (Fig. S6c and d†).

Fig. 5 (a) CO₂ consumption rate (s⁻¹) and (b) CH₃OH selectivity as a function of temperature on SA-In₂O₃ models (p = 50 bar, H₂/CO₂ ratio =5).

Compared to bare In₂O₃, the adsorption of a SA on the In₂O₃(111) surface results in lower CO₂ consumption rates (Fig. S6a†) and TOF_CH3OH for T < 300 °C (Fig. S6c†). In the same range, the TOF_CO decreases in the order Pd > Ni > bare In₂O₃ > Pt > Rh (Fig. S6d†).

The analysis of the surface state, apparent activation energy and reaction orders as function of temperature is given in Fig. 6. The coverages in Fig. 6 should be interpreted as the fraction of time the system spends in a particular state (i.e., the time average). According to the ergodicity principle in statistical thermodynamics, this equals the fraction of active sites that is in a particular state (i.e., the ensemble average).

Fig. 6 Surface state of the catalyst as a function of temperature for (a) Ni–, (b) Pd– (e) Pt– and (f) Rh-In₂O₃, respectively. Apparent activation energy (in kJ mol⁻¹) and reaction orders in CO₂ and H₂ as a function of temperature on (c) Ni–, (d) Pd– (g) Pt– and (h) Rh-In₂O₃, respectively.

On the Ni-In₂O₃ surface (Fig. 6a), the dominant surface state at lower temperatures features two adsorbed H species (2H). The reason behind this can be understood from Fig. 3a. Activation of H₂ and subsequent H₂O formation (steps 1–2) are facile (E_act = 19 kJ mol⁻¹ and E_act = 74 kJ mol⁻¹, respectively), whereas water desorption (step 3) features a high barrier (ΔE_des = 177 kJ mol⁻¹). The latter step is the most difficult one in the Ov formation pathway and limits the progress of the reaction. As water desorption is facilitated by a higher temperature, it leads to an increase in the number of oxygen vacancies and CO₂ adsorption. Consequently, the coverage of the 2H state decreases and the coverage of intermediate states in the formate pathway to methanol (CO₂ + 2H and HCO₂ + H) increases. At T > 400 °C, the latter states become dominant. On the Ni-In₂O₃ model, hydrogenation of CO₂ to HCO₂ (step 6) has an activation energy of 66 kJ mol⁻¹ and a reaction energy of −2 kJ mol⁻¹. In contrast, subsequent hydrogenation to H₂CO₂ (step 7) has a considerably higher activation energy of 114 kJ mol⁻¹ and is endothermic. As formation of H₂CO₂ is difficult and the previous two states are comparable in terms of energy, CO₂ + 2H and HCO₂ + H are the dominant states. At low temperature, the model predicts a reaction order of nearly zero in CO₂ and H₂ (Fig. 6c). Under these conditions, the reaction is limited by the rate of H₂O removal, and an increase of the partial pressure of CO₂ or H₂ does not affect this process. Thus, the corresponding reaction orders in CO₂ and H₂ must be zero. At higher temperature, however, a negative reaction order in H₂ is found. A higher partial pressure of H₂ would lead to higher rates of H adsorption facilitating hydrogenation reactions. In turn, this pushes the reaction away from the CO product leading to a decrease in the overall rate. On the Ni-In₂O₃ model, the apparent activation energy (E_app, Fig. 6c) has a constant value of at 175 kJ mol⁻¹ between 200 °C and 400 °C. This value corresponds to the desorption energy of H₂O. At higher temperature, a slight decrease in the E_app is observed, indicating a change in the rate-limiting step. This aspect will be clarified with a DRC analysis (vide infra).

On the Pd-In₂O₃ model, the COOH + H is the dominant state at T < 450 °C (Fig. 6b). Fig. 3b shows that hydrogenation of CO₂ to COOH (step 16) and its subsequent dissociation into CO and OH are facile (E_act = 61 kJ mol⁻¹ and E_act = 42 kJ mol⁻¹, respectively). However, the subsequent step of OH hydrogenation has a higher activation energy and is endothermic (E_act = 140 kJ mol⁻¹ and E_R = 58 kJ mol⁻¹, respectively). This makes the conversion of the COOH + H state difficult. At high temperature (T > 450 °C), the barrier for OH hydrogenation (step 17) can be more easily overcome, resulting in formation of oxygen vacancies. As a result, the dominant state becomes the Pd-In₂O₃ surface with an oxygen vacancy (Ov). The model predicts a reaction order of zero in CO₂ and H₂ for T < 450 °C (Fig. 6d). As the dominant working state already corresponds to intermediates derived from H₂ and CO₂, a change in the partial pressure of these reactants does not impact the activity. For T > 450 °C, the reaction orders in CO₂ and H₂ are positive, indicating that hydrogenation reactions are controlling the reaction rate. The E_app is approximately 150 kJ mol⁻¹ at 200 °C and decreases with increasing temperature (Fig. 6d).

On the Pt-In₂O₃ model, the HCO₂+H state is the dominant state for below 525 °C (Fig. 6e). Hydrogenation of CO₂ to HCO₂ (step 6 in Fig. 3c) is facile (E_act = 34 kJ mol⁻¹), while the subsequent step of HCO₂ hydrogenation to H₂CO₂ (step 7) is more difficult (E_act = 146 kJ mol⁻¹, E_R = +129 kJ mol⁻¹). This hampers the conversion of the HCO₂ + H state. At higher temperature, the surface coverage of the HCO₂ + H state decreases in favor of oxygen vacancies (Ov). The Ov state becomes dominant for T > 525 °C. The reaction order in CO₂ is zero at low temperature, because the surface is already covered with HCO₂+H resulting from CO₂ hydrogenation (Fig. 6g). At higher temperature (T > 475 °C), the reaction order in CO₂ becomes positive indicating that CO₂ activation reactions are controlling the rate. The reaction order in H₂ is negative, as also found for Ni– and Pd-In₂O₃. The E_app stays almost constant for temperature below 300 °C at 200 kJ mol⁻¹ and decreases with increasing temperature (Fig. 6g).

Within the temperature range of 200 °C and 600 °C, the dominant state of Rh-In₂O₃ is the COOH+H state (Fig. 6f). The reason becomes apparent from inspection of Fig. 3d. Hydrogenation of CO₂ to COOH and its further dissociation (step 16–17) are facile (E_act = 59 and E_act = 25 kJ mol⁻¹, respectively). However, the subsequent step of OH hydrogenation is more difficult (E_act = 186 kJ mol⁻¹) and hampers the conversion of the CO + OH state. With respect to the dominant COOH + H state, an overall barrier of 341 kJ mol⁻¹ must be overcome to complete the H-assisted rWGS cycle (Fig. 4d). As a result, only for temperatures far above 600 °C, states other than COOH + H observed. The reaction order in CO₂ is zero at low temperature, because the surface is already covered with COOH + H resulting from CO₂ hydrogenation (Fig. 6h). The reaction order in H₂ is negative, indicating that a higher partial pressure of H₂ would push the reaction away from the dominant CO formation pathway. The apparent activation energy decreases from 260 kJ mol⁻¹ to 220 kJ mol⁻¹ between 200 and 400 °C (Fig. 6h).

Sensitivity analysis

In this section, we analyse in more depth the steps that lead to CH₃OH and CO during CO₂ hydrogenation on the SA-promoted In₂O₃ models. We identify the elementary steps that control the overall CO₂ consumption rate and the CH₃OH selectivity and investigate how these steps change with reaction temperature. For this purpose, we conduct a sensitivity analysis based on the degree of rate control (DRC) analysis.⁵⁴ At zero CO₂ conversion, the sum of the DRC coefficients is conserved at unity.⁵⁵

The DRC analysis and reaction flux analysis for the Ni-In₂O₃ model are shown in Fig. 7a and b, respectively. At low temperature, the CO₂ consumption rate is mostly controlled by water desorption (step 3) preceding oxygen vacancy formation (Fig. 7a). This result is in line with the dominant surface state being the 2H state, which precedes H₂O formation (Fig. 6a). Furthermore, H₂O desorption has the highest barrier in the oxygen vacancy formation pathway (ΔE_des = 177 kJ mol⁻¹), making formation of an oxygen vacancy difficult. At higher temperature (T > 275 °C), water desorption results in oxygen vacancies and CO₂ adsorption. As a result, the DRC coefficient of the H₂O desorption step decreases, in favor of the CO₂ dissociation step (step 12). At higher temperature, the latter becomes the main rate-limiting step. An increased rate of CO₂ dissociation would result in a higher flux towards CO, benefiting the formation of this product. The dominant pathway (Fig. 7b) on the Ni-In₂O₃ model proceeds via formation of an oxygen vacancy followed by CO₂ adsorption and direct dissociation. The carboxyl pathway to CO is not taken because OH removal (step 18) is associated with a higher barrier than direct CO₂ dissociation (E_act = 140 kJ mol⁻¹ and E_act = 81 kJ mol⁻¹, respectively). Compared to the dominant pathway to CO, methanol formation pathways feature considerably lower molar fluxes. This is in keeping with our finding that Ni-In₂O₃ mainly produces CO. The formate pathway is not taken because the steps of H₂CO hydrogenation to H₂CO₂ and its further dissociation into CH₂O and O (steps 7–8) result in a combined reaction energy of +180 kJ mol⁻¹. Hydrogenation of CO to methanol does not proceed because CH₂O formation from HCO has a high activation energy and is endothermic (E_act = 177 kJ mol⁻¹ and E_r = 142 kJ mol⁻¹, respectively).

Fig. 7 Degree of rate control and molar fluxes analysis for (a and b) Ni-, (c and d) Pd, (e and f) Pt and (g and h) Rh-In₂O₃ models. The numbers in the arrows are molar reaction rates (s⁻¹) and are normalized with respect to the amount of adsorbed CO₂. The dominant pathway for each model is indicated by the colored pathways. The molar fluxes are reported at T = 200 °C.

For Pd-In₂O₃ (Fig. 7c), for T < 425 °C, the overall rate is mostly controlled by CO₂ dissociation (step 12, E_act = 65 kJ mol⁻¹). The dominant pathway on the Pd-In₂O₃ surface proceeds via formation of an oxygen vacancy and direct dissociation of CO₂ (Fig. 7d). At higher temperature (T > 400 °C), a decrease in the DRC coefficient of step 12 is observed in favor of the step of OH removal in the carboxyl pathway (step 18, E_act = 140 kJ mol⁻¹) which becomes the main rate limiting step. The change in rate limiting steps indicates a change in the dominant pathway. Indeed, for T > 425 C, the carboxyl pathway is the dominant pathway of CO formation (Fig. S7†). Furthermore, this is in line with the increasing reaction order in H₂ with increasing temperature (Fig. 6d). Compared to the dominant pathway to CO, methanol formation pathways have considerably lower molar fluxes. The formate pathway is not taken because the step of HCO₂ hydrogenation to H₂CO₂ has a high activation energy (E_act = 220 kJ mol⁻¹) and is endothermic (E_r = 167 kJ mol⁻¹). Furthermore, hydrogenation of CO to CH₃OH does not take place because the hydrogenation of CO to HCO has an activation energy of 148 kJ mol⁻¹, significantly higher than the CO desorption energy (31 kJ mol⁻¹). This makes further CO hydrogenation less favorable than its desorption.

The DRC analysis and reaction flux analysis for the Pt-In₂O₃ model are shown in Fig. 7e and f, respectively. At low temperature, the overall rate of is mainly controlled by CO desorption (step 20, ΔE_des = 158 kJ mol⁻¹). OH hydrogenation (step 18, E_act = 133 kJ mol⁻¹) controls the rate to a smaller extent as can be seen by its lower DRC coefficient. In line with Ni-In₂O₃, the dominant pathway proceeds via formation of an oxygen vacancy and direct dissociation of CO₂ (Fig. 7f). The carboxyl pathway is not taken because COOH formation and OH removal are associated with high barriers (E_act = 147 kJ mol⁻¹ and E_act = 133 kJ mol⁻¹, respectively). With increasing temperature, the DRC coefficient of CO desorption decreases in favor of CO₂ dissociation (step 12). This is in line with the positive reaction order in CO₂ calculated at higher temperature (Fig. 6e). Compared to the dominant pathway to CO, hydrogenation of CO₂ to methanol is associated with significantly lower molar fluxes. The formate pathway is not taken because the hydrogenation of HCO₂ to H₂CO₂ has a high activation energy and is endothermic (E_act = 146 kJ mol⁻¹, E_r = 129 kJ mol⁻¹).

On Rh-In₂O₃, for T < 300 °C, the rate of CO₂ consumption is mainly controlled by CO desorption (step 20, ΔE_des = 215 kJ mol⁻¹). At higher temperature, the DRC coefficient of this elementary step decreases in favour of CO₂ dissociation (step 12, E_act = 81 kJ mol⁻¹) that becomes the rate determining step (RDS) at T > 450 °C. In line with the other SA-In₂O₃ models, the dominant pathway proceeds via formation of an oxygen vacancy and direct dissociation of CO₂ (Fig. 7h). The high barrier associated with OH hydrogenation (step 18, E_act = 186 kJ mol⁻¹) hampers the carboxyl pathway. Compared to the dominant pathway to CO, methanol synthesis from CO₂ is associated with significantly lower molar fluxes. The formate pathway is not taken because the hydrogenation of HCO₂ to H₂CO₂ has a high activation energy and is endothermic (E_act = 155 kJ mol⁻¹, E_r = 136 kJ mol⁻¹).

General discussion

The present study clarifies the role of single atoms of Pt, Pd, Ni and Rh adsorbed on the In₂O₃(111) surface in CO₂ hydrogenation. The explored reaction mechanism includes a direct route for CO₂ hydrogenation to methanol (formate pathway), a pathway to methanol via CO hydrogenation and the competing reverse water-gas shift (rWGS) reaction. Earlier studies pointed out that SAs on the In₂O₃ surface do not catalyse CH₄ formation.^29,30,35 A key result of this study is that CO is the dominant reaction product for all SA-In₂O₃ models. CO is obtained via a redox mechanism involving oxygen vacancy formation followed by CO₂ adsorption and direct C–O bond cleavage. On all models, the CH₃OH selectivity was found to be negligibly small, because barriers associated with the hydrogenation of formate intermediates to methanol are significantly higher than those for the dominant redox pathway to CO.

Methanol formation on In₂O₃ involves the formation of oxygen vacancies and hydrogenation of formate intermediates.^63,66 On the unpromoted In₂O₃(111) surface, oxygen vacancy formation has an overall barrier of 67 kJ mol⁻¹ with respect to gas-phase H₂.^59,63 When a single Pd, Pt or Rh atom is adsorbed on In₂O₃, oxygen vacancy formation exhibits comparable overall barriers of respectively 62, 48 and 57 kJ mol⁻¹. In contrast, a much higher overall barrier of 150 kJ mol⁻¹ is computed for the single Ni atom on In₂O₃(111). Concerning the conversion of formate, HCO₂ hydrogenation to H₂CO₂ on an In₂O₃(110) surface was earlier found to exhibit an activation energy of 55 kJ mol⁻¹ and a reaction energy of +5 kJ mol⁻¹.⁶⁰ On In₂O₃(111), similar activation and reaction energies are found (98 and 15 kJ mol⁻¹, respectively).⁵⁷ All SA-In₂O₃ models show higher activation energies for this step (step 7 in Fig. 3), which is the main cause of the very low methanol formation rate. We performed a sensitivity analysis of the microkinetic simulations involving decreasing the activation energies of steps 7 by 10, 20, 30 and 40 kJ mol⁻¹. The results given in Fig. S8† show that decreasing barriers result in a higher methanol selectivity for Ni-, Pt- and Rh-In₂O₃. Herein, we consider methanol selectivity of at least 1% as significant. Methanol formation is least improved for Pd-In₂O₃, because the activation energy associated with HCO₂ hydrogenation (step 7) is very high (220 kJ mol⁻¹). In this case, a decrease of at least 80 kJ mol⁻¹ is required to obtain significant CH₃OH selectivity. We extended the sensitivity analysis to include H₂CO₂ dissociation (step 8) and reported the results in Fig. S9.† Also in this case, decreasing barriers results in a higher methanol selectivity for Ni-, Pt- and Rh-In₂O₃, whereas methanol formation is less improved for Pd-In₂O₃ because the previous step has a higher barrier (220 kJ mol⁻¹ and 62 kJ mol⁻¹ for steps 7–8, respectively).

Our microkinetic simulations show that all SA-In₂O₃ have CO as the main product obtained via a redox mechanism. The activity for the rWGS reaction decreases in the order Pd > Ni > bare In₂O₃ > Pt > Rh (Fig. S6d†). We now explain this finding based on the DFT-computed activation energies and the analysis of the rate limiting steps. On all SA-In₂O₃ models, direct CO₂ dissociation (step 12) and subsequent CO desorption (step 20) are the main rate-limiting steps towards CO formation. On Pd-In₂O₃, these steps have barriers of 65 and 31 kJ mol⁻¹, respectively. On Ni-In₂O₃, direct CO₂ dissociation and CO desorption are associated with barriers of 91 and 80 kJ mol⁻¹, respectively. On bare-In₂O₃, CO is obtained via carboxyl pathway involving a COOH intermediate. This pathway is mainly limited by the rate of COOH dissociation which is associated with a high barrier (E_act = 150 kJ mol⁻¹). Compared to Ni and Pd-In₂O₃, the barrier for CO₂ dissociation on Pt-In₂O₃ is similar (E_act = 65 kJ mol⁻¹), but CO desorption is associated with a higher barrier (ΔE_des = 158 kJ mol⁻¹). On Rh-In₂O₃, desorption of CO is associated with a high desorption energy (ΔE_des = 215 kJ mol⁻¹) hampering its formation. For this reason, this model shows the lowest activity to CO.

We know explain why the redox pathway is preferred over the carboxyl pathway for the formation of CO on our SA-In₂O₃ models. In the redox pathway, the C–O bond in CO₂ is broken directly resulting in CO binding to the SA and the O healing the vacancy. This process is associated with relatively low barriers on all models, in line with the strong SA–O bond (Table S14†). Notably, on a bare In₂O₃(110) surface, where CH₃OH is the main product, replenishing the oxygen vacancy by CO₂ dissociation is associated with a reaction energy of 1.4 eV.⁶⁰ This indicates that direct CO₂ dissociation would have a barrier of at least 1.4 eV, making this elementary reaction step unlikely. In contrast with SA-promoted In₂O₃, on bare-In₂O₃, there is no active site available for accepting the CO molecule resulting from CO₂ dissociation. On the SA-In₂O₃ models, breaking the C–O bond on a COOH intermediate is also possible, however the subsequent step of OH removal via H₂O formation is associated with high barriers because it requires the migration of a OH moiety to the SA.

A fundamental assumption in this study was that each kinetic pathway features one oxygen vacancy, in line with previous DFT calculations on In₂O₃.^19,21,22 To investigate whether the presence of a second oxygen vacancy could change the conclusions of this work, we calculated a few key elementary steps relevant to such a mechanism. The results are reported in Tables S15–S18† reveal that a second vacancy provides a lower activation energy for some elementary reaction steps. However, these deviations are not significant enough to impact the dominant pathway identified in our microkinetic simulations. Consequently, we deduce that the inclusion of a second oxygen vacancy does not change the conclusions put forth in this study.

We also compare our results to those recently reported in the literature. In line with the present work, Frei et al. concluded from DFT calculations that a Ni SA on In₂O₃ would be active for the rWGS, because it features lower barriers compared to CO₂ hydrogenation to methanol.²⁹ Furthermore, Millet et al. showed that isolated Ni atoms in a Ni₁–Mg_1−xO catalyst are active for the rWGS but are unable to catalyse CO₂ hydrogenation to CH₄ (or MeOH), for which Ni clusters are needed.⁶⁷ Indeed, Shen et al. deployed DFT to compare various pathways of CO₂ hydrogenation to methanol on a Ni₄/In₂O₃(111) model and found that methanol can be obtained via CO intermediate obtained by rWGS.⁶⁸ In an earlier computational work, we reported that single atoms of Ni either doped in or adsorbed on In₂O₃(111) surface would mainly catalyze CO formation, whereas small Ni clusters would mainly lead to methanol.²⁷ In line with this earlier work, we indicate here that the Ni-In₂O₃ model features low-barrier CO₂ dissociation resulting in a preference for CO formation over hydrogenation to methanol. Ye et al. reported that methanol is the main product of CO₂ hydrogenation via formate on a Pd₄/In₂O₃(111) model.⁶² Our Pd-In₂O₃ model features higher barriers for the formate pathway compared to direct CO₂ dissociation precluding methanol formation in favor of CO. Liu et al. reported that methanol formation on a Pt₄/In₂O₃(111) model takes place via hydrogenation of CO resulting from CO₂ dissociation.³² Our Pt-In₂O₃ model does not allow a TS for the subsequent hydrogenation steps of CO to HCO and H₂CO precluding methanol formation from CO hydrogenation. Clusters of Pd, Pt or Ni on the In₂O₃(111) surface catalyze the formation of methanol because hydrogenation of either formate or CO is favored over the rWGS reaction to CO. However, in the limit of a single atom of Pd, Pt or Ni on the In₂O₃(111) surface, pathways leading to methanol have higher barriers than the redox pathway to CO resulting in low methanol selectivity. In a recent theoretical study, Pinheiro Araújo et al. studied atomically dispersed metal atoms (Pd, Pt, Rh, Ni, Co, Au, Ir) on the In₂O₃ surface.²⁸ They speculated on the basis of energy diagrams that CO₂ hydrogenation to methanol via formate is favored over CO formation. Although this conclusion is at odds with our finding from microkinetic simulations that CO₂ hydrogenation on SA-promoted In₂O₃ would mainly lead to CO, it can be understood from the fact that the study of Pinheiro Araújo et al. only considered the COOH pathway to CO, i.e., they omitted the redox pathway.

Conclusions

Using microkinetic models and DFT calculations, we clarified the role of single atoms (SA) on In₂O₃ in CO₂ hydrogenation to CH₃OH and CO. We investigated the role of single atoms of Pt, Pd, Ni and Rh on the (111) surface of In₂O₃ for the formation of methanol via formate and CO as well as the competing rWGS reaction via redox and carboxyl pathways. Compared to the pristine In₂O₃(111) surface, adsorbing a single Pd, Pt or Rh atom results in comparable overall barriers for oxygen vacancy formation. In contrast, the presence of a single Ni atom on the In₂O₃ surface increases the overall barrier for oxygen vacancy formation significantly. Microkinetic simulations reveal that all SA-In₂O₃ models mainly catalyze CO via formation of an oxygen vacancy followed by CO₂ adsorption and direct C–O bond cleavage (redox pathway). The carboxyl pathway to CO is not taken, because OH removal resulting from COOH dissociation is associated with higher barriers compared to the dominant redox pathway. The model predicts that Pd is the preferred promoter for In₂O₃ model, achieving the highest activity among the transition metals considered due to a low barrier for CO₂ dissociation and a relatively weak adsorption strength of the CO product. For all models, the CH₃OH selectivity is negligibly low because high barriers are associated with the hydrogenation of HCO₂ to H₂CO₂ and its subsequent dissociation into CH₂O (formate pathway) compared to the dominant rWGS pathway.

Author contributions

Francesco Cannizzaro: DFT calculations, microkinetic modelling, writing and visualization. Sjoerd Kurstjens and Tom van den Berg: DFT calculations. Emiel J.M. Hensen and Ivo Filot resources, conceptualization, supervision.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge NWO and SurfSARA for providing access to computational resources used to carry out the DFT calculations reported in this work. This work was supported by the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation program funded by the Ministry of Education, Culture and Science of the government of the Netherlands. This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 801359.

Notes and references

1. A. González-Garay, M. S. Frei, A. Al-Qahtani, C. Mondelli, G. Guillén-Gosálbez and J. Pérez-Ramírez, Energy Environ. Sci., 2019, 12, 3425–3436.
2. E. V. Kondratenko, G. Mul, J. Baltrusaitis, G. O. Larrazábal and J. Pérez-Ramírez, Energy Environ. Sci., 2013, 6, 3112.
3. V. Dias, M. Pochet, F. Contino and H. Jeanmart, Front. Mech. Eng., 2020, 6, 21.
4. A. Álvarez, A. Bansode, A. Urakawa, A. V. Bavykina, T. A. Wezendonk, M. Makkee, J. Gascon and F. Kapteijn, Chem. Rev., 2017, 117, 9804–9838.
5. J. Zhong, X. Yang, Z. Wu, B. Liang, Y. Huang and T. Zhang, Chem. Soc. Rev., 2020, 49, 1385–1413.
6. X. Jiang, X. Nie, X. Guo, C. Song and J. G. Chen, Chem. Rev., 2020, 120, 7984–8034.
7. J. Wang, G. Zhang, J. Zhu, X. Zhang, F. Ding, A. Zhang, X. Guo and C. Song, ACS Catal., 2021, 2020, 1406–1423.
8. Y. A. Daza and J. N. Kuhn, RSC Adv., 2016, 6, 49675–49691.
9. S. Saeidi, N. A. S. Amin and M. R. Rahimpour, J. CO2 Util., 2014, 5, 66–81.
10. W. Li, H. Wang, X. Jiang, J. Zhu, Z. Liu, X. Guo and C. Song, RSC Adv., 2018, 8, 7651–7669.
11. S. Gumber and A. V. P. Gurumoorthy, Methanol Sci. Eng., 2018, pp. 661–674.
12. G. A. Olah, Angew. Chem., Int. Ed., 2005, 44, 2636–2639.
13. A. Goeppert, M. Czaun, J. P. Jones, G. K. Surya Prakash and G. A. Olah, Chem. Soc. Rev., 2014, 43, 7995–8048.
14. C. Baltes, S. Vukojević and F. Schüth, J. Catal., 2008, 258, 334–344.
15. M. Behrens, Angew. Chem., Int. Ed., 2016, 55, 14906–14908.
16. J. Nakamura, Y. Choi and T. Fujitani, Top. Catal., 2003, 22, 277–285.
17. M. B. Fichtl, D. Schlereth, N. Jacobsen, I. Kasatkin, J. Schumann, M. Behrens, R. Schlögl and O. Hinrichsen, Appl. Catal., A, 2015, 502, 262–270.
18. X. Jiang, X. Nie, X. Guo, C. Song and J. G. Chen, Chem. Rev., 2020, 120, 7984–8034.
19. J. Ye, C. Liu and Q. Ge, J. Phys. Chem. C, 2012, 116, 7817–7825.
20. M. S. Frei, C. Mondelli, A. Cesarini, F. Krumeich, R. Hauert, J. A. Stewart, D. Curulla Ferré and J. Pérez-Ramírez, ACS Catal., 2020, 10, 1133–1145.
21. Q. Sun, J. Ye, C. Liu and Q. Ge, Greenhouse Gases: Sci. Technol., 2014, 4, 140–144.
22. J. Ye, C. Liu, D. Mei and Q. Ge, ACS Catal., 2013, 3, 1296–1306.
23. A. Cao, Z. Wang, H. Li and J. K. Nørskov, ACS Catal., 2021, 11, 1780–1786.
24. J. Wang, G. Zhang, J. Zhu, X. Zhang, F. Ding, A. Zhang, X. Guo and C. Song, ACS Catal., 2021, 11, 1406–1423.
25. A. Posada-Borbón and H. Grönbeck, ACS Catal., 2021, 11, 9996–10006.
26. J. Zhu, F. Cannizzaro, L. Liu, H. Zhang, N. Kosinov, I. A. W. Filot, J. Rabeah, A. Brückner and E. J. M. Hensen, ACS Catal., 2021, 11, 11371–11384.
27. F. Cannizzaro, E. J. M. Hensen and I. A. W. Filot, ACS Catal., 2023, 13, 1875–1892.
28. T. Pinheiro Araújo, J. Morales-Vidal, T. Zou, R. García-Muelas, P. O. Willi, K. M. Engel, O. V. Safonova, D. Faust Akl, F. Krumeich, R. N. Grass, C. Mondelli, N. López and J. Pérez-Ramírez, Adv. Energy Mater., 2022, 12, 14.
29. M. S. Frei, C. Mondelli, R. García-Muelas, J. Morales-Vidal, M. Philipp, O. V. Safonova, N. López, J. A. Stewart, D. C. Ferré and J. Pérez-Ramírez, Nat. Commun., 2021, 12, 2–12.
30. M. S. Frei, C. Mondelli, R. García-Muelas, K. S. Kley, B. Puértolas, N. López, O. V. Safonova, J. A. Stewart, D. Curulla Ferré and J. Pérez-Ramírez, Nat. Commun., 2019, 10, 1–11.
31. Z. Han, C. Tang, J. Wang, L. Li and C. Li, J. Catal., 2021, 394, 236–244.
32. K. Sun, N. Rui, C. Shen and C. J. Liu, J. Phys. Chem. C, 2021, 125, 10926–10936.
33. J. Wang, K. Sun, X. Jia and C. J. Liu, Catal. Today, 2021, 365, 341–347.
34. N. H. M. Dostagir, C. Thompson, H. Kobayashi, A. M. Karim, A. Fukuoka and A. Shrotri, Catal. Sci. Technol., 2020, 10, 8196–8202.
35. N. H. M. Dostagir, C. Thompson, H. Kobayashi, A. M. Karim, A. Fukuoka and A. Shrotri, Catal. Sci. Technol., 2020, 10, 8196–8202.
36. X. Jia, K. Sun, J. Wang, C. Shen and C. J. Liu, J. Energy Chem., 2020, 50, 409–415.
37. J. Ye, C. J. Liu, D. Mei and Q. Ge, J. Catal., 2014, 317, 453.
38. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
39. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
40. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
41. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
42. G. Henkelman and H. Jónsson, J. Chem. Phys., 2000, 113, 9978–9985.
43. D. Heidrich and W. Quapp, Theor. Chim. Acta, 1986, 70, 89–98.
44. R. Nelson, C. Ertural, J. George, V. L. Deringer, G. Hautier and R. Dronskowski, J. Comput. Chem., 2020, 41, 1931–1940.
45. S. Maintz, V. L. Deringer, A. L. Tchougr and R. Dronskowski, J. Comput. Chem., 2016, 37, 11–12.
46. M. Yu and D. R. Trinkle, J. Chem. Phys., 2011, 134, 064111.
47. P. N. Brown, G. D. Byrne and A. C. Hindmarsh, SIAM Journal on Scientific and Statistical Computing, 2006, 10, 1038–1051.
48. G. D. Byrne and A. C. Hindmarsh, J. Comput. Phys., 1987, 70, 1–62.
49. G. D. Byrne and A. C. Hindmarsh, ACM Trans. Math. Softw., 1975, 1, 71–96.
50. C. H. Shomate, J. Phys. Chem., 2002, 58, 368–372.
51. NIST Chemistry WebBook, https://webbook.nist.gov/chemistry/, (accessed 1 February 2022).
52. S. Kozuch and S. Shaik, J. Am. Chem. Soc., 2006, 128, 3355–3365.
53. S. Kozuch and S. Shaik, J. Phys. Chem. A, 2008, 112, 6032–6041.
54. C. T. Campbell, ACS Catal., 2017, 7, 2770–2779.
55. A. W. Filot, Introduction to Microkinetic Modeling, Technische Universiteit Eindhoven, 2017.
56. C. Stegelmann, N. C. Schiødt, C. T. Campbell and P. Stoltze, J. Catal., 2004, 221, 630–649.
57. M. S. Frei, M. Capdevila-Cortada, R. García-Muelas, C. Mondelli, N. López, J. A. Stewart, D. Curulla Ferré and J. Pérez-Ramírez, J. Catal., 2018, 361, 313–321.
58. M. S. Frei, C. Mondelli, R. García-Muelas, J. Morales-Vidal, M. Philipp, O. V. Safonova, N. López, J. A. Stewart, D. C. Ferré and J. Pérez-Ramírez, Nat. Commun., 2021, 12, 1–9.
59. D. Albani, M. Capdevila-Cortada, G. Vilé, S. Mitchell, O. Martin, N. López and J. Pérez-Ramírez, Angew. Chem., Int. Ed., 2017, 56, 10755–10760.
60. J. Ye, C. Liu, D. Mei and Q. Ge, ACS Catal., 2013, 3, 1296–1306.
61. J. Ye, C. Liu and Q. Ge, J. Phys. Chem. C, 2012, 116, 7817–7825.
62. J. Ye, C. J. Liu, D. Mei and Q. Ge, J. Catal., 2014, 317, 44–53.
63. M. S. Frei, M. Capdevila-Cortada, R. García-Muelas, C. Mondelli, N. López, J. A. Stewart, D. Curulla Ferré and J. Pérez-Ramírez, J. Catal., 2018, 361, 313–321.
64. X. Liu, L. Sun and W. Q. Deng, J. Phys. Chem. C, 2018, 122, 8306–8314.
65. J. Ye, C. Liu, D. Mei and Q. Ge, ACS Catal., 2013, 3, 1296–1306.
66. J. Ye, C. Liu, D. Mei and Q. Ge, ACS Catal., 2013, 3, 1296–1306.
67. M. Millet, G. Algara-Siller, S. Wrabetz, A. Mazheika, F. Girgsdies, D. Teschner, F. Seitz, A. Tarasov, S. Levchenko, R. Schlögl and E. Frei, J. Am. Chem. Soc., 2019, 6, 2451–2461.
68. C. Shen, Q. Bao, W. Xue, K. Sun, Z. Zhang, X. Jia, D. Mei and C. J. Liu, J. Energy Chem., 2022, 65, 623–629.

---

Footnote

† Electronic supplementary information (ESI) available: Description of the structure of SA-In₂O₃ models and oxygen vacancy formation (section S1–S3), geometries of elementary reaction steps (section S4–S7) full DFT dataset (section S8), and further microkinetic and DFT results (section S9 and S10). See DOI: https://doi.org/10.1039/d3cy00222e

---