A computational study of direct CO 2 hydrogenation to methanol on Pd surfaces †

The reaction mechanism of direct CO 2 hydrogenation to methanol is investigated in detail on Pd (111), (100) and (110) surfaces using density functional theory (DFT), supporting investigations into emergent Pd-based catalysts. Hydrogen adsorption and surface mobility are firstly considered, with high-coordination surface sites having the largest adsorption energy and being connected by diﬀusion channels with low energy barriers. Surface chemisorption of CO 2 , forming a partially charged CO 2 d (cid:2) , is weakly endothermic on a Pd (111) whilst slightly exothermic on Pd (100) and (110), with adsorption enthalpies of 0.09, (cid:2) 0.09 and (cid:2) 0.19 eV, respectively; the low stability of CO 2 d (cid:2) on the Pd (111) surface is attributed to negative charge accumulating on the surface Pd atoms that interact directly with the CO 2 d (cid:2) adsorbate. Detailed consideration for sequential hydrogenation of the CO 2 shows that HCOOH hydrogenation to H 2 COOH would be the rate determining step in the conversion to methanol, for all surfaces, with activation barriers of 1.41, 1.51, and 0.84 eV on Pd (111), (100) and (110) facets, respectively. The Pd (110) surface exhibits overall lower activation energies than the most studied Pd (111) and (100) surfaces, and therefore should be considered in more detail in future Pd catalytic studies.


Introduction
Methanol synthesis by direct hydrogenation of CO 2 has been recognised as a potential route towards sustainable fuels for transport and a circular fuel economy. 1 The industrial synthesis of methanol involves syngas, which is a mixture of CO/CO 2 /H 2 commonly produced from coal gasification.Whilst methanol synthesis from fossil fuels is efficient and profitable, environmental pressures are urging the chemical industry to transfer from a linearoil economy to net zero emissions by 2050. 2 Bussche et al. built a steady-state kinetic model that revealed CO 2 , and not CO, is likely to be the main source of carbon in methanol synthesised from syngas; [2][3][4] such knowledge encourages consideration of direct CO 2 hydrogenation to methanol, using anthropogenic CO 2 from the atmosphere.However, a better understanding of the interaction of CO 2 with transition metal catalysts is required for the design of novel and effective catalytic systems.Many factors, such as the source of H 2 , affect the extent to which the process of direct methanol synthesis from CO 2 can be ''green''; however, the idea of using an atmospheric pollutant such as CO 2 for fuel synthesis, and/or also generating feedstock for further synthesis of chemical compounds, such as formic acid, is broadly appealing. 5 crucial step in the direct hydrogenation of CO 2 to methanol is the initial CO 2 activation.On a heterogeneous catalyst, the reverse water-gas shift (RWGS) reaction needs to be inhibited while maintaining a strong interaction between CO 2 and the catalytic surface. 5,6][9][10][11] In order to understand fully the Pdbased alloy reactivity, it is necessary first to know the nature of the interactions between CO 2 and Pd.3][14][15] Burghaus et al. reported that CO 2 reactivity on clean Pd surfaces is weak, not favouring dissociation to CO unless an alkali metal species is coadsorbed. 15The weak interaction is considered to be predominantly a van der Waals physisorption, based on the theoretical and experimental observations at the Pd (111) surface. 14,16,17CO 2 adsorption on Pd has been studied in the context of the RWGS reaction and utilisation of syngas, and desorption of CO 2 from the Pd (111) surface is reported as requiring 0.26 eV of energy. 14Solymosi et al. reported that CO 2 desorption from the Pd (100) surface also has a relatively low energy of 0.35 eV, which was, in contrast, associated with a chemisorption, involving a metal to empty CO 2 p* orbital electron transfer. 12Evidence of CO 2 chemisorption on Pd (110) in the presence of water was also reported by Brosseau et al. 18 Therefore, the character of the CO 2 interaction with Pd surfaces seems to depend on the surface structure.4]18 Complementary to these observations, the rate of catalytic hydrogenation of CO 2 on Pd increases greatly when the active species is paired with suitable metal oxide supports, such as TiO 2 and ZnO, as they facilitate CO 2 adsorption and activation. 9,19,20Ko et al. computed the adsorption of CO 2 on transition metal surfaces, using the dispersion-corrected PBE-D2 density functional, and reported two types of CO 2 adsorption on Pd (111): an exothermic physisorption (À0.33 eV) of undistorted CO 2 , parallel to the surface; and a less exothermic chemisorption (À0.18 eV) with CO 2 in a bent geometry, and having a partial negative charge. 21In contrast, Zhang et al. recently calculated the CO 2 chemisorption on Pd (111) to be endothermic (0.06 eV) using the PBE density functional, in agreement with Habas et al., who reported the adsorption energy of CO 2 to be 0.22 eV above the dissociation limit, using DFT with the B3LYP density functional. 22Liu et al. have also shown that, when using the PBE density functional, the inclusion of the DFT-D2 correction dramatically changes the adsorption energy of chemisorbed species on the Pd (111) surface, from 0.30 eV to À0.18 eV. 17 Although there is no consensus on the matter of the endo-or exothermicity of CO 2 chemisorption on Pd surfaces, the reported values are generally small, which agrees with the experimental reports of a weak interaction between CO 2 and Pd surfaces.
Direct CO 2 hydrogenation to methanol is proposed to proceed via a surface formate intermediate (HCOO*, where * indicates an adsorbed species), with Medford et al. having shown that HCOO* could act as a poison for other reaction pathways due to its high thermodynamic stability on the catalyst surface. 23Variations of the mechanism proposed by Grabow, which progresses via formic acid (HCOOH) as shown in Fig. 1, have been presented, such as an initial Eley-Rideal type mechanism on Cu-based catalytic systems, where CO 2 in the gas phase reacts with surface-bound hydrogen to yield formate. 24 Recently, Hus ˇet al. concluded that dioxymethylene (H 2 COO*) should be preferably considered over formic acid (HCOOH*) on Cu-based catalysts as the former is more strongly bound to the metallic surface and the activation energy towards hydroxymethoxy (H 2 COOH*) is lower. 25d-based catalysts supported on ZnO are potent alternative catalysts for this reaction, with their reactivity attributed to the Pd-Zn binary metallic phases and their stabilisation of the HCOO* intermediate, similar to the Cu-based catalyst.9,20 Zhang et al. have reported DFT studies of an alternative CO 2 formate mechanism that involves dissociation of HCOOH to HCO and OH, and subsequent hydrogenations of HCO to produce CH 3 OH.27 Furthermore, Brix et al. have recently considered the initial CO 2 hydrogenation on Pd (111) to proceed via carboxylic acid (COOH), instead of formate, in a DFT study using the dispersion-corrected PBE-D3 density functional, and they observed a high energy barrier of 2.23 eV for CO 2 hydrogenation to formate on Pd (111), in contrast to the barrier of 0.85 eV reported by Zhang et al. 23,24 Whilst binary metallic alloy catalysts may offer more desirable selectivity, stability, and tunability than their monometallic counterparts, the lack of basic understanding of the behaviour of monometallic materials hinders the design of emergent multi-component materials.To achieve the required insight in the context of CO 2 hydrogenation over Pd, we need to understand reactivity across all the prominent surface facets.Thus, we present here an in-depth investigation of CO 2 interaction with low energy Pd (111), ( 110) and (100) surfaces using DFT calculations, followed by investigation of the direct CO 2 hydrogenation to CH 3 OH, via the Grabow mechanism, on the Pd (111), ( 110) and (100) surfaces, in the context of rationalising CO 2 reactivity on Pd-based catalysts.

