Selectivity Trends and Role of Adsorbate-Adsorbate Interactions in CO Hydrogenation on Rhodium Catalysts

Predictive-quality computational modeling of heterogeneously catalyzed reactions has emerged as an important tool for the analysis and assessment of activity and activity trends. In contrast, more subtle selectivities and selectivity trends still pose a significant challenge to prevalent microkinetic modeling approaches that typically employ a mean-field approximation (MFA). Here, we focus on CO hydrogenation on Rh catalysts with the possible products methane, acetaldehyde, ethanol and water. This reaction has already been subject to a number of experimental and theoretical studies with conflicting views on the factors controlling activity and selectivity towards the more valuable higher oxygenates. Using accelerated first-principles kinetic Monte Carlo (KMC) simulations and explicitly and systematically accounting for adsorbate-adsorbate interactions through a cluster expansion approach, we model the reaction on the low-index Rh(111) and stepped Rh(211) surfaces. We find that the Rh(111) facet is selective towards methane, while the Rh(211) facet exhibits a similar selectivity towards methane and acetaldehyde. This is consistent with the experimental selectivity observed for larger, predominantly (111)-exposing Rh nanoparticles and resolves the discrepancy to earlier first-principles MFA microkinetic work that found the Rh(111) facet to be selective towards acetaldehyde. While the latter work tried to approximately account for lateral interactions through coverage-dependent rate expressions, our analysis demonstrates that this fails to sufficiently capture concomitant correlations among the adsorbed reaction intermediates that crucially determine the overall selectivity.


Introduction
The conversion of syngas (CO and H 2 ) into hydrocarbons and oxygenates is attractive as an alternative source of fuels and chemicals. However, selectivity towards the more useful higher oxygenates such as ethanol and acetaldehyde remains challenging, with methane being a common undesired product. 1,2 Many theoretical and experimental works have focused on understanding and tuning especially the selectivity of Rh catalysts, as Rh is generally recognized as one of the most promising elemental catalysts for the direct synthesis of higher oxygenates. It has now become clear that pure Rh catalysts are intrinsically selective towards primarily methane and acetaldehyde, 3,4 whereas ethanol synthesis requires promotors such as Fe or Mn. 3,[5][6][7][8][9] Recent experimental works have suggested that there is an inverse relationship between activity and selectivity for pure Rh catalysts, where an overall higher activity (CO conversion) correlates with a lower selectivity towards acetaldehyde. Explanations offered for this trend are, however, conflicting. Yang et al. have suggested, based on density-functional theory (DFT) and mean-field microkinetics, that it is the nature of the active sites exposed by the catalyst nanoparticles that is the deciding factor, with step sites being highly active and selective towards methane and terrace sites being less active and selective towards acetaldehyde. 3 This view was recently challenged by Schumann et al. 4 They synthesized Rh nanoparticles of different sizes and found that it is primarily small particles below 2 nm that exhibit high acetaldehyde selectivity and low activity, while the larger particles above 5 nm are the most active and selective towards methane. The surface fraction of edge/corner sites increases for smaller particles and was found to closely follow the selectivity trends. While this could indicate that step or corner sites are the active sites for acetaldehyde synthesis, Schumann et al. proposed instead that the smaller particles support a much higher local CO coverage at both terrace and step sites, which limits the activity due to poisoning, but increases the acetaldehyde se-lectivity by driving C + CO coupling reactions.
The conflicting views offered in the literature clearly call for more thorough theoretical investigations and improved microkinetic models that can account for the effects of high surface coverages and concomitantly increased lateral adsorbate-adsorbate interactions. In the meanfield study carried out by Yang et al. coveragedependent rate equations were parametrized and employed to mimic such interactions. However, it is well recognized in the literature that mean-field kinetics cannot properly account for effects of correlations and fluctuations, including fluctuations in the local coverage. [10][11][12][13] A more accurate approach is KMC simulations using a cluster expansion (CE) to treat lateral interactions. [14][15][16][17][18] Unfortunately, this can become very expensive for complex lateral interaction models and large disparities in the time scales of the different processes. In a recent work by Chen et al. on syngas conversion at the Rh(111) surface, some of these challenges were avoided by making the assumption that there is a clear time scale separation between fast diffusion processes and slow reactions such that the system can be solved by alternating between separate KMC and mean-field models. 19 This assumption may, however, not always be applicable as we will show in this work and as it has also been demonstrated in previous literature works. 10,20 In this work we revisit the question of the role played by step and terrace sites for activity and selectivity trends in syngas conversion over Rh catalysts. We consider the pristine Rh(211) and Rh(111) facets, which are representative of step and terrace sites found at larger nanoparticles where finite-size effects do no longer play a large role. Owing to recent methodological developments in our in-house KMC code kmos 21 concerning the efficient modeling of lateral interactions and the implementation of an acceleration algorithm 22,23 to tackle the time-scale disparity problem, we are able to carry out full-blown KMC simulations employing a CE model for lateral interactions. The results with and without account of lateral interactions are compared to corresponding mean-field kinetics.
The main finding of our work is that -in contrast to previously parametrized coverage-dependent mean-field models -KMC simulations that explicitly and systematically account for lateral interactions are able to correctly capture the experimental selectivity trends for large nanoparticles. We show and rationalize why in some cases the lateral interactions have a huge effect on the results (Rh(111) facet), whereas in other cases the effects are negligible (Rh(211) facet). The finding of a breakdown of mean-field kinetics is not restricted to models with lateral interactions. In fact, we show that also in the absence of lateral interactions, reaction-induced inhomogeneities and diffusion limitations can cause the meanfield-predicted activities to deviate substantially from the KMC results (Rh(211) facet). Finally, we show that in order to reach a quantitative agreement with both the selectivity and the activity trends observed in experiments, we need to correct for well-known errors in DFTpredicted adsorption energies such as the persistent overbinding of CO with standard semilocal DFT functionals. [24][25][26] Methods

