Study of the Rate-Determining Step of Rh Catalyzed CO 2 Reduction: Insight on the Hydrogen Assisted Molecular Dissociation

: In the context of climate change mitigation, CO 2 methanation is an important option for the production of synthetic carbon-neutral fuels and for atmospheric CO 2 recycling. While being highly exothermic, this reaction is kinetically unfavorable, requiring a catalyst to be efﬁciently activated. Recently Rh nanoparticles gained attention as effective photocatalyst, but the rate-determining step of this reaction on Rh surface has not been characterized yet. In this work, Density Functional Theory and Nudged Elastic Band calculations were performed to study the Rh-catalyzed rate-determining step of the CO 2 methanation, which concerns the hydrogen assisted cleavage of the CO* molecule and subsequent formation of CH* and O* (* marks adsorbed species), passing through the CHO* key intermediate. The conﬁgurations of the various adsorbates on the Rh (100) surface were investigated and the reaction mechanism was studied exploiting different exchange-correlation functionals (PBE, RPBE) and the PBE+U technique. The methanation rate-determining step consists of two subprocesses which subsequently generate and dissociate the CHO* species. The energetics and the dynamics of such processes are extensively studied and described. Interestingly, PBE and PBE+U calculated activation barriers are in good agreement with the available experimental data, while RPBE largely overestimate the CHO* dissociation barrier.


Introduction
The last Intergovernmental Panel on Climate Change (IPCC) summary report addressed the prominent impact that human industrial activities are having on global climate conditions [1]. The enormous quantity of greenhouse gases that have been released into the atmosphere since the mid of the 20th century are causing a sharp increase in the global average temperature, with dramatic consequences on many ecosystems [1,2]. Reducing carbon emissions and mitigating the impact of climate change is thus a challenging and fundamental task. In particular, actions leading to the replacement of fossil fuels with other environmentally sustainable sources and to the decrease of CO 2 atmospheric concentration are especially needed. The effective production of synthetic carbon-neutral fuels from greenhouse gases could help addressing both these issues, making it possible at the same time to recycle atmospheric pollutants while providing a new source of fuels. In this context, CO 2 methanation represents a possible route toward the production of carbon-neutral fuels [3,4]. This process is represented by the following reaction: CO 2 + 4H 2 → CH 4 + 2H 2 O ∆H 0 = −165 kJ/mol (298 • K) (1) separated CH* and O* species. In the second step, O* diffuses on the surface reaching its equilibrium configuration, overcoming another lower energy barrier and relaxing the stress accumulated in the distortion of Rh outermost surface layer. While all methodologies and functionals agree on the general dynamics of the process, when using the RPBE functional CH* and O* stabilities are lowered by c.a. 0.35 eV with respect to PBE-based methods, and the CHO* dissociation energy barrier is greatly overestimated compared to experimental data. We relate the divergence among the methods to the different treatment of the interface region in the RPBE functional [34], which yields lower adsorption energies, a tendency that is stronger for the products than for the reactant. The paper is organized as follows: first the computational details of the work are provided, then the results of the geometry optimizations and the NEB simulations are presented and discussed, finally conclusions are drawn.

Clean Rh (100) Slab Characterization
To determine the geometrical structure of the catalyst and to compare the performance of the different xc functionals, the physical and geometrical features of the clean Rh (100) surface were studied. In particular, the slabs were characterized on the basis of their surface energies σ, work function ϕ, and surface lattice relaxation parameter d ij %. These quantities were computed on the basis of the definition reported in Appendix A as Equations (A1)-(A3). The calculated data, together with theoretical and experimental results available in the literature, are collected in Table 1. Among all the available computational data, only those obtained with computational setups similar to the ones adopted in this work are reported here. As shown in Table 1, our PBE results are in good agreement with previous calculations. Differences on the absolute values of the analyzed parameters are very small and can be ascribed to the different codes used by the authors [36,37]. The difference of about 0.17 eV/atom with respect to the surface energy computed in ref. [36] can be related to the different method adopted in the surface energy calculation. Indeed, in ref. [36], the surface energy was obtained through the linear regression on the total energy of the slabs as a function of the number of layers of the slabs themselves. Moreover, as can be seen from Table 1, PBE results agree well with experimental data apart from the relaxation of the outermost layer, which is overestimated by all the calculations, except for the RPBE one. It is worth noting that the 0.19 eV difference with the experimental work function is well within the range of the accuracy of DFT-PBE for this quantity [41]. Regarding RPBE results, all values are comparable with the other theoretical estimations, with the exception of d 1,2 % , which interestingly is closer to the experimental data. This could suggest a better performance of RPBE in determining the physical properties of the Rh (100) slab. However, as mentioned in the computational details, RPBE predicts a larger (and less accurate) Rh bulk unit cell parameter. This is reflected in the larger denominator in Equation (A3) and the smaller relative surface layer relaxation. Globally, our analysis suggests that the two XC functionals give comparable results in determining the clean slab physical properties. The major difference resides in the accuracy of determining the bulk Rh-Rh bond length which is more accurate using PBE rather than RPBE.