Methodology
The Fritz Haber Institute ab initio molecular simulations (FHIaims) software package has been used for full potential allelectron DFT calculations, with the Pythonic Atomic Simulation Environment (ASE) used for management of calculation geometries. 28,29The default convergence criteria within FHIaims for self-consistent field (SCF) calculations were used, i.e. the changes between the current and previous SCF iterations in (via H 2 COO*, orange).* indicates a surface-bound species and d-indicates that CO 2 is partially charged (i.e.activated). 25,26harge density, sum of eigenvalues, and total energy were below N Â 1.67 Â 10 À5 e a 0 À3 , 10 À3 eV, and 10 À6 eV, respectively, where N is the number of atoms in the model.Scalar relativistic treatment of kinetic energy for all elements was achieved by the atomic zero-order regular approximation (ZORA), and a Gaussian-type broadening with width of 0.01 eV was applied to the occupation of electronic states.The Perdew-Burke-Ernzerhof exchange correlation (XC) density functional has been used unless explicitly stated otherwise, paired with the Tkatchenko-Scheffler van der Waals dispersion correction (PBE + vdW).A default ''light'' basis set (version: 2010) has been used for geometry optimisations, providing structural accuracy; 28,30,31 energy calculations were then performed with a ''tight'' basis set (version: 2010) on the optimised geometries, providing greater electronic accuracy and mitigation of basis set superposition error. 28For geometry optimisations, convergence was deemed complete when forces on all unconstrained atoms were less than 0.01 eV Å À1 .Due to the closed-shell electronic configuration of Pd ([Kr] 4d 10 ), spin-paired calculations were used in periodic calculations; gas-phase adsorbate structures were calculated both spin-paired and spin-unpaired, and the energy of the more stable configuration considered for reference in subsequent periodic calculations.The effect of the spin-paired approximation has been assessed towards the activation energies in relevant surface hydrogenation reactions in Section S3 of the ESI, † with a spin-paired treatment shown to introduce small error bars of AE0.05 eV.