Reaction model and parametrization
The reaction networks employed for the Rh(211) and Rh(111) facets are shown in Figure 1 and were inspired by the work of Yang et al. 3 As in our previous works on CO hydrogenation over stepped metal surfaces, 23,27 and extending over the work of Yang et al., we use a highly resolved active site representation for the Rh(211) facet consisting of a terrace site t, an upper step site s and a lower four-fold coordinated step site f . Only the two latter sites were considered by Yang et al. The full reaction networks we employed can be found in Section S1 of the Supplementary Information (SI). DFT data calculated at both low and high CO coverage were taken from Yang et al. and used to parametrize a CE model, which is used in connection with the KMC simulations. Some additional DFT calculations were carried out by us using the Quantum Espresso code 28 with the BEEF-vdW functional 29 and the exact same numerical settings as used by Yang et al., see Section S2 in the SI.
Adsorption energies are referenced to CO, H 2 O and CH 3 OH using the formation energy approach with the formation energy E x form of adsorbate x, the total energy of the adsorbed species on the slab E x slab+ads , that of the clean slab E slab , and the gas phase reference energy µ i of atom i as well as the occurrence n x i of atom i in the adsorbate x. Adsorbate-adsorbate interactions are incorporated using a CE which is terminated after the first term corresponding to pairwise nearest neighbor interactions with the occurrence n M x and formation energies E x form of the individual adsorbates x in structure M and the occurrence k M j and interaction energy ε j of the pairwise interaction j. The interaction energies (cf. Table S2 in the SI) are obtained by solving a set of linear equations for a number of different structures M . As the surface at the investigated conditions is mostly covered with CO, we consider only self-interactions between the CO molecules, interactions between CO and the remaining adsorbates (except H), as well as selected other interactions (e.g. between C and CH at the f site) following Yang et al. All other interactions are neglected due to their low probability of occurrence. Following common practice, 3,30 H is adsorbed at special "hydrogen reservoir" sites, which reflects the assumption that it can intercalate into sublattice sites because of its smaller size compared to the other adsorbates. In the KMC models, we implement this by treating H in a mean-field ansatz, following our previous work. 23 This means that it is not an actual species but only represented via an effective H coverage. We use different energetics (and thereby coverage) depending on whether the H atom is involved in reactions taking place a) CO CO at an s or a t site (for H we denote these reservoir sites as S and T sites both in the KMC and mean-field models). Our treatment thereby implicitly assumes that there are no spatial correlations in the distribution of H at the surface, which should be a good approximation when taking into account that H diffusion and adsorption/desorption are fast processes. Note, however, that the assumption that H does not block any surface sites, nor interacts with the other adsorbates, would most likely break down under high-pressure reaction conditions. For the rate constants of elementary steps corresponding to adsorption, desorption, reaction, and diffusion, we used standard expressions including zero-point energies and other enthalpy/entropy corrections within the harmonic approximation (adsorbates and transition states) or the ideal gas approximation (gas species) from the thermochemistry module of the Atomic Simulation Environment software package, 31,32 cf. our previous work. 23 Vibrational frequencies are taken from Yang et al. 3 In our previous work, 23 diffusion steps of CO, O, OH, CH, CH 2 , and CH 3 were found not to be rate-limiting and therefore we stick to approximate barriers calculated for Re(0001) from Hahn et al. 33 For same-site diffusion of C at f sites and CH at t sites we use values calculated for Rh(211) from our previous work. 23 Diffusion of all other species is neglected. To allow for comparison of our simulations to those of Yang et al., 3 we use the same reaction conditions of p CO = 13.33 bar, p H 2 = 6.66 bar, and with the energy shift ∆E FS (∆E IS ) of the final (initial) state due to lateral interactions, the zero-coverage barrier E a 0 , and the parameter α with values in the interval [0, 1] representing a reactant like (0) or product like (1) transition state. DFT-calculated energetics for determining the α parameters shown in Table S3 are taken from Yang et al. and Andersen et al. 3,23 Although BEP relations entail some energetic uncertainty, the magnitude is on the order of the error introduced by semi-local DFT. 14 The approach of combining CE with BEP relations is used across different KMC frameworks and effectively reduces the computational burden for the inclusion of lateral interactions. Previous studies employing this approach were able to quantitatively capture experimental observations without the need to explicitly calculate all possible lattice configurations, which is especially critical with an increasing number of sites and adsorbates in the reaction network under consideration. 14,36-38