Description of the Adsorbates
As mentioned in Section 1, the key molecular species involved in the rds of CO 2 reduction are CO*, H*, CHO*, CH* and O*. We thus studied the equilibrium configuration of these species and the energetics of the various systems. Three possible adsorption sites, Bridge, Top and Hollow were considered. A visual representation of these sites on the Rh (100) surface is shown in Figure 1. For each adsorbate, the most stable configuration was identified, and the adsorption energies were calculated using the following expression: where E slab+spec is the total energy of the system composed on the slab and the adsorbed specie, E slab is the energy of the slab itself, and E spec is the energy of the free species. It is important to underline that given the 3 × 3 supercell dimension adopted for these calculations, our results refer to low (Θ = 1/9 ML) adsorbate coverage of the surface. The adsorption energies of the investigated molecular species are collected in Table 2, while the most stable geometrical configurations of the different species, calculated at the PBE level, are shown in Figure 2.
data. This could suggest a better performance of RPBE in determining the physical properties of the Rh (100) slab. However, as mentioned in the computational details, RPBE predicts a larger (and less accurate) Rh bulk unit cell parameter. This is reflected in the larger denominator in Equation (A3) and the smaller relative surface layer relaxation. Globally, our analysis suggests that the two XC functionals give comparable results in determining the clean slab physical properties. The major difference resides in the accuracy of determining the bulk Rh-Rh bond length which is more accurate using PBE rather than RPBE.