Surface models
Using the optimised model of bulk FCC Pd, a surface supercell was created with dimensions of (3 Â 3 Â n), where n is the number of atomic layers in the z-direction perpendicular to the material surface.The x-and y-dimensions were chosen such that the adsorbates are significantly separated (7.5 Å), and a vacuum layer of 40 Å was added in the z-direction.The k-grid sampling was reduced appropriately for altered cell dimensions, with a k-grid of (3 Â 3 Â 1) applied.Due to the onesided nature of the slab models considered, a dipole-correction was used in all calculations.
The energy penalty for breaking chemical bonds at the surface of a material (E cleave ) is calculated as: where the DFT total energy of an unrelaxed surface slab model (E Unrelaxed Slab ), the bulk energy per atom (E bulk ), the number of atoms in the model (N), and the surface area (A), are needed.E cleave converges for the Pd (111), (100) and (110) facets when E cleave ceases to fluctuate as a function of slab thickness, as can be seen for n Z 5 in Fig. 2; thus, 5 layer models are used for all subsequent calculations.
To calculate the surface energy (E surf ), the energy of stabilisation provided by geometry relaxation (E relax ) needs to be obtained from the difference in total DFT energy of the optimised slab (E Relaxed

Slab
) and E Unrelaxed Slab : where a single-sided model of the surface is considered, hence the denominator is A only.Constraints were used to maintain the bulk structure for Pd atoms distant from the adsorption site, i.e. the bottom layers of the slab model.E Relaxed In summary, accurate results have been achieved herein with a 3 Â 3 Â 5 supercell surface model, with the bottom two layers constrained to their bulk positions and the three top surface layers unconstrained.The surface energies (E surf ) can subsequently be calculated as follows: with these settings, and are presented in Table 1.The calculated Pd (111), ( 100) and (110) E surf match previous computation and experiments, thus supporting the validity of our approach.

Surface adsorption
For catalytic surface reactions, the adsorption energy (E ads ) measures the interaction between a surface and reactant, and is deduced from comparison of the energies of the optimised gas-phase adsorbate (E A ), optimised surface (E S ) and the combined system (E A-S ).
where a negative value indicates favourable adsorption.Due to basis set incompleteness when using an atom centred basis, a Boys-Bernardi counterpoise correction is necessary for surfaceadsorbate interactions to account for the basis set superposition error (BSSE). 43In our work, the BSSE for CO 2 adsorbed on Pd (111) was assessed on an aperiodic model with all Pd atoms within 7.0 Å of the adsorbed CO 2 molecule included (i.e.all atoms within the distance of the atom-centred basis cut-off, including those in neighbour cells).The energy of the CO 2 in the presence and absence of Pd basis functions (E A(A-S) and E A(A) , respectively) were compared, and the equivalent comparison of the energy of the slab model in the presence and absence of the basis functions of the CO 2 adsorbate (E S(A-S) and E S(S) , respectively) was also performed. 43The BSSE energy (E BSSE ) was then calculated as: 43 A more negative E BSSE indicates a greater overbinding error; however, by subtracting E BSSE from E ads , the counterpoise corrected adsorption energy can be established (E CP ads ) as: With the ''light'' basis set, E BSSE is À0.08 eV for CO 2 on Pd (111), but E BSSE was reduced to À0.02 eV with the ''tight'' basis set.Considering the low BSSE with the ''tight'' basis, which is used subsequently throughout this work, the E BSSE contribution to E ads was deemed negligible and was not subsequently calculated for species other than CO 2 .

Transition state structures
For kinetic studies, we have used a machine learning nudged elastic band (MLNEB) method to identify saddle points and minimum energy paths (MEPs). 44,45A spring constant of 0.05 eV Å À1 has been used throughout; the convergence criterion of forces on all unconstrained atoms of below 0.05 eV Å À1 , with energy uncertainty below 0.03 eV, was deemed sufficiently accurate for CO 2 adsorption.Comparison with a more stringent force criterion of 0.01 eV Å À1 altered the activation energy for CO 2 adsorption on FCC (111) surface by 5 meV only (Section S2, ESI †).

Vibrations and Thermodynamics
Vibration frequency calculations have been performed for structures in the reaction pathway using the ASE Vibrations module and the Frederiksen method. 46To reduce the computational cost, the number of vibrating Pd atoms in the structures was tested as a function of distance from the adsorbate, from the nearest Pd neighbours outwards.The vibrational energy converged at first Pd neighbours of the adsorbate, and this approximation has been applied throughout.The vibration frequencies were also calculated for transition state structures, and the saddle points were confirmed by the presence of an imaginary frequency vibration.Due to the use of a finite difference approach, occasional additional imaginary frequencies did occur due to numerical noise, but were always below 10i meV and thus did not impact on the results.Free energies of gas components were calculated using an ideal gas-phase approximation, whereas the energies of periodic surface structures with adsorbates were calculated using a harmonic approximation.Once the zero-point energy (ZPE) was established, the contribution of ZPE was added to E ads yielding enthalpy of adsorption, H ads .