Kinetic Monte Carlo
The KMC simulation technique allows for a numerical solution to the time evolution of the spatial distribution of the adsorbates on the coarse-grained sites of the catalytic surface, along with related properties such as the occurrence of individual reaction pathways and the catalytic activity and selectivity in terms of the turn-over-frequencies (TOFs) for the formation of the reaction products. Since individual elementary steps are executed step-by-step with probabilities reflected by their rate constants, KMC can suffer from a performance bottleneck in case of processes that occur on disparate timescales, e.g. fast diffusion and slow reactions. This timescale disparity challenge is here tackled through an acceleration algorithm developed by Dybeck et al. 22 In a recent work 23 we implemented this algorithm in the kmos code 21 and used it to study trends in CO methanation activity over stepped transition metals. While some challenging cases where observed, the algorithm was found to work well for the mainly CO-covered Rh(211) facet. Our previous work disregarded lateral interactions, however, and these can play an important role on the outcome of a reaction as shown in this work. Both repulsive and attractive interactions with neigh-boring adsorbates can alter the energetics of elementary steps and lead to changes in their individual rates. In this work we used the acceleration algorithm in connection with the recently developed on-the-fly backend in kmos. 39,40 This backend features increased performance and reduced memory requirements for models with many lateral interactions compared to the original backend, since the rates of processes affected by lateral interactions are calculated at runtime according to the above-mentioned CE model and BEP relations. The error bars for the TOFs are obtained from a Bayesian error analysis. Reactions in KMC follow the Poisson distribution where P (n|v) is the probability for observing n turnovers during a fixed simulation time t, given the expected value v. Here, we are rather interested in calculating the probability distribution for v, given a (possibly small or even zero) number of observed turnover events n. This posterior distribution can be obtained via Bayes theorem with the prior probability distribution P 0 (v), which we assume to be a constant C as there is no a priori information on the TOF. To normalize P 0 (v) we set a v max , i.e. C = 1/v max , to a value high enough that the following approximation holds This gives the posterior probability distribution which is also a Poisson distribution. Note that this also gives us the probability distribution for the TOF, since the TOF is simply v divided by the simulation time t. The most likely value of the TOF is determined from the maximum of the posterior distribution, and the upper and lower bounds of the error bars are obtained as the smallest credible interval that contains 99% of the total probability mass. Sensitivity analysis was carried out by calculating the degree of rate control (DRC) proposed by Campbell and coworkers 41 for each TS i (X RC,i ) with the rate r, the free energy G i of TS i, the universal gas constant R and the absolute temperature T . In KMC the derivative was approximated by the finite-difference expression The pair probabilities of second order processes are obtained by storing the timeintegrated counts of lattice configurations in which the process can be executed over the entire simulation, and the counts are then normalized by the total simulation time and the size of the system. All simulations were run for 5×10 7 steps to reach a steady-state and subsequently until 11 (Rh(111)) or 26 turnovers (Rh(211)) of acetaldehyde were observed or the total KMC simulation time reached one week. The simulations were carried out using a (10×10) lattice for Rh(211) and a (25×25) lattice for Rh(111) with periodic boundary conditions. Further details about the KMC simulation settings and the convergence tests carried out for both the lattice size as well as the parameters employed in the acceleration algorithm can be found in Section S3 of the SI.