Description of the Adsorbates
As mentioned in Section 1, the key molecular species involved in the rds of CO2 reduction are CO*, H*, CHO*, CH* and O*. We thus studied the equilibrium configuration of these species and the energetics of the various systems. Three possible adsorption sites, Bridge, Top and Hollow were considered. A visual representation of these sites on the Rh (100) surface is shown in Figure 1. For each adsorbate, the most stable configuration was identified, and the adsorption energies were calculated using the following expression: where is the total energy of the system composed on the slab and the adsorbed specie, is the energy of the slab itself, and is the energy of the free species. It is important to underline that given the 3 × 3 supercell dimension adopted for these calculations, our results refer to low (Θ = 1/9 ML) adsorbate coverage of the surface. The adsorption energies of the investigated molecular species are collected in Table 2, while the most stable geometrical configurations of the different species, calculated at the PBE level, are shown in Figure 2. A first fundamental step of the methanation reaction is the generation of CO*. This adsorbed species has to satisfy two main requirements in order to represent a reliable intermediate for the reaction mechanisms namely thermodynamic stability and geometrical accessibility. The first requirement ensures that the molecular intermediate is more stable when adsorbed on the catalyst surface compared to the gas phase form, ensuring the stability of such species which would otherwise desorb, preventing further reaction steps. On the other hand, the geometrical accessibility is mandatory to allow interactions with co-adsorbed hydrogen atoms and thus to proceed towards the reaction products. As visible from the data collected in Table 2 and the picture in Figure 2 panel A, both these criteria are satisfied. Regardless of the calculation level, the adsorbed configurations are always more stable than the isolated systems, indicating that the adsorption of CO is an exothermic (and thus favorable) process. Moreover, within the relaxed geometry, the CO molecule is adsorbed perpendicularly to the surface, with the carbon interacting with the rhodium atoms and the oxygen pointing in the opposite direction. This suggests that the C atom can easily interact with hydrogen, because no steric hindrance is present along the directions parallel to the surface plane. From our analysis, the most stable configuration is obtained when CO is adsorbed on a Bridge site. The adsorption energies are close to the ones of the Top adsorbed geometries (maximum energy difference of about 70 meV), and A first fundamental step of the methanation reaction is the generation of CO*. This adsorbed species has to satisfy two main requirements in order to represent a reliable intermediate for the reaction mechanisms namely thermodynamic stability and geometrical accessibility. The first requirement ensures that the molecular intermediate is more stable when adsorbed on the catalyst surface compared to the gas phase form, ensuring the stability of such species which would otherwise desorb, preventing further reaction steps. On the other hand, the geometrical accessibility is mandatory to allow interactions with co-adsorbed hydrogen atoms and thus to proceed towards the reaction products. As visible from the data collected in Table 2 and the picture in Figure 2 panel A, both these criteria are satisfied. Regardless of the calculation level, the adsorbed configurations are always more stable than the isolated systems, indicating that the adsorption of CO is an exothermic (and thus favorable) process. Moreover, within the relaxed geometry, the CO molecule is adsorbed perpendicularly to the surface, with the carbon interacting with the rhodium atoms and the oxygen pointing in the opposite direction. This suggests that the C atom can easily interact with hydrogen, because no steric hindrance is present along the directions parallel to the surface plane. From our analysis, the most stable configuration is obtained when CO is adsorbed on a Bridge site. The adsorption energies are close to the ones of the Top adsorbed geometries (maximum energy difference of about 70 meV), and this justifies the co-existence of these two configurations observed experimentally [42][43][44]. Despite the qualitative agreement between theory and experiments concerning the adsorption geometries, experimental values for CO adsorption energy range between 1.35 ÷ 1.45 eV [43][44][45][46][47], while our DFT PBE calculations predict a much larger adsorption energy (ca. 2 eV), in agreement with previous theoretical results [48][49][50][51]. Indeed, despite PBE properly predicts properties like Rh bonds length and electronic structures [52], it was shown that it overestimates the interaction among the Rh d-type orbitals and the 2π* level of CO, leading to an overestimation of the CO adsorption energy on rhodium surfaces [53,54]. To overcome this problem, on one hand, we employed the PBE+U method, which allowed us to properly shift the CO 2π* level in order to obtain more accurate results [55,56]; on the other hand, we used the RPBE functional, which was specifically designed for yielding more accurate adsorption energies compared to PBE functional [34] and has been demonstrated to give more reliable results in modeling the CO-metal interactions [52]. As shown in Table 2, the two methods yield a much better agreement with the experimental adsorption energy values. Interestingly, these results are not associated with significant geometrical changes; indeed, the computed C-O bond length for the Bridge adsorption sites are 1.17 Å and 1.18 Å for the PBE-based and RPBE calculations respectively, while the distance between C and Rh atoms resulted in be 2.02 Å and 2.07 Å for the PBE-based and RPBE calculations, respectively. All these data are coherent with previous calculations [49,51,52].  The other reactant that must adsorb on the Rh surface is hydrogen. Previous investigations have shown that H 2 spontaneously dissociates and diffuses on the catalyst surface as soon as it is chemisorbed [57,59]. The adsorption energies reported in Table 2 indicate that the most stable adsorption sites for H* are Bridge and Hollow (the most stable geometry is shown in Figure 2 panel (B). These configurations are practically isoenergetic, as already reported in previous studies [60,61]  Given the atomic nature of the adsorbed species, the values reported in Table 2, computed with respect to the total energy of the isolated hydrogen atom, cannot be directly compared with experimental data which refer to the H2 dissociative chemisorption. To be able to compare directly with such experimental data and check the reliability of our results, we repeated the calculations for a slab with an H atom adsorbed on each side, this time computing the adsorption energy with respect to the total energy of the isolated H2 molecule; the results are listed in parenthesis in Table 2. Additionally in this case, the favorite chemisorption sites are Bridge and Hollow with an adsorption energy calculated with PBE functional of about 0.52 eV, in line with what has previously been observed [57,60]. Using the RPBE functional gives an adsorption energy of about 0.65 eV, slightly larger than the PBE case, and in worse agreement with the experimental values. Given the atomic nature of the adsorbed species, the values reported in Table 2, computed with respect to the total energy of the isolated hydrogen atom, cannot be directly compared with experimental data which refer to the H 2 dissociative chemisorption. To be able to compare directly with such experimental data and check the reliability of our results, we repeated the calculations for a slab with an H atom adsorbed on each side, this time computing the adsorption energy with respect to the total energy of the isolated H 2 molecule; the results are listed in parenthesis in Table 2. Additionally in this case, the favorite chemisorption sites are Bridge and Hollow with an adsorption energy calculated with PBE functional of about 0.52 eV, in line with what has previously been observed [57,60]. Using the RPBE functional gives an adsorption energy of about 0.65 eV, slightly larger than the PBE case, and in worse agreement with the experimental values.
The stability of the adsorbed CHO* complex, involved in the rds, compared to its gasphase counterpart, is confirmed by the adsorption energies reported in Table 2. The most stable geometry, reported in Figure 2 panel C, was obtained starting from both Bridge and Hollow configurations, using all three computational methods. As shown in Figure 2 panel C, the C-O bond is basically parallel to the surface, with a C-O distance, C-H bond length and an HCO angle of about 1.35 Å, 1.11 Å and 113 • , respectively, with these being the same for PBE, PBE+U and RPBE calculations. Such adsorbate geometry can be viewed as a Hollow configuration since the atoms' coordinations are comparable to those of a fourfold Hollow site, as visible from Figure S2 in the Supplementary Materials. This configuration represents a reliable starting point to study the C-O bond cleavage, which can qualitatively occur without requiring a large amount of energy, because the proximity of the Rh surface atoms can partially stabilize the highly energetical intermediates that are generated during the bond dissociation. Since this molecular fragment is not stable in the gas phase (it would exist only as a radical species), no comparison with experimental data could be carried out. However, calculations performed on similar systems revealed the same geometrical configuration and an adsorption energy on the Rh (100) surface of about 3.62 eV, which is compatible with our results [30]. As shown in Table 2, the PBE-based and RPBE CHO* adsorption energies are very different, the latter being 0.75 eV less stable than the former. This effect is related to the way the molecule-slab interface regions are treated by the functionals [34,62] and is even more pronounced for the CH* and O* systems, which are the products of the CHO* dissociation. As presented in Table 2, CH* adsorption energy changes of about 1.18 eV changing the functional from PBE to RPBE. This indicates that the choice of the functional can affect the calculated stability of the adsorbed species, and this is reflected in the results of the NEB calculations that we will discuss in the next section. Despite these differences, all calculations provide as equilibrium structure for the CH* species the Hollow adsorbed configuration, with the C-H bond directed perpendicularly to the Rh surface, as visible form Figure 2 panel D. All other adsorption sites give unstable structures that naturally relax into the Hollow configuration during the geometry optimization. The prominent stability of this configuration compared to the other adsorption geometries can be related to the undercoordination of the C atom in the CH species, which needs to be saturated by the highest number of surface metal atoms, occurring in Hollow sites (4-fold coordination). Due to the intrinsic instability of the CH species, in the gas phase, it is not possible to compare the adsorption energies with experimental data. Moreover, most of the available computational results on the CH* specie refer to different Rh surfaces [63][64][65] or to systems where CH* is co-adsorbed with hydrogen atoms [61]. Despite these differences, these studies in any case indicate the Hollow configuration to be the most stable one and report geometrical parameters comparable with the computed C-H bond length of 1.10 Å and 1.11 Å for the PBE-based and RPBE calculations respectively, and C-Rh distance of 2.09 Å and 2.11 Å for the PBE-based and RPBE calculations respectively.
Additionally, in the case of O*, adsorption energies are strongly dependent on the choice of the XC functional. Here RPBE underbinds the O* species with respect to PBE of about 0.95 eV for the Hollow adsorption site and of 1.05 eV for the Bridge configuration. As happened in the previous cases, these adsorption energy differences seem to be uncorrelated with the adsorption sites stability trend and almost identical geometrical configurations were obtained in all the three calculations. The most stable adsorbed configuration for O* is on the Hollow site, as reported in previous theoretical studies [66][67][68]. The Bridge-to-Hollow energy difference ranges between 0.12 and 0.22 eV, depending on the calculation level. As reported in Figure 2 panel E, the O-Rh distance is 2.15 Å for the PBE-based calculations, while it increases to 2.20 Å for the RPBE calculation. This bond length difference is compatible with the lower affinity between oxygen and the rhodium surface using the RPBE functional compared to the PBE-based schemes. To be able to compare with the experimental data, as in the H* case, we repeated the calculations referring to O 2 dissociative chemisorption. As visible from the data in parenthesis in Table 2, the favorite chemisorption site is still the Hollow spot with an adsorption energy calculated with PBE functional of about 4.53 eV, which overestimates the available experimental data for dissociative chemisorption of about 0.53 eV [58,69]. On the other hand, changing the functional from PBE to RPBE decreases the adsorption energy to 3.04 eV, and strongly underestimates the O 2 dissociative chemisorption energy (by 0.96 eV) with respect to the experiment.
Besides considerations on single molecule adsorption, we searched for the most stable co-adsorbate configurations for both the CO*/H* and CH*/O* pairs, which represent the initial and final stage of all our processes, respectively. Our analysis revealed that the less energetic conformation in the case of the CO*/H* pair can be obtained when these species are co-adsorbed in adjacent Bridge sites, as shown in the Supporting Movie (SM) and in Figure S2 attached as Supplementary Material. Regarding the CH*/O* pair, we found that the most stable configuration is the one where the two products species are located on the opposite corner of the square identified by four Hollow sites, as shown in the Supplementary Materials. This arrangement represents the target state of our reaction energy profile analysis.

Study of the Methanation Rate-Determining Step
As already mentioned, the rate-determining step of the Rh catalyzed CO 2 methanation reaction was identified in the CO* + H* → CH* + O* step through the formation of CHO* species [5,6,10,13,14]. Previous DFT simulations concerning the study of Ru (0001) catalyzed CO 2 reduction, showed that the CHO* formation and its dissociation have comparable energy barriers [13,29]. A kinetic study performed on top of such calculations points to the CHO* dissociation as the rds of the whole reduction process [13]. To the best of our knowledge, no theoretical studies are currently available concerning the CO 2 reduction over rhodium surfaces; what is known from the study of CH 3 OH dehydrogenation on Rh (111) is that the CHO* formation has an energy barrier similar to that found on the Ru (0001) surface [13,29,70]. Therefore, we calculated the energy landscape of the whole process generating the CH*/O* co-adsorbed pair from CO* and H* through CHO*, at the PBE, PBE+U and RPBE levels. We will discuss in the next paragraph the effect of the different functionals, while we focus here on the reaction mechanism itself, studied at the PBE+U level. The calculated energy profile is shown in Figure 3. chemisorption site is still the Hollow spot with an adsorption energy calculated with PBE functional of about 4.53 eV, which overestimates the available experimental data for dissociative chemisorption of about 0.53 eV [58,69]. On the other hand, changing the functional from PBE to RPBE decreases the adsorption energy to 3.04 eV, and strongly underestimates the O2 dissociative chemisorption energy (by 0.96 eV) with respect to the experiment.
Besides considerations on single molecule adsorption, we searched for the most stable co-adsorbate configurations for both the CO*/H* and CH*/O* pairs, which represent the initial and final stage of all our processes, respectively. Our analysis revealed that the less energetic conformation in the case of the CO*/H* pair can be obtained when these species are co-adsorbed in adjacent Bridge sites, as shown in the Supporting Movie (SM) and in Figure S2 attached as Supplementary Material. Regarding the CH*/O* pair, we found that the most stable configuration is the one where the two products species are located on the opposite corner of the square identified by four Hollow sites, as shown in the Supplementary Materials. This arrangement represents the target state of our reaction energy profile analysis.

Study of the Methanation Rate-Determining Step
As already mentioned, the rate-determining step of the Rh catalyzed CO2 methanation reaction was identified in the CO* + H* → CH* + O* step through the formation of CHO* species [5,6,10,13,14]. Previous DFT simulations concerning the study of Ru (0001) catalyzed CO2 reduction, showed that the CHO* formation and its dissociation have comparable energy barriers [13,29]. A kinetic study performed on top of such calculations points to the CHO* dissociation as the rds of the whole reduction process [13]. To the best of our knowledge, no theoretical studies are currently available concerning the CO2 reduction over rhodium surfaces; what is known from the study of CH3OH dehydrogenation on Rh (111) is that the CHO* formation has an energy barrier similar to that found on the Ru (0001) surface [13,29,70]. Therefore, we calculated the energy landscape of the whole process generating the CH*/O* co-adsorbed pair from CO* and H* through CHO*, at the PBE, PBE+U and RPBE levels. We will discuss in the next paragraph the effect of the different functionals, while we focus here on the reaction mechanism itself, studied at the PBE+U level. The calculated energy profile is shown in Figure 3.  As noticeable, all processes are endothermic and require a certain amount of energy to occur. The green line in Figure 3 separates the two subprocesses in which we split the calculations: the CHO* generation starting from CO* and H*, occurring in a single-step, and the CHO* dissociation, occurring in two sub-steps. During the CHO* formation, the species moves and rotate on the Rh (100) surface in order to generate CHO* and the Rh outermost layer gets slightly distorted. This step requires ca. 1.17 eV to take place, while the reaction energy is around 0.67 eV, as written in Figure 3. Such results compare well with simulations performed on Rh (111) and Ru (0001) surfaces (although different XC-functionals were used) [29,70]. A movie representing this first reaction step is available as SM, while the meaningful snapshots of the process are reported in Figure S2 in the Supplementary Material section. Differently from the simple landscape of this first process, the CHO* dissociation exhibits a more complex energy profile. We focus now on this process whose PB energy profile and related geometries are shown in Figure 4. The initial configuration of the dissociation is shown in Figure 2 panel C (side view) and in the inset A of Figure 4 (top view), while the final configuration is shown in inset C of Figure 4. The values of the energy profiles at the most relevant points of the reaction are explicitly indicated in Figure 4. As noticeable, all processes are endothermic and require a certain amount of energy to occur. The green line in Figure 3 separates the two subprocesses in which we split the calculations: the CHO* generation starting from CO* and H*, occurring in a single-step, and the CHO* dissociation, occurring in two sub-steps. During the CHO* formation, the species moves and rotate on the Rh (100) surface in order to generate CHO* and the Rh outermost layer gets slightly distorted. This step requires ca. 1.17 eV to take place, while the reaction energy is around 0.67 eV, as written in Figure 3. Such results compare well with simulations performed on Rh (111) and Ru (0001) surfaces (although different XCfunctionals were used) [29,70]. A movie representing this first reaction step is available as SM, while the meaningful snapshots of the process are reported in Figure S2 in the Supplementary Material section. Differently from the simple landscape of this first process, the CHO* dissociation exhibits a more complex energy profile. We focus now on this process whose PB energy profile and related geometries are shown in Figure 4. The initial configuration of the dissociation is shown in Figure 2 panel C (side view) and in the inset A of Figure 4 (top view), while the final configuration is shown in inset C of Figure 4. The values of the energy profiles at the most relevant points of the reaction are explicitly indicated in Figure 4. Following the energy profiles, the whole reaction can be roughly divided into two sub-processes: the first one comprises the C-O bond dissociation and O* migration reaching a first local energy minimum corresponding to the configuration B in Figure 4; in the following step, the O* further migrates to reach the products most stable configuration (inset C). In the first subprocess, starting from the CHO* equilibrium configuration (inset Following the energy profiles, the whole reaction can be roughly divided into two subprocesses: the first one comprises the C-O bond dissociation and O* migration reaching a first local energy minimum corresponding to the configuration B in Figure 4; in the following step, the O* further migrates to reach the products most stable configuration (inset C). In the first subprocess, starting from the CHO* equilibrium configuration (inset A in Figure 4), there is an initial energy growth, corresponding to the elongation of the C-O bond and consequent adjustment of the outermost Rh layer, which actively takes part in this first dissociation step. After this initial growth, the profiles reach their absolute maximums, which correspond to a state where the CH* and O* species are separated but the system is highly unstable, as the two fragments are located into the same 2 × 2 Rh cell, in mixed Bridge-Hollow sites. This configuration constitutes the transition state for the CHO* bond dissociation, as confirmed by the vibrational analysis which reveals the presence of an imaginary frequency associated with a normal mode pointing along the reaction coordinate. After this sharp maximum, the profiles exhibit a shoulder, corresponding to a configuration in which both species are in Bridge adsorption sites. This is a metastable configuration where the adsorbates are located in a relatively stable point but are not in their minimum energy configuration. In the following steps O* starts to migrate towards the Hollow site of the adjacent Rh 2 × 2 cell, while the CH* makes structural adjustments and falls into its energy minimum structure, a Hollow adsorption site, as previously discussed. Within this subprocess, the outermost Rh layer plays a primary role in supporting all structural changes. This can be noticed for example by looking at the inset B of Figure 4, where despite the whole system being in an energy minimum configuration, the Rh atoms of the outermost layer are still distorted, due to the long-range interactions acting among the two adsorbed species. In the second subprocess of the reaction, O* diffuses towards the Hollow site corresponding to its position in the target state (inset C of Figure 4), passing through a Bridge configuration. Once the final Hollow spot is reached, the interactions among the adsorbed species are minimized, the surface atoms can relax all the tension accumulated along these processes and the system reaches its global minimum. These two subprocesses can be viewed in the movie SM attached as Supplementary Material, while in Figure S2 the snapshots of the molecular configurations in the main points of interest are collected.
As shown in Figure 4, the process that requires more energy is the first dissociation of the C-O bond, with an activation energy barrier of about 1.12 eV to occur. These activation energy barriers are compatible with those of theoretical studies carried out on similar systems: in particular, Zhang et al. [29] estimated that CHO* dissociation can take place on Ru (100) with an activation barrier of about 1.48 eV, while Avanesian et al. [13] showed that on Rh (211), this step can occur with a free energy barrier of about 1.43 eV. At the same time, the computed barriers overestimate experimental values for the whole CO 2 methanation reaction that range from 0.7 to 1.0 eV depending on the support [16,19,[71][72][73]. The energy barrier overestimation with respect to experiments can be linked to contributions arising from the presence of co-adsorbed species. Indeed, the calculated energy profiles refer to the dissociation of an isolate CHO* species without accounting for the presence of other co-absorbed species. However, it was proved that after the first dissociative adsorption of CO 2 , the dissociated oxygen atoms remain absorbed on the surface and their presence reduces the activation barrier of the whole methanation reaction [19,66,73]. This explains at the same time why the methanation activation energy is lower for CO 2 than that of CO of about 0.15 ÷ 0.35 eV, depending on the experimental conditions [16,19,21,73], and the overestimation of the reaction barrier presented in this work.
Considering the high similarities in the physico-chemical features between the systems, our calculations suggests that the Rh (100) catalyzed process occur as in the case of Ru (0001), being the energy landscapes profiles and the activation barriers comparable for CHO* formation and dissociation [13,29]. Taking into account the role of kinetics, it was found that on the Ru (0001), the slowest process is the CHO* dissociation, even if this step is less energetically demanding compared to the CHO* generation [13]. In analogy with these results, the CHO* cleavage could be the rds of the whole CO 2 reduction also in this case, although further investigation on the kinetics of this reaction exploiting, e.g., variational transition state theory [74], are needed to fully validate this hypothesis.
Finally, all the results presented here refer to systems in zero-absolute temperature conditions and do not include any entropic contribution. Indeed, previous studies concerning this reaction on similar systems showed entropic corrections to the free energy on the order of 0.1 eV at room temperature [13,28,29], smaller than the difference related to the use of various computational setups, as shown in the next section.

Comparison among Different Computational Techniques
To highlight the different performances given by the various computational techniques, we report in Figure 5 the calculated NEB profile for the whole process, where for a better comparison, the three profiles are aligned by setting the energy of the initial configuration to zero.
Finally, all the results presented here refer to systems in zero-absolute temperature conditions and do not include any entropic contribution. Indeed, previous studies concerning this reaction on similar systems showed entropic corrections to the free energy on the order of 0.1 eV at room temperature [13,28,29], smaller than the difference related to the use of various computational setups, as shown in the next section.

Comparison among Different Computational Techniques
To highlight the different performances given by the various computational techniques, we report in Figure 5 the calculated NEB profile for the whole process, where for a better comparison, the three profiles are aligned by setting the energy of the initial configuration to zero. PBE, PBE+U and RPBE yield approximately the same results concerning the generation of CHO*. The activation energy is estimated in ca. 1.20 eV, while the reaction energy is around 0.75 eV for PBE and RPBE calculations, and 0.67 eV for the PBE+U one, as reported in Figure 5. At odds with these results, RPBE description of the CHO* dissociation step deviates strongly from the PBE and PBE+U ones. As is clearly visible from Figure 6, starting from a certain point along the reaction coordinate prior to the first maximum, the RPBE energy diagram is shifted to higher energies when compared to the PBE and PBE+U cases. This upshift slightly varies with the reaction coordinate until the minimum corresponding to configuration B is reached, while it is almost constant, ranging between 0.36 eV and 0.39 eV, during the second part of the process. PBE, PBE+U and RPBE yield approximately the same results concerning the generation of CHO*. The activation energy is estimated in ca. 1.20 eV, while the reaction energy is around 0.75 eV for PBE and RPBE calculations, and 0.67 eV for the PBE+U one, as reported in Figure 5. At odds with these results, RPBE description of the CHO* dissociation step deviates strongly from the PBE and PBE+U ones. As is clearly visible from Figure 6, starting from a certain point along the reaction coordinate prior to the first maximum, the RPBE energy diagram is shifted to higher energies when compared to the PBE and PBE+U cases. This upshift slightly varies with the reaction coordinate until the minimum corresponding to configuration B is reached, while it is almost constant, ranging between 0.36 eV and 0.39 eV, during the second part of the process. To understand the nature of such energy difference, we decomposed the energetics of the system on the basis of energy data coming from simpler systems. As shown in Figure 7,. we were able to decompose the system energetics in three main contributions: the adsorption energies of single CHO*, CH*, O* species, the energy of the free CHO, CH and O (and thus the CH-O bond strength), and the CH*-O*co-adsorption energy. The latter can be recovered on the basis of the dissociated CH*/O* system (the last image of NEB calculation) and the isolated ones. More details on these calculations can be found in the Supplementary Material. To understand the nature of such energy difference, we decomposed the energetics of the system on the basis of energy data coming from simpler systems. As shown in Figure 7, we were able to decompose the system energetics in three main contributions: the adsorption energies of single CHO*, CH*, O* species, the energy of the free CHO, CH and O (and thus the CH-O bond strength), and the CH*-O*co-adsorption energy. The latter can be recovered on the basis of the dissociated CH*/O* system (the last image of NEB calculation) and the isolated ones. More details on these calculations can be found in the Supplementary Material.  The analysis shows that, in line with the data reported in Table 2, the main difference between PBE and RPBE functionals can be ascribed to the way the free species and the adsorbed species energetics are modeled, while a minor contribution is given by the estimation of the co-adsorption energy. As shown in Table 2, RPBE has a tendency toward lower binding energies with respect to PBE schemes, and this tendency is particularly strong for CH* and O* where the PBE to RPBE adsorption energy differences are about twice the CHO* PBE to RPBE adsorption energy difference. These considerations can motivate the shift observed in the NEB profiles and indeed it can be noticed that the profiles in Figure 4 begin to diverge once the C-O bond is split, i.e., when the substrate-adsorbate interaction increases.

Materials and Methods
All simulations were carried out at the DFT level using a plane-waves basis set, as implemented in Quantum Espresso version 6.5. [75] The Perdew-Burke-Ernzerhof (PBE) [33] and the Revised Perdew-Burke-Ernzerhof (RPBE) [34] exchange-correlation (XC) functionals were employed to describe the electron-electron interactions. PBE+U calculations were also performed on systems containing C and O atoms, using as corrective parameter U = 0.65 eV both for C and O species. By using this U value, in line with the one adopted in similar previous calculations [55,56], the computed CO adsorption energies on Rh (100) surface match the experimental data. As the main process in this reaction concerned the generation and cleavage of covalent bonds, we did not include any dispersion corrections. The energy cut-off for the plane waves expansion was set to 1080 eV and 820 eV for PBE and RPBE calculations respectively. Core-electrons were modeled with the XC-compatible Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotentials. Bulk calculations for the Rh crystals were performed on a face centered cubic (fcc) cell, using an 18 × 18 × 18 k-point mesh, yielding an equilibrium lattice constant of 3.833 Å and 3.922 Å for the PBE and RPBE calculations respectively, in close agreement (especially in the case of PBE) with the experimental value [76] extrapolated at 0 • K of 3.797 Å. The PBE and RPBE calculated bulk cell parameters were used to build the 8-layers Rh (100) slabs which represent the catalyst surface. Building the supercells, a 20 Å vacuum layer was added along the direction orthogonal to the surface in order to avoid any interaction between consecutive replicas. Keeping the two inner layers fixed at their ideal bulk positions, all the remaining atomic positions were relaxed until the residual forces were lower than 10 −4 eV/Å, with a convergence threshold on the ground state total energy of 10 −9 eV. All relaxation runs were carried out using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton algorithm [77].
Geometry optimizations involving adsorbed species were performed on single molecules adsorbed on a 3 × 3 eight-layer-thick supercell. This supercell dimension was chosen in order to avoid spurious interactions among molecules belonging to nearby unit cells within the same surface. The calculations on these systems were performed using a 6 × 6 × 1 k-point mesh, the same Brillouin zone sampling of the converged bulk calculations ensuring adsorption energies accurate up to ±0.01 eV. Dipole corrections were applied along the direction orthogonal to the slab surface to prevent artifacts related to the presence of an artificial net electric field in the case of asymmetric slabs [78]. Calculations on the adsorbed species in gas phase were carried out by fully relaxing the respective geometries within fcc cells with a 40 Å lattice parameter. The study of the methanation rds was conducted using the Climbing Image Nudged Elastic Band technique (CI-NEB) [32]. Transition state geometries and energy barriers for the process were obtained by imposing a convergence threshold on the forces of 0.05 eV/A. The CI-NEB simulations was split into two separate calculations, namely the formation of the intermediate CHO* and its dissociation. The former calculation was conducted using 20 images, while the latter exploited a total of 40 images. Transition states were identified as the maxima in the NEB energy profile in analogy with previous study concerning CO 2 reduction over various metal surfaces [13,29]. Furthermore, to confirm the nature of the transition state associated with the rate-determining step, we performed vibrational frequency and normal mode calculations on the geometry corresponding to the maximum of the rds NEB profile, by computing the Hessian matrix on the solely atoms of the adsorbate, while the atoms of the slab were kept fixed [75]. The starting point of the CHO* generation simulation was the most stable co-adsorbed CO*/H* configuration, while the final image was the relaxed CHO* complex adsorbed on the Rh (100) surface. CHO* dissociation calculations started from this last image and ended in the configuration in which the complex was dissociated, and the CH* and O* relative positions minimized the total energy of the system.

Conclusions
In conclusion, we studied the rate-determining step of the rhodium catalyzed CO 2 reduction, which is the hydrogen-assisted CO* dissociation through formation and subsequent dissociation of the CHO* intermediate [4,5,13,15]. In the present work, this reaction step was studied by means of DFT and CI-NEB techniques on the Rh (100) surface. Within this framework, the PBE, RBPE XC functionals and PBE+U methods were employed. Overall, all the methods predicted the same energetics for the CHO* generation from CO* and H*, which was found to be a single step process. On the contrary, the CHO* dissociation was found to be happening in two separate steps: first the elongation and breaking of the C-O bond, and second the migration of the O* species to its equilibrium site, with the concomitant relaxation of surface tension. Although finding the same process and geometries, RPBE results differ from PBE based ones when taking a closer look at the energies which are involved in the CHO* → CH* + O* step: the products being less stable by~0.35 eV and the C-O bond break energy barrier strongly overestimated for RPBE. The origin of this difference resides in the way RPBE was designed in order to correct the adsorption energetics of LDA and PBE functionals [34,62]. The tendency of RPBE to underbind with respect to PBE is especially enhanced in the case of the product species. For example, our results show that RPBE strongly underbinds the O* species with respect to experimental data, largely overcorrecting PBE overbinding. The larger underbinding in the O* and CH* species compared to CHO* explains both the deviation of the RPBE profile with respect to the PBE based one in the moment the C-O bond breaks and the product species are formed, and the subsequent almost rigid shift of the profile. A similar picture to that found on Ru (0001) concerning the energetics of the two subprocesses was found, suggesting the CHO* cleavage as the rate-determining step of the whole Rh-catalyzed methanation reaction.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/catal11050538/s1, Video SM: CI-NEB calculated mechanism for the CHO* generation and dissociation; Text: details on the calculations of the interaction energy among co-adsorbed species, snapshots of the meaningful involved in the PBE+U calculated NEB trajectory.