Hydrogen adsorption
Prior to investigating the reaction steps in CO 2 hydrogenation, an understanding of hydrogen behaviour on Pd is crucial as the interaction of H with CO 2 is integral to the reaction profile.Thus, a survey was conducted of E ads for H, E ads (H), on the Pd surfaces; the H atom was positioned at various locations on the surface and optimised, with constraints in the xy-plane.E ads (H) on the Pd (111), (100), and (110) surfaces was calculated with  41 Experiment 2.00 --Boer et al. 42 Experiment 2.01 --respect to the gas-phase diatomic hydrogen molecule and is plotted across the xy-plane in Fig. 3(i), 4(i), and 5(i), respectively.The ZPE-correction was not applied in these scans and all energies presented in this section are based on differences in total electronic energy.The most stable adsorption site for H on the Pd (111) surface is the HCP hollow position, which is site B in Fig. 3(ii), with E ads (H) of À0.67 eV; similar stability over high coordination sites is observed on the (100) surface, where the hollow site (Fig. 4(ii), site C) has lower energy with E ads (H) of À0.54 eV, and the ''FCC'' site (Fig. 5(ii), site B) is of lower energy in the case of the (110) surface, with E ads (H) of À0.56 eV.The least stable adsorption site for H on all surfaces is atop, with E ads (H) of À0.12, À0.08 and 0.00 eV on the (111), ( 100) and (110) Pd surfaces, respectively.The typical reactant feed for CO 2 hydrogenation is between 1 : 3 to 1 : 9 molar ratio of CO 2 and H 2 , and thus dissociated hydrogen would be readily available on the catalyst surface. 9,20High hydrogen mobility can be deduced from Fig. 3(i), 4(i) and 5(i), as differences in favourable E ads are low along specific channels, highlighted in red.110) surfaces, respectively. 47There is also good agreement with previous theoretical research, as Herron et al. calculated atomic hydrogen adsorption energies on Pd (111) using PW91 and reported them to be À0.59,À0.56, and 0.00 eV on FCC hollow, HCP hollow, and atop positions, respectively. 48Similarly, Fonseca et al. used PBE in their DFT study of hydrogen adsorption on Pd (111) and observed À0.66, À0.61, and À0.50 eV adsorption energies of hydrogen atom on FCC hollow, HCP hollow, and bridge position, respectively. 49

CO 2 adsorption
The adsorption energies and structures for CO 2 on the Pd (111), (100) and (110) surfaces are reported in Table 2.The undistorted CO 2 is most stable with a C-Pd bond distance, d(C-Pd), of 3.454 Å, which agrees with the physisorbed species observed by Habas et al. 22 H ads (CO 2 ) is strongest on the close-packed (111) surface, and is found to relate linearly with the number of Pd atoms that neighbour the surface adsorption site; when H ads (CO 2 ) is plotted as a function of surface atom coordination number, which are 9, 8, and 7 for the Pd (111), (100), and (110) surfaces, respectively, a linear fit returns R 2 = 0.998.
The stronger physisorption, rather than chemisorption, observed for CO 2 on the Pd (111) surface (E ads (CO 2 ) = À0.21eV) was reported previously by Ko et al. 22  ) is noted as increasingly negative (i.e.strengthens) with increasing E surf for the Pd facets, and the energy difference between surface-bound CO 2 and CO 2 dÀ also decreases; these observations agree with experimental data that show an absence of chemisorption on the Pd (111) surface, and both physisorption and chemisorption on the Pd (100) 3][14] Despite differences in H ads (CO 2 dÀ ) on the surfaces examined, the adsorbed geometries of CO 2 and CO 2 dÀ are consistent across all surfaces (Table 2); only a small difference in angles (0.41) is calculated for either the physisorbed or chemisorbed geometries when compared across the three facets.The impact of steric interactions for adsorbed CO 2 can be quantified via the distortion energy, i.e., the gas-phase energy of the adsorbed bent CO 2 geometry relative to the preferred linear CO 2 configuration, which is 1.35, 1.33 and 1.33 eV for Pd (111), Pd (110), and (100) facets, respectively.Given that the overall adsorption energies are exothermic on Pd (100) and ( 110) surface facets, it can be concluded that the   For CO 2 physisorption on the Pd (111) surface (Fig. 7a), the charge of the carbon (q C ) is +0.47 e, very similar to the gas phase CO 2 (q C = +0.48e), and only small changes are observed on the surface Pd.For CO 2 dÀ on the Pd (111) surface (Fig. 7b), negatively charged Pd atoms bond to an oxygen and carbon (q Pd 1 = À0.30e, q Pd 2 = À0.32 e).The distance d(C-Pd 2 ) is 2.85 Å, and there is a direct electronic interaction between Pd 2 and the carbon atom of CO 2 .The average charge on the second layer of Pd atoms, q Pdsublayer , decreases from +0.03 to +0.02 e upon physisorption and decreases further to À0.02 e upon chemisorption.The average charge on the first surface atomic layer of Pd, q Pdsurf , is À0.03 e, 0.00 e and À0.07 e for pristine Pd (111) surface, Pd (111) slab with CO 2 , and Pd (111) slab with CO 2 dÀ , respectively, suggesting that the electron density has been pulled to the first two layers of Pd, and to the CO 2 dÀ adsorbate via Pd 1 and Pd 2 .q C has decreased from +0.47 e to +0.39 e, indicating some metal (Pd 1 ) to empty CO 2 p* orbital electron transfer. 12The negatively charged oxygen close to the negative q Pd 1 and q Pd 2 will result in electrostatic repulsion, and thus are likely to contribute in the decreased stability of CO 2 dÀ on the Pd (111) surface. 22,51n contrast, for CO 2 dÀ on Pd (110) the q C reduction is by 0.07 e upon chemisorption, with higher electron density on the oxygens and much less charge redistribution on the surface Pd atoms compared to the Pd (111) surface (q Pd 1 of À0.10 e and À0.30 e, respectively), which all contribute to the overall stability (i.e.lower H ads ).(110) surfaces are more suitable for CO 2 activation than the most stable Pd (111) surface.