Mean-field microkinetic modeling
For comparison to the KMC simulations, we also carried out microkinetic simulations using the CatMAP software package. 42 CatMAP employs the MFA, meaning that the spatial distribution of the adsorbates is further coarsegrained into a mean coverage of each site type, thereby also neglecting coverage fluctuations around the mean. This results in a set of coupled rate equations, which are solved at steady state. All of our MFA simulations were performed without lateral interactions, as the parametrization of coverage-dependent rate equations is not a topic of this work. We emphasize that the MFA simulations used exactly the same reaction network and adsorption energetics as our KMC simulations, so that differences between both simulation approaches are entirely due to the MFA employed in the prior.

Results and discussion
We begin by presenting MFA and KMC results for the Rh(211) facet in Figure 2a. The KMC results are reported at three different temperatures without (left offset) and with (right offset) lateral interactions. Without lateral interactions, a comparison of our MFA results (solid lines) and KMC results (left offset) reveals that the activities are slightly overestimated in the MFA model by about a factor of 2.7-6.0 and that both models show similar selectivity trends towards methane rather than acetaldehyde. Ethanol turnovers are not observed in the KMC simulations -hence the large error bars -which is consistent with the very low TOFs obtained in the MFA model. We will come back to the reason for the differences in actual TOF prediction between MFA and KMC below.
Turning to the effect of lateral interactions, we can observe from the KMC results that they do not influence the TOFs much. Our MFA results are only presented without lateral interactions as the parametrization of these are not available for our modified active site representation. However, in the original Rh(211) MFA model from Yang et al. (see Figure S1 in the SI), lateral interactions do influence the TOFs, especially at lower temperatures. At 523 K the difference ranges from a factor of about 300 for ethanol to about 4000 for acetaldehyde.
Before diving into a deeper analysis of the differences between the different models for the Rh(211) facet, we present also the TOFs for the Rh(111) facet, see Figure 2b. This facet contains only a terrace site, and thus the active site representation and reaction network we employ (cf. Section S1.3 in the SI) is in this case completely identical to that of Yang et al. Comparing the MFA results with and without interactions (dotted and solid lines) we can see that for this facet the influence of the lateral interactions is huge -at all temperatures the catalyst is essentially inactive if lateral interactions are not taken into account. Furthermore, it is surprising to note that the lateral interactions in the MFA model, which were parametrized using a second-order expansion in the coverage, give rise to a change in the selectivity trends. That is, without interactions the catalyst is always selective towards methane, but with interactions the catalyst is selective towards acetaldehyde at the experimentally relevant temperature below about 600 K. We recall that the MFA results for Rh(211) and Rh(111) with interactions were used by Yang et al. to propose that the activity-selectivity trends obtained with experimentally synthesized Rh nanoparticles are caused by varying amounts of step and terrace sites at these particles, where step (terrace) sites would then be the active sites for methane (acetaldehyde) formation, respectively. 3 However, as discussed in the introduction, this explanation was recently challenged by Schumann et al. 4 based on new and more detailed experiments that showed that large particles above 5 nm, which are expected to predominantly expose the Rh (111)  we find that the Rh(111) facet is indeed selective towards methane, see Figure 2b. Without interactions we do not observe any meaningful TOFs for any of the products (not shown), which is consistent with the very low TOFs obtained in the MFA model without interactions. Thus, our KMC results confirm that lateral interactions are of paramount importance for the Rh(111) facet, but the CE used to parametrize these interactions in KMC gives qualitatively different results to the MFA-parametrized lateral interaction model. Importantly, only the explicit site-resolving KMC results are able to reproduce the experimentally observed selectivity trends for large particles. The case of small particles is outside the scope of the present work as the low-index Rh(211) and Rh(111) surfaces mainly represent step and terrace sites at larger particles free from finite-size effects.
We believe that the shortcomings of the MFA model for both Rh(211) and Rh(111) are caused by the well-known problems of these models with accounting for effects of correlations and fluctuations, including fluctuations in the local coverage. [10][11][12][13] In the following we exemplarily analyse this for the Rh(211) surface without interactions. This analysis starts by performing a sensitivity analysis to identify the rate-limiting steps (RLSs) of the two models. As shown in Figure 3, the RLS for both the methane and the acetaldehyde pathways is mainly water formation at the t site in the KMC model, whereas for the MFA model it is mainly methane formation and CH-CO coupling at the s site. The sensitivity analysis for the KMC model with interactions (see Figure S3 in the SI) is very similar to that of the KMC model without interactions.
Since a breakdown of the MFA is typically associated with second-order reaction steps where correlations in the spatial distribution of the two reacting species at the surface occur, 23 we next evaluate the pair probabilities for selected elementary steps from our KMC simulations. In Figure 4 we plot the ratio between the KMCsimulated pair probability to find the reacting species A and B at neighboring sites and the MFA-assumed probability equal to c[A] [B], where c is a geometric factor (the site connec-

tivity) and [A] ([B]) is the surface coverage of species A (B)
. It is seen that the MFA indeed breaks down for several of these steps, i.e. the ratio is significantly different from one. In particular, we note that the MFA breaks down for the reaction between C at the f site and OH at the t site to form CO at the t site and H. This is the reverse of the CO dissociation step, which is executed just before the RLS in the KMC model where OH at the t site is further hydrogenated to form water. The probability of finding the two reactant species, C * f and OH * t , next to each other is about 57 times larger in KMC than in MFA. This is caused by the high barrier for C to diffuse along the step and the high CO coverage at the t site, which effectively hinders OH diffusion at the t sites, cf. our previous work on the comparison of MFA and KMC models for CO hydrogenation at stepped metals. 23 Figure 4: Ratio between KMC and MFA probabilities to find the pair of reacting species in selected second-order reaction steps at neighboring sites in the Rh(211) model without lateral interactions. The analysis is performed at 650 K.
As a consequence of this persistent correlation, it becomes more probable for C * f and OH * t to react back to form CO * t and H in the KMC model compared to the MFA model. In a free energy diagram, cf. orange and black lines in Figure 5a, this change in pair probability can be represented as an increase in the free energy barrier for the reverse reaction step of C * f + OH * t in the MFA model (orange line) caused by a decrease in the free energy of the C * f + OH * t state (i.e. a higher configurational entropy in the MFA-assumed well-mixed state). Since the free energy barrier for further reaction out of the C * f + OH * t state (i.e. water formation at the t site) is unchanged, the effective barrier in the MFA free energy landscape (difference between H-OH * t transition state and CO * t state) is smaller than the effective barrier in the KMC free energy landscape. This explains the larger TOF obtained in the MFA model and the change of the RLS from water formation at the t site to methane formation and CH-CO coupling at the s site (the latter steps are not shown in Figure 5a).
From the free energy diagram in Figure 5a Figure S5 where each CO has three neighboring CO molecules, we instead determined the energy levels and barriers specifically for this surface configuration. The analyses are performed at 650 K.
we can also explain why the KMC TOF and RLS for both methane and acetaldehyde formation at the Rh(211) facet does not change much upon inclusion of lateral interactions. As seen by comparing the black curve (without interactions) and the blue curve (with interactions), the free energies of the states that determine the effective barrier (H-OH * t transition state and CO * t state) move up in energy by a similar amount as a consequence of interactions, causing the effective barrier and TOF to remain at a similar value. The detailed lateral interaction parameters at play in the two states are illustrated in Figure S4 in the SI. While unimportant for the Rh(211) model, we already mentioned above that lateral interactions do have a huge effect on the Rh(111) model, and we will next analyze why this is the case. The free energy diagram with and without lateral interactions for the Rh(111) facet is shown in Figure 5b. The biggest difference compared to the Rh(211) facet is that here lateral interactions actually push the free energy of the CO * t state up close to the free energy of gas-phase CO. Thereby, the average CO coverage decreases from around 100% without interactions to about 62% with interactions at the analyzed temperature of 650 K, cf. Section S3.3 in the SI. In Figure 5b the CO * t state is split into two levels since CO desorption primarily happens from a local high-coverage state where on average 5.1 out of 6 neighboring sites are occupied by CO, whereas hydrogenation of CO primarily happens from a lower-coverage state where only 3 neighboring sites are occupied by CO. The average coverage pattern with about 2 /3 CO coverage resembles a honeycomb lattice of CO molecules as illustrated in Figure S5 in the SI. The RLS for both methane and acetaldehyde formation in the model with interactions is primarily dissociation of CHOH to form CH and OH (and for acetaldehyde also the CH 2 CO hydrogenation step), cf. Figure S6 in the SI.
For Rh(111), the profound impact of lateral interactions on the simulation results cannot be explained by a change to the effective free energy barrier (i.e. the difference between the rate-limiting CH-OH * t TS and the CO * t state in Figure 5b), since this barrier decreases with only around 132 meV upon inclusion of lateral interactions. Rather, the important point is that for CHOH to dissociate, a neighboring t vacancy is required. Without lateral interactions, the surface is completely poisoned with an average CO coverage close to 100%. Thus, the probability for CHOH dissociation to oc-cur is negligible, which explains why no product formation could be observed in the KMC simulations. With interactions, however, a formed CHOH will on average have 3 CO and 3 vacancy neighbors (see Figure S5 in the SI), which makes the dissociation reaction possible. For the Rh(211) facet a similar CO poisoning is not observed, since here CO dissociation relies on a vacant f site, and the vacancy coverage of the f sites remains high even without lateral interactions, cf. Table S4.
Up until now, we have shown that the MFA can break down both in the absence of lateral interactions (analyzed for the Rh(211) facet where the TOFs are overestimated in the MFA model compared to KMC) and in the presence of lateral interactions (analyzed for the Rh(111) facet where wrong selectivity trends are obtained in the MFA model compared to experiments). However, even if the KMC model does recover the experimental selectivity trends, the absolute TOF predicted by both the KMC model and the MFA model is not in good agreement with the experiments.
In Figure 2a and Figure 2b we mark with a blue diamond the total TOF (primarily methane) measured for the largest Rh nanoparticles investigated by Schumann et al. 4 It is larger by about five orders of magnitude compared to the Rh(211) KMC simulation with lateral interactions (cf. Figure 2a) and by about 5-6 orders of magnitude for the Rh(111) KMC and MFA simulations with lateral interactions (cf. Figure 2b). The consistent, gross underestimation of the TOF by both KMC and MFA and for both investigated facets leads us to suspect that this discrepancy has its origin in the DFT calculations used to parametrize the kinetic models. In particular, it is well-known that generalized gradient approximation (GGA) functionals tend to overbind CO on transition metals, 26 e.g. the CO adsorption energy obtained with the here employed BEEF-vdW functional for Rh(111) (−1.7 eV) is by 0.25 eV stronger than the experimental value (−1.45 eV 25,43 ). Since the state with CO adsorbed at a terrace site directly influences the effective barrier of the reaction, cf. Figure 5, this error directly influences the obtained TOF. In order to assess the possible implications of this error, we show in Figure 6 the TOFs and selectivities obtained for the Rh(211) and Rh(111) facets for both the hitherto discussed base model and a model where the CO adsorption energies are increased by 0.25 eV (E ads CO ↑). This already improves the agreement with experiments significantly, i.e. the TOF predicted for the Rh(211) facet (Rh(111) facet) is now only about two (three) orders of magnitude lower than the experiment. Of course, this analysis neglects the fact that there could also be errors associated with other rate-controlling parameters in our kinetic model, e.g. the transition state for water formation at the terrace site or lateral interaction parameters affecting either the latter transition state or the CO adsorption energy. Since we do not know in which direction such other errors might point, we show as an example in Figure 6a and Figure 6b also results where the interaction parameter for selfinteractions between the CO molecules at t and s sites are increased (ε CO ↑) or decreased (ε CO ↓) by 50 meV from their base values. An error of this magnitude seems quite reasonable since we use an approximate CE that is terminated after pairwise nearest neighbor interactions. Taking these results as plausible "error bars" on the theoretical results, we can now reconcile theory with experiment to the extent that the experimental TOF lies within (admittedly rough) error bars on the theoretical values for the Rh(211) facet. In reality, the experimentally used nanoparticles of course contain both step and terrace sites, and the actual TOF measured might be the result of an interplay between reaction steps taking place at both site types. Such bifunctional effects have been previously discussed in the literature, 44,45 and have indeed recently been demonstrated theoretically for e.g. CO oxidation at Pt nanoparticles 46 and hydrogen evolution reaction at jagged Pt nanowires. 47 The finding here that the TOF is largest, and in best agreement with experiment, on the Rh(211) facet is related to the facile hydrogen-assisted CO dissociation at the step sites exposed by this facet. However, the subsequent rate-controlling water formation step takes place at the terrace site, thus, both types of sites are required for a high activity. An analysis of how this would play out at realistic nanoparticle geometries will be an intriguing topic for future work.
Comparing the experimental selectivities shown as dashed horizontal lines in Figures 6c  and 6d with the theoretical selectivities reveals that for both facets variations in the adsorption energies and interaction energies (as well as combinations) bring the selectivities closer to the experimental measurements. The observed trend of an increase in the acetaldehyde selectivity with a destabilization of adsorbed CO for both Rh(211) and Rh(111) is caused by a decrease of the effective barriers in the acetaldehyde pathways involving the RLSs of CH-CO coupling (for Rh(211) as shown in Figure 3) and CH 2 CO hydrogenation (for Rh(111) as shown in Figure S6). Nevertheless, for all parameter modifications, our results robustly indicate that the Rh(111) facet is selective to methane, whereas the Rh(211) facet exhibits a similar selectivity to methane and acetaldehyde. Given that realistic 5 nm nanoparticles will dominantly expose the Rh(111) facet, as well as step and corner sites to a smaller extent, the sum of the theoretical selectivity trends should overall result in a methane selectivity in quite good agreement with the experimental measurements.