Interactions of intermediates with Pd surfaces.
Reaction intermediates from the Grabow mechanism, as introduced in Section 1.2, have been optimised on the pristine Pd ( 111), ( 100) and ( 110) surfaces, in each case starting from an atop position, which ensured that adsorbates were starting at a proximity allowing metal-adsorbate bond formation during the geometry optimisation process.For example, in the case of CO 2 , the chemisorbed species could easily be missed starting from the gas phase due to an energy barrier for the chemisorption of CO 2 on Pd ( 111) and (100) surfaces.
The calculated values of H ads are presented in Fig. 8 and tabulated in Section S4 of the ESI.† For the intermediates considered, the average difference between the highest and lowest H ads across the three surfaces is 0.22 eV; the smallest difference is for the CO 2 molecule (0.05 eV), and the largest for H 2 CO, H 2 COOH, and CO 2 dÀ (0.36, 0.33, and 0.29 eV, respectively).Plotting the surface energy (E surf ) of the low-index Pd surfaces against the adsorption enthalpy (H ads ) of these intermediates on the corresponding surfaces (Fig. 9) illustrates where surface properties associate with these observations.In particular, H ads of CO 2 , CO 2 dÀ , H 2 COOH, and H 2 CO present clear linear correlations with the stability of the surface facets, giving R 2 of 0.988, 0.997, 0.987, and 1.000, respectively.HCOO, HCOOH, H 3 CO, and CH 3 OH give a poor linear fit, which indicates that other factors, such as steric effects, should be considered for rationalising the strength of these adsorbate interactions with the Pd surfaces.For example, due to additional space on the long-bridge site on the Pd (110) surface, the HCOOH can be accommodated in a different orientation from that on the Pd ( 111) and ( 100) surfaces (i.e.C-H atoms facing down, rather than up), which makes the resulting structures more difficult to compare directly.

Transition states and reaction profile
To gain further insight into reaction mechanisms, activation energies were calculated.Here, for each reaction step that involved a hydrogenation, it was necessary to set the transition state (TS) starting geometry such that a hydrogen atom was positioned near to the intermediate; optimisation of these starting models with proximal hydrogen in some instances led to instability of the intermediate adsorption structures, which, however, did not cause problems when directly applied to the TS calculations.Once a transition state structure was confirmed, the ZPE-corrected activation energy (H act ) was calculated for each reaction step as: H act = H TS (species) À H ads (species).( 7) Fig. 8 H ads of the intermediates in the direct CO 2 hydrogenation to methanol, as studied on the low-index Pd surfaces, presented in order of increasing E surf : (111), ( 100) and (110), 26 in blue, orange and grey, respectively.Error bars of AE0.05 eV are provided to account for the spin-paired approximation applied to the adsorbed species, as described in Section 2.3.
This journal is © the Owner Societies 2022 Phys.Chem.Chem.Phys.
The reference initial state for the calculation of H act is the species adsorbed on the surface, and for hydrogenation steps a hydrogen atom is also adsorbed, but the two adsorbates have a very limited interaction between them.The resulting H act are presented in Fig. 10 and structures tabulated in Section S6 in the ESI.† Here, CO 2 dÀ was considered as the starting point, i.e.
proceeding via a Langmuir-Hinshelwood mechanism, and not a physisorbed CO 2 . 2626]55 All TS have also been validated by vibrational analysis, displaying only one imaginary frequency each.The elementary step towards which each energy barrier refers are presented in Table 4.Given that the adsorption energy of HCOOH on the Pd (111) surface is calculated as À0.58 eV, and desorption is considered as the reverse process, the high H act (HCOOH) observed (1.41 eV) for the Pd (111) surface suggests that HCOOH is more likely to desorb than react further.The high activation barrier for HCOOH hydrogenation agrees with work by Hus ˇet al. on Cu-based catalysts; however, formic acid is not amongst the product stream observed when using Pd catalysts experimentally, with CH 3 OH, CO and trace to significant amounts of CH 4 reported. 5,9,20,25,56Thus, another intermediate, such as H 2 COO, might be of importance in leading to the experimental products, as was determined for Cu-based catalysts. 25In our work, the H act (HCOOH) on the Pd (110) surface is about 40% lower than on Pd (111) and almost 45% lower than on the Pd (100).The reduction of H act (HCOOH) might stem from lower stability of the HCOOH, and reduced stability of the hydrogen atom on Pd (110), which translates into a more accessible transition state.
The activation enthalpy for dissociation of H 2 COOH species, H act (H 2 COOH), is highest on the Pd (100) surface, where the H 2 COOH intermediate is stabilised.Brix et al. reported a high H act (H 2 COOH) of 2.01 eV on Pd (111), while we calculate H act (H 2 COOH) to be only 0.40 eV; the significant discrepancy of 1.61 eV arises from a considerable difference in the transition state geometry, i.e., our transition state involves breaking of a single C-O bond, whereas both C-O bonds were broken in the transition state found by Brix et al.For hydrogenation of formaldehyde, H act (H 2 CO) is similarly low (0.67-0.74 eV) on the three surfaces; however, on the (111) surface it is higher than H ads of H 2 CO (À0.58 eV), whilst on Pd (100) and (110) surfaces, H 2 CO is stabilised more (À0.83and À0.94 eV) than on Pd (111).The stronger H ads on (100) and ( 110) surfaces means that H 2 CO desorption is less likely, and reactivity favoured, whilst desorption would be a competitive process on the (111) surface.Desorption of H 2 CO during CO 2 hydrogenation to methanol on Pd catalysts is a major concern in experiment, 5 and thus the Pd (100) and (110) surfaces may be preferable in catalyst design.
In most hydrogenation steps examined on the three Pd surfaces, the reaction pathway favoured migration of the hydrogen atom towards the least stable on-top site before bonding to Fig. 11 The ZPE-corrected energy profile of CO 2 hydrogenation to methanol, via the formate pathway, on Pd (111), (100), and (110) surfaces, plotted in blue, orange, and grey, respectively, relative to the energy of pristine surface and gas phase reactants. 26Energies of intermediate structures and transition state geometries have been stoichiometrically balanced with energies of gas phase reactants; * indicates surface bound species.
the intermediate.Therefore, the relative stability of the gen adsorption sites, as shown in Section 3.1, has a major impact on the H act for most hydrogenation reactions on the Pd (111), ( 100) and (110) surfaces.Reducing the difference in stability for hydrogen atoms on the possible surface sites might be an important factor in the design of catalysts for CO 2 hydrogenation to methanol, as it could lead to reduction of H act for species reacting on a Pd-based catalyst.
A reaction profile based on the energy of initial, TS and final geometries, relative to the energy of isolated Pd (111), (100) and (110) surfaces and gas-phase reactants, is plotted in Fig. 11, with each individual step balanced stoichiometrically by energies of gas-phase molecules.
Based on total electronic energy for all the surfaces, which is presented in the ESI, † the reaction energy for the conversion of CO 2 to methanol is exothermic (À1.05 eV) relative to gas phase reactants, which agrees with reaction energy (À1.17 eV) derived from atomization energies. 5,57,58The gas-phase reaction enthalpy presented above (À0.26eV) is underestimated by 0.26 eV with respect to experimental values reported in literature, 5 and the magnitude of the error in gas-phase energies of molecules is typical for gradient corrected functionals, such as PBE. 59The highest H act in the CO 2 hydrogenation reaction across the Pd (100), ( 111) and (110) surfaces is H act (HCOOH), with values of 1.51, 1.41, and 0.84 eV, respectively, which is therefore a likely rate determining step (RDS) for the reaction.However, in an experimental study by Aas et al., the decomposition of HCOOH on Pd (110) was shown to require 0.42 eV, which is only 50% of the H act (HCOOH) on Pd (110), and therefore much more likely. 60An important feature of the reaction energy profile is that TS1 remains endothermic on all three surfaces, with respect to the gas phase reactants, which would inevitably influence the rate of the reaction.All transition states on Pd (110) remain either below net zero energy of the reaction or significantly lower than Pd (111) and Pd (100) when above, which indicates that Pd (110) is the most active among the surfaces investigated here.As highlighted in Section 3.1, the hydrogen atoms are stabilised strongly on the Pd (111), (100), and (110) surfaces; the binding energies of intermediates with a neighbouring hydrogen atom do not vary significantly from the sum of binding energies of the adsorbates calculated separately, which suggests that the presence of hydrogen neither stabilises nor destabilises the intermediates at the low 1/9 ML coverage of hydrogen considered. 26However, presence of a hydrogen atom at the nearest neighbouring site to chemisorbed CO 2 was observed to result in CO 2 desorption during geometry optimisation, to form a linear physisorbed species, which may indicate a lower stability of chemisorbed CO 2 with increasing hydrogen ML coverage.Experimentally, the presence of H 2 appears to induce a larger CO 2 intake both at increased temperature and/or pressure, but this phenomenon has been linked to CO 2 dissociation. 61

Gibbs free energy analysis
Gibbs free energy changes (DG, eV) between the reaction steps a-v in Table 5 across Pd (111), ( 110) and (100) were calculated and are shown in Fig. 12 and 13.The pressure (p) used was 1 atm (1013125 Pa), which was applied to gas components.The temperatures (T) considered were 0 K (i.e., enthalpy), 300 K for ambient conditions, and 500 K as typical experimental conditions for CO 2 hydrogenation over Pd catalysts. 20Graphs of the effect of temperature are shown for Pd (111) in Fig. 12, with other surfaces presented in Section S7 in the ESI.† A comparison of results for 500 K presented in Fig. 13.
As T increases (Fig. 12), formation of species on Pd (111) from respective TS structure shows more negative (favourable) Fig. 12 The Gibbs free energy changes between reaction steps in CO 2 hydrogenation reaction via formate on Pd (111) at p of 1 atm and T of 0 K, 300 K and 500 K; reaction steps a-v are explained in Table 5.
DG for formation of H Large positive DG is observed for all reaction steps involving breaking of Pd-H bonds and attaching of the hydrogen to the adsorbates, i.e. formation of TS1 (d), TS2 (g), TS3 (j), TS5 (r) and TS6 (u), which indicates that the very strong Pd-H interaction at 1/9 ML hydrogen coverage impedes the reaction on Pd (111), (100) and ( 110) surfaces even at the first hydrogenation step to formate.
The Pd (110) surface has the lowest DG values for TS formations, except for dissociation of H 2 COOH (m), which is more favourable on the Pd (111) surface.The DG associated with TS formation are not strongly affected by T, which suggests that the flat low-index surfaces of Pd are not the likely source of methanol formation in supported metallic Pd catalysts.The conclusion is in agreement with experiment, showing that pure unsupported Pd does not produce methanol at 463 K and atmospheric pressure. 9Moreover, changing T was shown to have a very limited effect on formation of intermediates in CO 2 hydrogenation on Pd catalysts.Adsorbing hydrogen on Pd surfaces is less favourable at higher T, but formation of TS structures remains unaffected, thus rendering the process less feasible at high T; however, elevated T is necessary to activate CO 2 on Pd (100), showing that low T CO 2 activation is key for CO 2 hydrogenation to be kinetically viable.