Conclusions
We investigated the CO hydrogenation reaction on Rh(111) and Rh(211) using accelerated firstprinciples KMC simulations with and without lateral interactions parametrized from a cluster expansion model. The results are compared to MFA simulations from Yang et al. 3 with and without coverage dependence in the rate expressions. The coverage-dependent MFA model predicts acetaldehyde selectivity for the Rh(111) facet below ∼580 K, which is however in contrast to the KMC simulations with lateral interactions that predict methane selectivity at all temperatures. Importantly, only the KMC results are in agreement with recent detailed experimental investigations on selectivity trends of Rh nanoparticles from Schumann et al. 4 The inclusion of lateral interactions are found to have a huge effect on the Rh(111) simulations. We explain this from the fact that the RLS in this model (CHOH dissociation) requires a neighboring vacant terrace site, which is a probable lattice configuration only when taking into account the reduction of the surface CO coverage caused by repulsive CO-CO interactions. In contrast, for the Rh(211) model we find that the TOF and RLS (water formation at the terrace site) are only weakly affected by lateral interactions. This is explained from similar surface coverages with and without lateral interactions and from the fact that interactions shift the two states controlling the effective barrier in the dominant reaction pathway by similar amounts, leading to net similar effective barriers with and without interactions. For Rh(211) we furthermore show that the MFA breaks down also in the absence of interactions due to diffusion limitations and reactioninduced lattice inhomogeneities. This leads to higher TOFs and changes in the RLSs in the MFA model compared to the KMC model.
Finally, we compare the overall activities (methane and acetaldehyde TOFs) to the experimental measurements from Schumann et al. We find that we need to correct for the wellknown CO overbinding of GGA functionals (about 0.25 eV for the here applied BEEF-vdW functional) in order to approach agreement with the experiments. This correction also improves the quantitative agreement with the carbon selectivites obtained by Schumann et al. In particular, it leads to a larger increase in the acetaldehyde TOFs compared to the methane TOFs for both Rh(211) and Rh(111) since the effective barrier of the acetaldehyde pathway involves CO both in the CO activation step leading to CH x species at the surface (i.e. the step shared with the methane pathway) and in the subsequent step where the formed CH species reacts with another CO.
Overall, the insights obtained in this work could be relevant for further tailoring heterogeneous catalysts to improve their selectivity towards the desired higher oxygenates acetaldehyde and ethanol. The methodological ad- . Results are shown for CH 4 (blue), CH 3 CHO (orange), and CH 3 CH 2 OH (black). "base" represents the unmodified results also shown in Figures 2a and 2b. The dashed horizontal lines are the experimentally measured total TOF (primarily methane, (a) and (b)) and selectivities ((c) and (d)) of the largest nanoparticles with diameters above 5 nm from Schumann et al. 4 Both the theoretical and the experimental results are obtained at 523 K.
vances demonstrated here, i.e. the combination of acceleration algorithms in KMC with efficient modeling of lateral interactions, open up possibilities for treating also complex reaction networks at a level of detail beyond the hitherto applied approximate MFA models. This is important in order to achieve reliable mechanistic insights as a solid basis for the rational design of selective catalysts.

Supporting Information Available
Supporting Information. Additional details on the compared microkinetic models, DFT, and KMC. This material is available free of charge via the internet at http://pubs.acs.org. Input files for the kmos simulations for the different models together with additionally calculated DFT geometries are available at https://github.com/mdeimel/CO hydrogenation.git.