Summary and conclusions
Direct hydrogenation of CO 2 to methanol on transition metal catalysts is a promising approach for green energy storage, and new and more efficient catalysts need to be designed to make the technology viable; here, we have investigated the CO 2 hydrogenation reaction via the formate pathway on Pd (111), (100), and (110) surfaces to support this endeavour.
Firstly, we investigated the stability of H on the Pd surfaces, showing that high coordination sites have the largest adsorption energy, and these high stability sites are interlinked via channels with low diffusion barriers; we also show for CO 2 adsorption that the preference of physical or chemical adsorption is dependent on the stability of the Pd surface facet.For the CO 2 hydrogenation reaction, the transition state for CO 2 dÀ hydrogenation (TS1), to form formate, is endothermic, which will influence the overall rate of the reaction.H act (HCOOH) is the highest energy step in the reaction profile on the Pd (111), (100), and (110) surfaces (TS3), and it can be considered as the likely rate determining step of this reaction on the surfaces examined.Based on thermodynamical analysis, the formation of TS1 has a high barrier that is not significantly influenced by reaction conditions, which indicates that flat surfaces of metallic Pd are unlikely to be the source of formate, and subsequently, methanol in product streams of reactions using supported monometallic Pd catalysts.Increased temperature was found to facilitate CO 2 chemisorption on Pd (100) and ( 110), but has an adverse effect on multiple other reaction steps, including the free energy of hydrogen adsorption.Novel Pd-based polymetallic nanoparticle catalysts for direct CO 2 hydrogenation to methanol could be designed to lower the barrier to initial CO 2 hydrogenation, TS1, and lower the barrier for formic acid hydrogenation (TS3) or facilitate a mechanism that proceeds via an alternative intermediate, such as H 2 COO.Importantly, a Pd-based CO 2 hydrogenation catalyst should have lower Pd-H binding strength to facilitate the reaction.
Overall, we show that the most stable geometry of CO 2 adsorbed on Pd surface varies across (111), (100) and (110) facets and future studies should not be limited to consideration of the lowest energy (111) surface facet.Inclusion of zero-point energy has shown the Pd (100) surface to be unsuitable for CO 2 hydrogenation, highlighting that consideration of enthalpy is important for accuracy in computational catalysis.Low-index Pd surfaces are therefore unlikely to be the source of methanol formation on supported monometallic Pd catalysts, which indicates the potential importance of low-coordination metallic sites and metal-support interfacial sites.[64]

Conflicts of interest
There are no conflicts of interest to declare.13 The Gibbs free energy changes between reaction steps in CO 2 hydrogenation reaction via formate on Pd (111), (110), and (100) at p of 1 atm and T of 500 K; reaction steps a-v are explained in Table 5.

Fig. 1
Fig. 1 Formate pathway of direct CO 2 hydrogenation to methanol on metallic surfaces, as proposed by Grabow et al. (via HCOOH*, blue) and Hus ˇet al.(via H 2 COO*, orange).* indicates a surface-bound species and d-indicates that CO 2 is partially charged (i.e.activated).25,26

Fig. 4
Fig. 4 (i) Adsorption energy (E ads ) of hydrogen atom on Pd (100) surface calculated as a function of x-and y-coordinate; the H atom remained constrained in the xy-plane during each geometry optimisation.A key is provided for the adsorption energies, in units of eV.(ii) Top-down view of the FCC Pd (100) surface with a 3 Â 3 Â 5 atoms simulation cell.Blue spheres represent Pd atoms and yellow circles represent unique adsorption sites: (a) atop, (b) bridge (c) hollow.Black lines represent the x-and y-direction cell boundaries.

Fig. 5
Fig. 5 (i) Adsorption energy (E ads ) of hydrogen atom on Pd (110) surface calculated as a function of x-and y-coordinate; the H atom remained constrained in the xy-plane during each geometry optimisation.A key is provided for the adsorption energies, in units of eV.(ii) Top-down view of the FCC Pd (110) surface with a 3 Â 3 Â 5 atoms simulation cell.Blue spheres represent Pd atoms and yellow circles represent unique adsorption sites: (a) short bridge, (b) ''FCC'' (c) atop, (d) long bridge, (e) hollow.Black lines represent the x-and y-direction cell boundaries.
For CO 2 dÀ on a Pd (100) facet, the charges calculated are intermediary between the results on the Pd (111) and (110) surfaces, and H ads also falls between the values observed for Pd (111) and (110) surfaces.The overall charge transfer from the metal to CO 2 dÀ is À0.04, À0.10 and À0.11 e on Pd (111), (100), and (110) surfaces, respectively, which is small but correlates directly with adsorption strength.In the previous literature, Bader charge analysis has been considered for CO 2 chemisorption on Pd (111) surfaces, and the transfer to CO 2 dÀ reported as À0.28 and À0.43 e by Tang et al. and Habas et al., respectively; 22,52 the direction of charge transfer is consistent with our own observations, with the quantitative difference attributed to methodological differences, i.e.Mulliken charge analysis has a stronger basis set dependency than Bader analysis. 53,54Importantly, we show qualitatively that the charge transfer to CO 2 dÀ increases over Pd (111), (100), and (110) surfaces, indicating that Pd (100) and

Fig. 6
Fig. 6 Side-and top-view of CO 2 chemisorbed on the FCC Pd (100) surface, illustrating notations used for Mulliken analysis.Blue, red, and grey spheres represent Pd, O, and C atoms, respectively.Black crosses mark constrained bulk Pd layers.

Fig. 7 A
Fig. 7 A red-white-blue (negative-neutral-positive charges) color-coded visualisation of the net Mulliken charge on atoms for (a) CO 2 physisorbed and (b) CO 2 dÀ chemisorbed on the Pd (111) surface.

Fig. 9 E
Fig. 9 E surf of the Pd (111), (100) and (110) surfaces plotted against H ads of intermediates in the mechanism of CO 2 hydrogenation to methanol.The red dashed line is the linear fit of the data points, and R 2 is the linear coefficient of determination showing the quality of the fit.

Fig. 10
Fig.10The ZPE-corrected activation energies (H act , eV) of reaction steps in the pathway for CO 2 hydrogenation to methanol, presented for Pd surfaces in order of increasing E surf , i.e. (111), (100), and (110), given in blue, orange, and grey, respectively.In each case, the label refers to the initial reaction species, and the activation energies are also tabulated.

Fig.
Fig.13The Gibbs free energy changes between reaction steps in CO 2 hydrogenation reaction via formate on Pd (111), (110), and (100) at p of 1 atm and T of 500 K; reaction steps a-v are explained in Table5.

Table 1
Pd FCC (111), (100) and (110) surface energies calculated using the outlined settings.Literature and experiment are provided for comparison 250.33 eV); they also identify a chemisorbed state CO 2dÀ with E ads = À0.16eV,17whichcompares with our observation of H ads (CO 2 dÀ ) = 0.09 eV.Similarly, Hus ˇet al. observed that on Cu catalysts, CO 2 binds to the metal surface in a bent geometry, where one of the oxygens binds to a secondary metal atom and the carbon binds to a metal atom underneath.25Highamet al. observed an endothermic CO 2 chemisorption on Cu (100) and (110) surfaces, similarly to our result on Pd (111).50 Energy differences between our results and those of Ko et al. are probably due to the choice of van der Waals correction; 21 never-the-less, the observed trends are very similar, and the stability of the physisorbed CO 2 implies that there is an energy barrier on the Pd (111) surface for the activation of CO 2 .H ads (CO 2 27) is endothermic (0.09 eV) on the Pd (111) surface, matching the work of Zhang et al.,27and exothermic (À0.09 and À0.19 eV) on the Pd (100) and (110) surfaces, respectively.27Reduction of the size of the model surface, such that 1/4 monolayer (ML) coverage of CO 2 is achieved on Pd (111), (100) and (110) surfaces, results in H ads (CO 2 dÀ ) of 0.12 eV, À0.03 eV and À0.16 eV, respectively.The higher (less favourable) H ads (CO 2 dÀ ) values for 1/4 ML coverage, when compared to the 1/9 ML coverage presented in Table 2, are intuitively linked to unfavourable interactions between neighbouring adsorbates.H ads (CO 2 dÀ

Table 2
Geometric and energetic observations for CO 2 and CO 2 Pd 1 ) is the distance between the carbon and the nearest neighbouring Pd given in Å, and + O-C-O is the angle between the oxygen, carbon, and oxygen, given in 1 dÀ physisorbed and chemisorbed species on low-index Pd surfaces, respectively; CO TS 2 is the transition state geometry between these stable local minima.H ads is the ZPE-corrected species adsorption energy, given in eV; d (C-adsorbed on the Pd (111), (110), and (100) facets provides insight into the electronic charge of the surface species, and the data acquired are reported inTable 3. The notation used for describing charges on atoms of interest is shown in Fig. 6: O 1 and O 2 are oxygen atoms on CO 2 molecule; the two closest Pd atoms interacting with CO 2 are labelled Pd 1 and Pd 2 , where Pd 1 is closest to O 1 and Pd 2 is closest to O 2 ; and Pd surf , Pd sublayer , and Pd slab refer to the first, second and all layers of Pd atoms in the model, respectively.

Table 3
Net Mulliken charges, in units of e, on relevant atoms for CO 2 physisorption and chemisorption on the Pd (111), (110) and (100) surfaces; the charges (q) over Pd atoms have been averaged in the first surface layer surface (Pd surf ) and sublayer (Pd sublayer ), and summated over the whole slab (Pd slab )

Table 4
Detailed reaction steps towards which the H act (species) abbreviations refer to in Fig.10