{'fraction_non_alphanumeric': 0.02667589936360339, 'fraction_numerical': 0.01167276684142843, 'mean_word_length': 4.595995940584925, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamic Microkinetic Modeling for Heterogeneously Catalyzed Hydrogenation Reactions: a Coverage-Oriented View

In most studies, the microkinetics for multistep reactions are numerically solved due to their complexity; the obtained numerical results are only valid under given reaction conditions at a static point. In this work, the microkinetics of heterogeneously catalyzed hydrogenation reactions are analytically solved as a function of three coupled physical parameters, which are energy, reaction rate, and coverage. The results correlate the surface reactions and the gaseous-phased reactant/product by energy and thus provide a dynamic view over the whole reaction process rather than at a static point. The analytical expressions are given for a simple hydrogenation reaction and three more complicated hypothetical hydrogenation reactions with side products, side reaction paths, or even multiple active sites. Compared with the numerical solution, the analytical solution is valid under all reaction conditions in practice and can provide more guidance to optimize the overall outcome or catalyst development.


INTRODUCTION
Synthesis of a number of valuable chemicals over catalysts has been industrialized and plays an important role in fields such as climate change and new energy. 1−3 Understanding the detailed mechanics of these catalyzed chemical reactions is of great importance.
Most important catalyzed chemical reactions involve multiple elementary steps; therefore, they can be modeled on two different levels. At the level regarding a single elementary step, the modeling is based on density function theory (DFT) and mainly addresses the atomic configuration change during the specified reaction step and the corresponding energy landscape. 4,5 On the other hand, the modeling can also be done on the microkinetic level, where the complicated correlations among all reactants, intermediates, and elementary reaction steps are the central problem.
There are some key concepts regarding the microkinetics, based on which some modeling methods are proposed. The most important and fundamental concept is energy due to its tight correlation to the reaction rate. According to transition state theory, the reaction rate of an individual step is dependent on its activation energy. There is also an apparent activation energy for a multistep reaction, which is determined by the activation energies of all elementary steps. The energetic span model (ESM) is proposed to conceptualize microkinetics from the energy aspect of view. 6−8 The model interprets the reaction rate [or referred as turnover frequency (TOF)] dependence on all the energetic parameters characterizing elementary steps; the "energetic span" is approximated by the energy difference between a pair of specific transition state and intermediate, which are referred as the TOF-determining transition state (TDTS) and TOF-determining intermediate (TDI), respectively. However, the ESM simplifies all steps with the reactants being converted to the products as a whole (e.g., A ↔ B), while most realistic reaction steps involve multiple reactants and/or products (e.g., A + B ↔ C + D). Besides, multiple reaction paths may exist and interfere with each other instead of a single reaction path (e.g., A ↔ B or C). In all these cases, the "energetic span" is very difficult to define, 9 as the microkinetics cannot be interpreted ignoring the complicated coupling among intermediates and paths.
A widely applied practical tool to analyze complex reactions is the "degree of rate control" (DRC). In this approach, the complete system is then treated as a black box with a number of tunable knobs, each knob representing a specific energetic parameter. The outcome of the black box is the TOF of the desired product, which can be numerically solved under given conditions, and varies whenever the knobs are being adjusted. A DRC parameter can be defined for each energetic parameter by how sensitive the TOF is while tuning the corresponding knob (similar to a partial derivative). 10,11 Under given reaction conditions, usually, only a single transition state and intermediate have a significant DRC (∼1), indicating that the TOF is very sensitive to their energetic parameters. The step corresponding to the transition state is referred to as the "ratedetermining step" (RDS), and the intermediate is referred to as the "rate-determining intermediate" (RDI) in the DRC approach, which are similar concepts to the TDTS and TDI in the ESM. The RDS and RDI can also be identified by a maximum rate-analyzing scheme. 12 As a numerical method, the DRC can be applied to complicated reactions and still reveals the relative significance of certain energetic parameters. However, it is still important to look into the black box and theoretically interpret the detailed mechanisms and their correlations to the DRC parameters. 13−15 The chemical reaction rate is dependent not only on the activation energy but also on the quantities of reactants/ products. In a catalyzed reaction, most intermediates require the catalyst molecules (site) to be formed. Due to the fact that the catalyst applied in industrial production is limited in amount compared to reactants, the total amount of intermediates is limited as well. As a result, the qualities of intermediates are determined by the catalyst partition among these intermediate states and affect the reaction rates of the corresponding reaction steps. For such reasons, catalyst partition should also be considered as a critical factor for the catalysis performance: in some cases, some intermediates may take up most catalyst, which is referred to as "poisoning" and can hinder the overall reaction. A reasonable partition among intermediate states is thus preferred. 16 Catalyst partition is also correlated with the selectivity among main/side reactions or paths; if intermediates of a side reaction take up most catalyst, it would hinder the main reaction and lead to low selectivity. 17 Besides, the energies of intermediates are also dependent on their occupancy of the catalyst through various mechanisms, such as configurational entropy 18 and adsorbate−adsorbate interactions. 19−23 Due to its correlation to many important variables, taking the catalyst partition into microkinetic modeling is critical and has various unique advantages, such as easily correlating with DFT results, capable of dealing with complicated reactions, and providing clear chemical physics insights. However, such an approach was only considered auxiliary. The primary reason is that in realistic reaction systems, most catalyst molecules/sites end up bound to a single intermediate, which is referred to as the TDI or RDI or the "most abundant reaction intermediate" (MARI). 24,25 Based on the idea "one TDI and one TDTS determine the overall catalysis process", further study on the catalyst partition after identifying the MARI was not considered to have a significant effect on the overall outcome.
In this work, a microkinetic model is built based on three coupled key concepts, which are energy, reaction rate, and coverage. Here, coverage refers to the catalyst partition among intermediates for heterogeneously catalyzed reactions. The reactions under study are heterogeneously catalyzed hydrogenation reactions, as the category includes a number of economical industrialized reactions. 1−3 Analytical solutions are given in a general expression, which can interpret the dynamic microkinetics through the complete reaction rather than only at a static point. Based on the analytical solution, the coupling among reaction rates, energies, and coverage in the steady state is revealed, and their expressions are given. In addition to the traditional energy-oriented view, the microkinetic results are also presented in a coverage-oriented way. Some important practical parameters such as the TOF, convert rate, and selectivity under varying different reaction conditions are also discussed in this framework. This paper is arranged in the following way: Section 1 briefly goes through the basic assumptions. Section 2 shows the microkinetic model for a simple hydrogenation reaction and how to interpret the microkinetics in the coverage-oriented way. Besides, the TOF and convert rate of the reaction are discussed. From Section 3 to Section 5, the approach is applied on three more complex reaction models, each extended from the model in Section 2 in a distinct way, to show its generality. Some important parameters and phenomena are analyzed and discussed in these sections as well. For each model, theoretical analysis and simulation are done for a given hypothetical reaction. The simulation result is compared with the theoretical result to provide necessary supports. Results and Discussions are presented separately for each model.

THEORETICAL FUNDAMENTALS
There are some general assumptions applied in the models throughout the work: First, it is assumed that all adsorbates involved occupancy of one active site. Second, it is assumed that the desorption energy barrier of a molecule can be approximated by its binding energy. These are common assumptions applied in microkinetic studies. Besides, the TOF is estimated under the steady state, where the coverage of all adsorbates does not change with time.
Also, it is assumed that all surface reaction steps follow the Langmuir−Hinshelwood mechanism, where the energy parameters are independent of the adsorbate−adsorbate interactions. In the case when such factors should be taken into consideration, all the adsorbate free energies applied must be corrected (by Bragg−William approximation, e.g.) to be consistent with the final coverage result. The logic here seems circular: in order to obtain the coverage result, free energies must be determined, which are again dependent on the coverage result. In practice, as long as the coupling relation can be estimated, this issue can be solved by iterating the simulation a few times to correct the results of energy and coverage until they are consistent. This mathematical approach is also adopted in DFT simulations to solve the coupling problem between the electronic density and potential. 5 The differences between the two cases are mainly in the DFT simulation part, so in the rest of the paper, the Langmuir−Hinshelwood mechanism is applied for the sake of simplification.
The analytical expressions of reaction rates as a function of coverage and free energies can be found in the Supporting Information.
2.1. Pure Hydrogenation Model. In this section, we discuss a pure hydrogenation reaction over a heterogeneous catalyst surface and show how the energy, coverage, and reaction steps are coupled within the microkinetic model. Here, "pure hydrogenation reaction" refers to a reaction with only hydrogenation elementary steps, except for the adsorption and desorption step. Such a reaction has the following form For example, hydrogenation of acetone, acetylene, and benzene can be classified into this reaction category. 19,22,23,26−28 The elementary steps of the reaction can be written in the form shown in Figure 1. In the figure, "*" indicates adsorbates or an empty site, and the activation energies refer to the forward reaction.
During such a reaction, each site on the surface is either empty or covered with hydrogen or intermediates R n (n = 0, 1, ..., N). Let the coverage of the abovementioned be written as θ # , θ H , and θ n . Through theoretical analysis (details regarding the deduction can be found in the Supporting Information), the coverage at the steady state can be expressed by four equations characterized by eight key energetic parameters, as given below Ä In the abovementioned equation, h and k are Plank's constant and Boltzmann's constant, respectively. T and r refer to temperature and TOF, respectively. E R i and E P i are the free energy of reactant R/product P in the transition state of the adsorption/desorption step, respectively; E R g and E P g are the free energy of the gaseous reactant and product, respectively; E l and E h refer to the lowest and highest total free energy level among all intermediate states, respectively, where the hydrogen species are counted in the gaseous form; E A refers to the overall free energy barrier for hydrogenation steps. E H refers to the free energy difference between hydrogen gas and the adsorbate per atom. γ is a coefficient regarding net reverse reaction with the expression given in (eq 5). The precise definition of the key energetic parameters is given in the Supporting Information; they are also marked in Figure 2a.
In previous studies, microkinetics of a multistep reaction are usually modeled and presented through an energy-oriented way. Figure 2a gives the typical energetic curve of the hydrogenation reaction, and the key energy levels in eqs 1−4 are marked. Alternatively, based on eqs 1−4, the microkinetics can be presented in a coverage-oriented way, as shown in Figure 2b: the catalyst surface is divided into six parts; each corresponds to 1−2 terms in eqs 1−4. It must be noted that the "part" here is not classified by intermediates, as in part 3, 4, and 5, there are multiple different intermediates (R 0 through R N ). The "part" here is rather classified by the reaction mechanism, as discussed below: Part 1 corresponds to the reactant adsorption. This part indicates the empty sites required to provide a reactant adsorption rate of r. This part only has one term, which is dependent on the adsorption energy barrier, where the entropy contribution is of great importance.
Part 2 corresponds to the reactant desorption from the reverse reaction. The adsorbates may be converted back to reactants through reverse reactions and desorb from the surface, thus reducing the net adsorption rate. Some extra empty sites are required in order to compensate the contribution from reverse reactions and reactant desorption. This part is made of two terms, which are corresponding to the intermediates represented by part 3 and 5. The intermediate represented by part 4 does not result in a separate term but causes all empty site expressions to have the coefficient γ.
Part 3 corresponds to the hydrogenation reaction. This part indicates the coverage of all intermediates required for all the hydrogenation steps to have a forward rate of r. It has only one term, which is dependent on the largest energy barrier E A among the series of hydrogenation reactions. It is also dependent on the hydrogen coverage for obvious reasons.
Part 4 corresponds to the saturation caused by the reverse reaction or, in other words, product readsorption. At the beginning of the reaction in the absence of the product, this part is 0. As the reaction carries on, gaseous products become  significant, and they can readsorb to the surface and undergo a reverse reaction, thus reducing the overall desorption rate. Some extra intermediates are required to compensate this part of readsorption and the surface reverse reaction. It is mostly determined by the free energy of the most stable intermediate and gaseous product, while also being dependent on the empty site coverage, as it also determines the readsorption rate of the product. It plays a very important role, which will be discussed later.
Part 5 corresponds to the product desorption. Similar to the reactant adsorption part, it only indicates the coverage of all intermediates required in order to have a net forward desorption rate of r.
Part 6 is the hydrogen coverage. This part has only one term; it is dependent on the empty site coverage and the free energy difference between hydrogen gas and the hydrogen adsorbate.
In Figure 2b, the six parts look comparable to each other in size; in practice, we will usually find one or few parts dominant due to the nature of their exponential form. The dominant part(s) has a coverage of ∼1 and is rate-determining. We will show the general procedure of applying this microkinetic model in Section 2.3.
2.2. Comparison with the Energetic Span Model and DRC Approach: Theory. Due to the exponential expression, usually in Figure 2b, one of the six parts is dominant over others, and its coverage is ∼1. By approximating that dominating part as 1, the expression of θ H and θ # is substituted if necessary, and its expression can be written in the following form In the expression given above, C refers to possible coefficient γ when the dominating part is either part 1, 2, 4, or 6. Note that although part 2 has two terms, one usually dominates over the other. The abovementioned expression can be rewritten as Expression (eq 7) is essentially the result given by the ESM, with E representing the energetic span or the apparent activation energy in experimental studies. The ESM has a more detailed expression, 29,30 and similar results to that expression can be obtained in this coverage-oriented framework by substituting (eq 1) through (eq 3) to (eq 4), and r can be solved. However, when there are multiple comparable parts, it is difficult to get an expression similar to the form of (eq 7); this is the disadvantage of the ESM.
To compare with the DRC approach, the unity condition (eq 4) is rewritten according to Figure 2b hr where E X is the index of any exponential term in Figure 2b; note that in the case that there is the term θ H or θ # , it has to be written in its explicit form. C X is the possible coefficient γ for the term. The abovementioned expression can be rewritten in the following form The DRC of a reaction, according to its definition, can be written as Ä where k i is the reaction rate of the specified reaction as in "krepresentation". Now (eq 9) is substituted to (eq 10), and we shall have Ä where θ X is the corresponding coverage of part X (X = 1, ..., 6). Equation 11 correlates the coverage-oriented model and the DRC approach. If one part is dominating (∼1), according to (eq 11), reaction steps directly related to the index expression of the dominating part (partial derivative is 1) in Figure 2b are the RDS, as marked in red in Figure 2a. Similar results can be obtained for any intermediate, which is also suggested in a previous study Ä Once again, in the case of a single dominating part, from (eq 11) to (eq 12), we see that all reaction steps and intermediates end up with a DRC of either 1 or 0. However, when there are multiple comparable parts, the DRC may end up with some other values.
Beside the transition state and intermediates, the abovementioned approach can also be applied to determine the order of reaction by taking the derivative to gaseous reactant free energy.
2.3. Comparison with the Energetic Span Model and DRC Approach: Numerical. Consider that the following hypothetical hydrogenation reaction takes place at kT = 0.04 eV: H 2 + R ↔ P. The detailed thermodynamic parameters characterizing the reaction are given in the Supporting Information, and their energetic-oriented presentation is given in Figure 3a. On top of that, microkinetics can be presented in the coverage-oriented way, as shown in Figure 3b.
The coverage of each part can be solved by simple algebra by substituting the energetic parameters as shown in Figure 3b. Part 3 is dominant in this case, and by approximating it as 1, the solutions are θ H ∼ e −2.5 = 0.082 and r ∼ kT/h e −27.5 ∼ 10 s −1 . Besides, the energetic span of part 3 (∼27.5 kT = 1.1 eV) can be considered as the apparent activation energy of the complete reaction.
The numerically solved solutions are θ H = 0.066 and r = 7 s −1 ; the coverage of some other intermediates can be found in Figure  4c in a later section. These results are within the same order of magnitude with those from the analytical approach, and the ∼30% error is due to the rough approximation of taking the coverage of part 3 as ∼1, while taking part 2 and 6 as ∼0. If we normalize their coverage according to Figure 3b with more accuracy, their coverage should be 0.86, 0.07, and 0.07, respectively, and a more accurate TOF result can be obtained.
For the transition state from R 1 to R 2 , its free energy is directly correlated to both E A and E h , as shown in Figure 3a, so it is directly correlated with all the significant parts (part 2, 3, and 6 in this example). The derivatives of the exponential indexes of these parts to the free energy of the transition state are all ∼1. Note that while taking the derivative of part 3, the hydrogen coverage must be written in a detailed form. According to eq 11, the DRC of the transition state is the net coverage of the three parts, which is ∼1.
For the intermediate R 0 , its free energy is directly correlated to E A and E l . The derivatives of the exponential indexes of parts 2, 3, and 6 to the free energy of R 0 * are approximately 1, −2, and 1, respectively. To interpret the −2 here, stabilizing R 0 not only increases the overall hydrogenation energy barrier but also decreases the hydrogen coverage. The coverage of part 2, 3, and 6 is 0.07, 0.86, and 0.07, respectively, according to the numerical result. According to (eq 12), its DRC is thus 0.07 − 2 × 0.86 + 0.07 = −1.58. By numerical simulation, the TOF as a function of energy variance on the two free energies is plotted in Figure 3c, and the DRC is 0.9996 and −1.5937, respectively, a good match with the analytical result.
2.4. Results and Discussion: Maximum TOF. Maximizing the TOF of the desired reaction is the key function of a catalyst.
With the analytical solution, we are able to predict the maximum TOF for a given reaction.
There are a few additional assumptions in the following discussion: First, it is assumed that the free energy of the reactants can be adjusted freely and independently by tuning the pressure. Also, the maximum TOF is only achieved before saturation becomes significant; we shall neglect the saturation part in this section; this corresponds to the beginning of a reaction with the absence of the product. Finally, the maximum TOF is evaluated at a given temperature, as increasing the temperature can usually raise the TOF.
Based on the coverage-oriented view given in Figure 2b, in order to increase the TOF, the exponential indexes of the dominating parts should be reduced: indexes of part 1 and 2 can be decreased as much as we want by increasing the free energy of reactants. The index of part 6 can be decreased by reducing the free energy of hydrogen gas. By neglecting the saturation, part 4 is absent as well. The exponential indexes of part 3 and 5 can be decreased by increasing the free energy of hydrogen gas, so the index corresponds to the activation energy of a single step, given as 14) and the coverage of the corresponding part can be written as Because the coverage of neither part can exceed 1, the maximum TOF of a reaction can be estimated by The insights of eq 17 are pretty straightforward: the maximum TOF of the reaction is determined by the step with the highest energy barrier. If the step is a hydrogenation step, then, the maximum TOF is achieved when the intermediate regarding the step and hydrogen each occupies ∼0.5 of the surface, leading to 0.25 coefficient for the first term in eq 17. On the other hand, if the step with the highest activation energy is the desorption step, the maximum TOF is achieved when the product adsorbate dominates the surface. It should be noted that eq 17 only indicates a good upper limit of TOF; the net coverage of hydrogen and another RDI can be close to but does not reach 1; a minimum amount of coverage for the empty site is always required for the reactants to adsorb. Also, the result of eq 17 does not take the adsorbate−adsorbate interaction into account; in the case that the coverage of intermediate has a significant impact on key energetic parameters, eq 17 may require further corrections. 19−23 The hypothetical reaction defined in Section 2.4 is taken as an example, and the free energy of reactants (both R and hydrogen) is tuned. TOF as a function of G H and G R is plotted and shown in Figure 4. Here, G H refers to the free energy of hydrogen per atom in the gaseous form, and G R refers to the free energy of reactant R in the gas form.
In Figure 4a, it can be observed that the TOF can vary a few orders of magnitude depending on the free energy of the reactants. Four points are marked in Figure 4a as examples. At point 1, the free energy of reactants is low (E A < E h − E R g ); thus, the reactant adsorption is relatively slow, while desorption is fast. Due to the low coverage of intermediates, the overall TOF is only 2.5 s −1 . At point 2, the free energy of reactant R is increased compared to point 1; thus, the intermediate coverage increases as well. According to the definition, the net energy barrier E A for the hydrogenation steps is 1.0 eV. While the RDS is the hydrogenation from R 1 * to R 2 *, the energy barrier actually corresponds to two hydrogenation steps from R 0 * to R 2 *. As a result, R 0 * is the rate-determining intermediate with the highest coverage (∼0.8). The hydrogen coverage is decreased due to the abundant intermediates compared to point 1, and the TOF is ∼7 s −1 . At point 3, by tuning G H , the net energy barrier for the hydrogenation steps is lowered to a minimum value of 0.9 eV and corresponds to the hydrogenation from R 1 * to R 2 *. R 1 * becomes the dominant intermediate along with hydrogen, and the TOF rises significantly to 370 s −1 . This TOF is also close to the maximum TOF estimated with eq 17 (∼420 s −1 ). At point 4, the dominant factor is still the hydrogenation reaction with the same overall barrier as point 3. The RDS is still the hydrogenation from R 1 * to R 2 * with a barrier of ∼0.9 eV. However, due to the increased G R , R 1 * occupies most of the surface and suppresses the hydrogenation coverage, which hinders the hydrogenation reaction, and the TOF is only 9 s −1 .
Note that the abovementioned discussion assumes that we can tune the free energy of gaseous reactants as much as we want. In practice, this is likely very unfeasible: The free energy of a gaseous reactant is exponentially correlated with its partial pressure, and the temperature of all reactants is the same. Therefore, the free energy of gaseous reactants can only be tuned in a limited range. Whether the optimal reaction conditions are viable in practice should be considered as a critical factor while developing the catalyst base.
In some previous studies regarding hydrogenation of acetylene or benzene, 22,23 it is reported that a high TOF is achieved when the corresponding intermediate and hydrogen have a comparable coverage over the surface, which agrees with the abovementioned discussion. However, in these studies, such an optimized coverage distribution is achieved mainly by adsorbate−adsorbate interactions and tuning the reaction temperature, so the free energies of the reactant and hydrogen gas are not tuned independently. Also, increasing temperature likely accelerates the reaction even if the coverage is off from the optimal distribution.
2.5. Results and Discussion: Saturation and Convert Rate. In principle, the reactants are constantly converted to products, until the free energies of gaseous reactants and products are equal. Ideally, the maximum TOF is sustained until such equilibrium is achieved. However, in practice, the TOF can decrease rapidly far before reaching equilibrium. In this part, we discuss how the TOF changes as the reaction progresses from the coverage point of view.
In Figure 2b, the factor regarding the overall reverse reaction is γ. According to its definition, at the beginning of the reaction without any product γ = 1, as the reaction gets close to the equilibrium state, γ → ∞, and thus, r → 0. However, γ is only significant when the reaction is close to equilibrium. Taking kT = 0.04 eV as an example, when the free energy of the product is 0.1 eV lower than that of reactants, γ ∼ 1.09 and does not significantly change the total TOF.
In the presence of the product, part 4 in Figure 2b should be taken into consideration. According to the previous discussion, let us assume that we already tune the free energy of reactants such that in the absence of the product, the TOF is maximized. In the case that the hydrogenation reactions are ratedetermining, at maximum TOF, part 3 and part 6 should be the dominant parts. Assuming an initial hydrogenation coverage of ∼0.5, the coverage of part 3 and part 4 can be written as Comparing the two coverages, under the following conditions, part 4 is dominant over part 3 and thus becomes the ratedetermining factor and significantly affects the TOF In the case that the desorption of the product is ratedetermining, at the maximum TOF, part 5 should be dominant. Through a similar approach, part 4 is dominant over part 5 and becomes rate-determining under the following conditions In the abovementioned expressions, some constants that do not change the order of magnitude are neglected. The abovementioned expressions indicate that when the free energy difference between reactants and products is lower than a threshold, saturation becomes dominant and determines the overall TOF. Note that here, the overall reverse reaction rate may still be trivial compared to the forward rate. However, the products saturate the surface by readsorption and are trapped as certain stable intermediates. This is the key difference of this mechanism from the overall reverse reaction.
To verify the abovementioned conclusions, we recalculate the TOF for point 3 in Figure 4 with various product free energies; the result is plotted in Figure 5. Remind that at point 3 in Figure  4, the dominant factor is the hydrogenation reaction in the absence of the product; according to eq 20 and the energetic parameters, coverage part 4 becomes the dominant factor when E R g − E g g < 0.6 eV. In Figure 5, we do see a relatively stable TOF until around E R g − E P g = 0.6 eV; then, the TOF rapidly decreases. Again, the factor γ ∼ 1 at this point, so this dramatic TOF drop is not due to the overall reverse reaction. In this specific case, the product saturates the surface by forming R 2 * and hinders the overall TOF.
Due to saturation, the reaction may have a good initial TOF but may end up with a poor convert rate in some cases. Here, the convert rate refers to the percentage of the reactant converted to the product. It is a critical parameter in industrial production and is usually estimated by the ratio of partial pressures between gasphase reactants and products. In order to interpret the convert rate via microkinetics, we estimate the convert rate by the free energy difference, as pressure can be derived from free energy under given reaction conditions. Obviously, the upper limit of the convert rate is achieved when the product has the same free energy as that of the reactant under reaction conditions.
In practice, we can assume a cutoff rate r co < r max : when the TOF drops to r co , the reaction is too slow and can be considered as stopped. Now, we can substitute the cutoff rate into Figure 2b and solve for the free energy of the product when the reaction is "stopped" by assuming the saturation part (part 4) being dominant (∼1). The results are given as If the reaction is initially carried out with optimized conditions to achieve a maximum TOF, hydrogen coverage can be written as a function of E A . If we further assume that γ is still close to 1 when the reaction "stops", we can rewrite eq 2 as If the convert rate refers to the rate from R to P, we may rewrite eq 23 as By substituting the corresponding parameters (temperature, enthalpy, entropy, etc.) to the left side of eq 23 or eq 24, the partial pressures of reactant and product can be solved, and the convert rate can be estimated. Note that eq 23 and eq 24 are only valid for a reaction that has a relatively high TOF but almost stops before the overall reverse reaction becomes significant and gives the maximum possible convert rate based on optimized reaction conditions. Sabatier's principle indicates that the catalyst should bind the intermediates neither too strong nor too weak. As the second and third terms in abovementioned expressions are related to the most stable/unstable states in the reaction, it actually indicates that in order to increase the convert rate, the catalyst should also minimize the free energy variances among intermediates and transition states.
In a previous model, under the reaction conditions at point 3 in Figure 4 and by taking the cutoff rate as 1 s −1 , according to eq 23, the free energy difference between the reactant and product is estimated as ∼0.45 eV, which agrees with Figure 5.
2.6. Hydrogenation Reaction Model with a Decomposition Step. A number of practical reactions of interest, such as reversed water−gas shift 32,33 and methanol synthesis, 34,36 include at least one surface decomposition reaction step, which leads to a side product. This introduces more coupling correlations in the microkinetic model, which should be addressed to expand the generality of the method. In this section, we modify the hydrogenation model to make it compatible with the reaction with one decomposition step and show that the previous results are still applicable with necessary changes.
Consider the following reaction The elementary steps involve three series of hydrogenation reactions (marked as A, B, and C), linked by a surface decomposition reaction, which is shown in Figure 6. Due to the decomposition step, there is a main product P and a side product S involved in this reaction.
Following a similar approach as we did in the previous section, we may build the microkinetic model and present it in a coverage-oriented way, as shown in Figure 7. In Figure 7, most  Microkinetics of a hydrogenation reaction with the decomposition step present in a coverage-oriented way; intermediates within the three hydrogenation series are separated. Part 5 is correlated with the decomposition step. Note that some terms in 3A are not directly related to the hydrogenation reaction. However, here, we do not distinguish them further. energetic parameters have the similar expression as in Figure 2b. The upnotes A, B, or C indicate that the corresponding energetic parameters are within the specific series of hydrogenation steps (e.g., E l A refers to the lowest free energy of all states in hydrogen series A). E d i refers to the free energy of the transition state for the decomposition step. Detailed definition for energetic parameters in Figure 7 is given in the Supporting Information.
There are two important differences in Figure 7 compared with Figure 2b: First, an extra part corresponding to the decomposition step is involved (part 5); the part is correlated with the empty site coverage because the decomposition reaction requires one empty site. In the case that the decomposition step takes place simultaneously with hydrogenation (H* + A A * ↔ B 0 * + C 0 *), we only need to replace the θ # with θ H for coverage part 2, 3, and 5 in Figure 7. Second, due to the decomposition reaction, there are more terms in the expression for the hydrogenation reaction series A. Most of these terms correspond to different combinations of one part regarding θ n B and the one regarding θ n C , and some combinations result in second-order terms. These terms also propagate to the expression of the reactant desorption part for empty site coverage (part 2).
Following a similar approach discussed in Section 2.4, the maximum TOF of the hydrogenation reaction with the additional decomposition step taken into consideration can be roughly estimated as The first three terms have similar meanings as those in eq 17, except for that there are two products, so both their desorption barriers must be considered. The last term in the abovementioned expression corresponds to when the decomposition step is rate-determining. In this case, part 5 in Figure 7 is dominant. The coefficient 0.25 indicates that the maximum TOF in this case is achieved when the corresponding intermediate and empty sites both have a coverage close to 0.5. Also, in the case that adsorbate−adsorbate interactions are significant and must be considered, eq 25 may subject to corrections.
For saturation, as shown in Figure 7, three parts (4A, 4B, and 4C) are involved. Depending on the properties of the products and the reaction conditions, they should be carefully examined  Figure  5. The black area is where the free energy of products is higher than that of reactants, and therefore, the TOF is negative.  with the free energy of both P and S to estimate the significance of saturation.
For the convert rate, due to the existence of two products, a number of solutions can be obtained when the cutoff TOF is substituted. Keep in mind that in practice, the free energy of the two products should be correlated in some way. For example, the free energies of product P and S at 1 bar and reaction temperature are 0.5 and 0.4 eV, respectively. In the case that the reaction begins without P or S, during the reaction, their free energy should always satisfy G P − G S = 0.1 eV. With the correlation, the convert rate can be estimated in a similar way as discussed in Section 2.5.

RESULTS AND DISCUSSION
Consider the following hypothetical hydrogenation reaction with a decomposition step. Again, the reaction is carried out with kT = 0.04 eV The energetic curve of the reaction is presented in Figure 8a, and the detailed energetic parameters of the reaction are listed in the Supporting Information. As discussed in Section 2.5, in order to maximize the TOF, the hydrogen coverage should be ∼0.5 and dominant over the empty site; thus, we set E H to be 0.1 eV to first guarantee its dominance over the empty site, and G R is to be found later to achieve a balance coverage of hydrogen and RDI.
Based on these parameters of elementary steps, we see that the step with the highest activation energy is the hydrogenation from B 0 * to B 1 *. According to a previous discussion, the maximum TOF at the beginning of the reaction is achieved when B 0 * and H* both have a coverage of ∼0.5. By substituting the abovementioned conditions and solving the coverage figure in Figure 7, G H = 0.1 eV and G R = −0.2 eV result in a nearoptimized TOF, which is estimated to be ∼420 s −1 .
The abovementioned TOF can be maintained until any term that is related to product free energy becomes significant. By substituting the abovementioned results back to the terms with product free energy, we get the conditions under which the abovementioned TOF can be maintained, which are expressed as In order to verify the abovementioned results, we simulated the TOF of the abovementioned reaction with G H = 0.1 eV, G R = −0.2 eV, and various product free energies. In the absence of products, the TOF turns out to be ∼370 s −1 , which is close to our estimated number. In the presence of products, the results are plotted in Figure 8b, which agree with the abovementioned conditions as well.
Based on the abovementioned results, we see that the analytical approach for a pure hydrogenation reaction is still valid with a decomposition step as long as we deal with the extra step and the free energy of two products. The approach can be applied to hydrogenation reactions with more decomposition steps. In the following sections, we will stick to the pure hydrogenation reaction due to its simplicity.

Branched Hydrogenation Reaction Model.
In previous models, it is assumed that only a single reaction is taking place over the catalyst surface with a specified product. In reality, a reaction can have multiple branches, which correspond to different reaction paths and may lead to side product(s). For example, the hydrogenation of acetylene may lead to either ethylene or ethane; 23 the hydrogenation of benzene has various reaction paths with different intermediates; 22 carbon dioxide hydrogenation to methanol also involves at least two reaction paths. 34,36 An immediate and important problem is then the selectivity of paths or products. In this section, we set up a branched pure hydrogenation model, discuss how multiple branches are coupled, and illustrate the factors determining the selectivity with the coverage-oriented approach.
Consider the following two reactions The two reactions share some elementary steps, which are shown in Figure 9.
In this model, the TOF values for P A and P B are different; let them be r A and r B , respectively. Following a similar approach as discussed in Section 2.1, the microkinetics can be expressed in the coverage-oriented way, as shown in Figure 10a. The energetic parameters applied in Figure 10a are defined in a similar way as those in Figures 2b and 7.
In Figure 10a, part 4A has a more complicated expression compared with Figure 2b, as the three terms have the coverage required to compensate reverse reactions of hydrogenation steps, desorption step, and saturation. Note that here, we do not classify them further but put them in a single part, as the two equivalent expressions are derived from the two paths separately. The parameters characterizing the two paths can be correlated by the following Ä To evaluate the selectivity at a given point, we can directly solve the microkinetics in Figure 10a by plugging in the reaction conditions at the point. In the case that we need to evaluate the final product selectivity, the TOF of both products should be evaluated from the beginning of the reaction. According to the TOF and the properties of the reactants/products, their free energy shall be updated as the reaction carries on. By updating the TOF and free energy alternately until both the TOFs is below a cutoff TOF, the amount of each product can thus be determined.
A special case of the abovementioned model is that if the two reactions are actually two reaction paths that lead to the same product, according to (eq 26), in the case that the hydrogenation steps are rate-determining for both paths, we can immediately write down the TOF ratio between the two reaction paths without loss of generality, let us assume that r A dominates over r B , so the net TOF can be approximated as r A ; then, the coverage figure can be plotted as Figure 10b. Some interesting phenomena due to the coupling among the branches/paths may be observed, which we will be discussing in the next section. In this model, the approach to model one reaction with multiple branches is introduced. A similar model can be built for multiple independent reactions. Similar to the coupling based on the coverage of part 4A, in case of multiple independent reactions, they can be coupled by the coverage of part 2.
3.2. Results and Discussion: Poisoning and Phantom Path. Consider the following hydrogenation reaction The microkinetics of the reaction is presented in Figure 11a. From the figure, there are actually two reaction paths (A and B). The free energy of gaseous reactants is G R = −0.2 eV and G H = 0.1 eV, which is close to the optimal TOF. Also, kT = 0.04 eV.
Detailed thermodynamic parameters can be found in the Supporting Information.
The simulated TOF from both paths as a function of the product free energy is plotted in Figure 11b. As predicted by eq 27, the ratio of the TOF from the two paths is a constant (∼e 2.5 ), and the overall TOF can be approximated by that contributed by path A. Furthermore, the TOF decays rapidly as the product free energy goes above −0.3 eV; this is due to the saturation, as discussed in Section 2.4, and in this case, it is mainly due to B 2 * involved in part 4B of Figure 10b. Now let us consider two interesting cases, and before that, we repeated the abovementioned simulation with only one path; the results are plotted as dashed lines in Figure 11b. Besides, the simulated coverage of each intermediate at a specific point is plotted as well.
In the first case, assume that through a theoretical study, path A is identified, but path B is ignored. In Figure 11, we see that the TOF with only path A is significantly higher. The reason is that by ignoring path B, we ignored the coverage part 3B, 4B, and 5B Figure 11. (a) Energy curve of the branched hydrogenation reaction. (b) TOF dependence on the product free energy. The TOF of path A and path B, the case with only path A or path B, and the net TOF are plotted by different lines. Note that the TOF with both path A and B almost overlaps with that with only path A in the graph because the contribution from path B is negligible. (c) Intermediate coverage at G P = −0.4 eV when considering both reaction paths or only path B. The coverage distribution is almost the same in the two cases, due to the fact that intermediates in path A are relatively unstable. Therefore, the existence of path A can hardly be confirmed by identifying intermediates experimentally.
in Figure 10b. In this example, intermediate B 1 is too stable and blocks the overall reaction even in the absence of products. This indicates that some reaction paths may poison the overall reaction, and emphasizing one path with high activity while ignoring other paths can lead to TOF overestimation.
In the second case, on the other hand, path B is identified but not path A. If we repeat the simulation with only the intermediates and steps in path B, the intermediate distribution is very similar to that with both paths, due to the very low coverage of intermediates in path A. In practice, this may make path A very difficult to identify experimentally. However, as shown in Figure 11, path A is actually the main contributing path.
In previous studies, usually, the coverage and contribution of a reaction path are correlated: if intermediates of a specific path have high coverage, the path is considered as the primary one. On the contrary, if intermediates associated with a path all have low coverage, the path is considered insignificant. The coverage of intermediates can be estimated by experimental methods (spectroscopy, e.g.). However, based on abovementioned discussion, we see that a path can be contributing but with low-coverage intermediates, which we refer to as a "phantom" path; it can dominate the surface coverage but contribute little to the overall TOF, and we refer to it as a "poisoning" path. Due to the possible existence of "phantom" and "poisoning" paths, predicting and identifying possible reaction paths through a theoretical approach are of critical importance.
3.3. Pure Hydrogenation Reaction over Multiple Strongly Correlated Sites. In all the previous models, all reactions are carried out over a single type of site. However, in industrial production, complex catalysts are usually applied; 31−35 thus, various sites can participate in the reaction process. Even when a simple catalyst is applied, it still involves various sites composed of different facets, edges, or even defects. 37,38 In the case that the multiple sites function independently, the abovementioned approach can be applied for each type of sites. Otherwise, some necessary modifications are required to take the coupling of different sites into consideration.
In this section, the coverage distribution over multiple strongly correlated sites for a pure hydrogenation reaction is discussed. Here, the term "strongly correlated" means that the sites are well mixed, and intermediate exchange among them is much faster than reaction steps, such that the diffusion among sites is considered to be at equilibrium. It is an extreme case as opposed to the case where sites of different types are completely "isolated".
In order to interpret the microkinetics over strongly correlated sites, the strong correlation in both coverage and reaction rate must first be discussed.
Assume that there are two types of sites θ and φ, which are well mixed on the catalyst surface. The percentages of the two sites are P θ and P φ , so P θ + P φ = 1. The coverage of intermediates on the two sites is represented by θ and φ.
Due to the assumption that the diffusion among different sites is at equilibrium, the coverage of any intermediate X on site θ and φ must satisfy the following (detailed deduction can be found in the Supporting Information).
In the abovementioned equation, G X with upnote θ or φ indicates the free energy of adsorbate X on either site θ or φ, respectively. The coefficient λ is the ratio of the empty site on site θ and φ; its chemical physics insight is the free energy contribution difference on the two sites due to configurational entropy contribution. The insight of eq 28 is that the free energy of any intermediate on different strongly correlated sites must be equal, considering the configurational entropy contribution. 16,18 The hydrogenation reaction in Section 2.1 but this time on the correlated sites θ and φ with a total TOF of r is discussed. The elementary steps are shown in Figure 12. Following a similar approach while taking eq 28 into consideration, we can find that the net desorption rate from the two sites satisfies r r / ds ds λ = θ ϕ (29) In the abovementioned equation, r ds θ is the net desorption rate from site θ.
For the hydrogenation steps, due to the two different sites, each step now has four possible site combinations characterized by different energetic parameters, as shown in Figure 12. For any specified step, the net contribution from the site combination should meet  (30) In eq 30, E TS,n is the energy of the transition state for the hydrogenation step; its upnote indicates its site combination, as in Figure 12. The insights of eq 29 and eq 30 are also straightforward: when there are multiple site combinations for an elementary step, the contribution from each combination is associated with the free energy of the transition state, where the configurational entropy contribution must be taken into account. The combination with the lowest barrier contributes more than others. Detailed deduction for eq 29 and eq 30 can be found in the Supporting Information.
Based on eq 30, we may define the equivalent transition-state energy E TS,n for any hydrogenation step as Now that the strong correlation of coverage and reaction rates is given by eqs 28−30, we are able to derive the expression for the microkinetics and presented it in Figure 13. Compared with Figure 2b, an extra parameter λ is involved, which can be solved with one more restraint regarding the percentage of the two Figure 12. Elementary steps of a hydrogenation reaction over two sites θ and φ. Marked activation energies refer to a forward reaction. Note that diffusion among different site types is not shown here.
sites. While important parameters such as maximum TOF and convert rate cannot be written in a simple form due to the complicated coupling, we are still able to interpret some interesting mechanisms with a coverage-oriented microkinetic model, as discussed in the next section.
3.4. Results and Discussion: the Assembly Line Reaction Mode. Consider the reaction discussed in Section 2.4, but it takes place on the catalyst with two different sites θ and φ this time, and each type accounts for a half of the total sites. The energetic parameters for the reaction step on site θ are the same as Section 2.4.
The energetic curve of the reaction is given in Figure 14, and detailed information can be found in the Supporting Information. In the figure, the intermediates on two sites and the hydrogenation with the hydrogen initially adsorbing on the two sites are marked with different colors, for example, the red curve connecting red flat lines R 0 * + H 2 and R 1 * + 1/2H 2 indicates the reaction step with R 0 and hydrogen both initially on site θ. For the free energy of the reactants, these values from point 3 in Figure 4a are applied, as they give high TOF on site θ alone (G H = 0.2 eV and G R = 0.1 eV). By substituting these parameters and solving the reaction numerically, the TOF turns out to be r ∼ 27,000 s −1 and λ ∼ 0.013.
The TOF in this case is around two orders larger than that in Section 2.4, and it is not due to the activity of site φ alone. As we can see, the reaction step from R 0 to R 1 on site φ has an activation energy of ∼0.9 eV; according to eq 17, the maximum TOF with only site φ is also around 420 s −1 . The TOF can only be interpreted by the synergy of the two coupled sites.
From Figure 14, we see that for the hydrogenation from R 1 to R 2 , the combination of R 1 × and H* has the lowest transition-state free energy, lower than that of the combination of R 1 * and H*.  Because this step is the RDS when there is only site θ, lowering its reaction energy barrier can significantly increase the TOF, which accounts for the TOF difference from that in Section 2.4. If we further examine the contribution for each step according to eq 29 and eq 30, the adsorption step is mainly done by site φ, as it has more empty sites and has strong binding to R 0 ; the first hydrogenation step (R 0 to R 1 ) is mostly done on site θ alone, while the second step is mainly carried out by R 1 on site φ and hydrogen on site θ; eventually, the desorption is mainly carried out on site φ. In between these steps, intermediates diffuse among sites θ and φ, reaction steps with high energy are greatly avoided, and the overall TOF is significantly increased. Here, we refer to this mechanism as "the assembly line reaction mode", as it is similar to the corresponding industrial production mode, where the product has to be transferred among multiple sites to achieve high efficiency. Such a mechanism is also reported in a few previous studies. 39 Next, we repeat the simulation in the presence of products. The TOF as a function of product free energy is plotted in Figure  15. For comparison, the TOFs with only site θ or φ are plotted as well. From Figure 15, it can be observed that the TOF on both sites has the following trend: initially, it is relatively stable, until G P ∼ −0.1 eV, where it begins to decrease rapidly. This is due to the saturation on site θ (part 4θ in Figure 13), which reduces the activity of the site. As the product free energy keeps increasing, site θ is hardly contributing due to saturation. The assembly line reaction mode breaks, and the reaction takes place on site φ alone, as if site θ is absent. Eventually, site φ gets saturated as well, and the TOF begins to decrease rapidly again. The two inflection points in the curve match the saturation point of individual sites. The abovementioned result indicates that while the assembly line reaction mode can potentially increase the maximum TOF, due to the saturation, the convert rate is still mostly dependent on the site which yields the highest convert rate when applied alone.

CONCLUSIONS
In this work, the microkinetic model for heterogeneously catalyzed hydrogenation reactions is built, which is compatible with the industrial production where multiple main and side products and complex site features are involved.
The approach is applied on four different hydrogenation reaction models to comprehend industrial synthesis. In Section 2, we show how to evaluate the maximum TOF of the desired reaction and how to achieve this TOF by optimizing the coverage of the corresponding intermediates instead of leaving it in a single "MARI" state. Besides, it indicates that the catalyst should be designed in a way such that its optimal reaction conditions are feasible in practice. The theoretical convert rate of the desired reaction on a given catalyst can also be estimated analytically. These are key parameters characterizing the production outcome. In Section 3, we extend the model to make it compatible with more reactions and show that the abovementioned conclusions hold true with necessary modifications. In Section 4, we show how to determine the selectivity or relative activity of multiple reaction paths and how a reaction path can be either "poisoning" or "phantom". This is also the reason why microkinetics are difficult to be interpreted through only experimental methods. In Section 5, the microkinetic model is extended to multiple strongly correlated sites. It is shown that involving multiple different active sites can potentially increase the TOF of a reaction through the "assembly line reaction model", but the convert rate is limited by the single type of site with the best performance. This helps to optimize the production and develop more complicated catalytic systems.
Although the reactions discussed in this work are all heterogeneously catalyzed hydrogenation reactions, the framework can be extended to simple oxidation reactions, such as the oxidization of CO, due to their similarity to simple hydrogenation reactions. Note that most oxidation reaction systems have more complicated mechanisms, and further adjustments are necessary should the framework be applied. Besides, by replacing "coverage" with its equivalent in homogeneous catalytic system, the framework may help in interpreting homogeneously catalyzed reactions as well.

SIMULATION METHOD
The steady state of a reaction refers to the state where the amount of each adsorbate on the catalyst surface does not change with time, but the overall reaction is not at equilibrium, so the reactant is still being converted to the product. In the case that the amount of reactants overwhelms that of the catalyst, the system generally reaches the steady state quickly compared to the reaction process, making the concept very useful for estimating the rate of a reaction. Mathematically, given all the energetic parameters, the coverage of all adsorbates and the reaction rates of all steps can be solved in the steady state.
According to its definition, at the steady state, for each adsorbate, its generation rate r + is equal to its consumption rate r − . Here, the generation rate is the summation of all the rates regarding an elementary step that generates the intermediate, similar for the consumption rate.
The consumption rate of an intermediate is related to its coverage: the higher its coverage, the faster it is consumed. Based on the correlation, the following algorithm is processed in a self-written C++ program to find the steady state of a reaction: Figure 15. TOF of the hydrogenation reaction as a function of product free energy on both site θ and φ and only θ or φ. G H = 0.2 eV, G R = 0.1 eV, and kT = 0.04 eV. It can be observed that the TOF rapidly decreases as site θ gets saturated and then becomes relatively stable after θ is almost deactivated. Another rapid decrease in TOF is observed as site φ gets saturated. The TOF trend at high product free energy is similar to that on site φ only and so is the convert rate.
All the energetic parameters from DFT results and reaction conditions are set. Then, the process is begun with an initial arbitrary guess for the coverage of all intermediates. The only requirements for the initial guess are that all intermediates (including unoccupied free sites) have non-zero coverage and the net coverage is unity. In each iteration step, the generation rates and consumption rates for all intermediates at the current state are first calculated. For each intermediate i, a parameter α is defined to evaluate how close it is close to the equilibrium state Apparently, if the coverage of an intermediate does not change with time, eq 6 is substituted into eq 7, and its α equals to 0. If α for all intermediates is below a threshold (10 −10 in this work), the current state is considered as the steady state, and the iteration stops. Otherwise, a random intermediate i is picked, and its coverage is modified by multiplying it with factor c. In this work, factor c is calculated in the following way c r r ( / ) In other words, if the generation rate is higher than the consumption rate for the picked intermediate, its coverage is increased to balance the generation rate and vice versa. The index 1/2 is arbitrarily chosen, as it gives a good convergence performance. The coverage of all intermediates is then normalized to keep the net coverage unity. The algorithm then proceeds to the next iteration step, until the steady state is found. The result is verified as α for all intermediates is close enough to 0 (less than 10 −10 in this work), and the consumption rates of gaseous reactants match the generation rates of gaseous products. Funding from the "Liquid Sunshine" energy system project (grant no. 211211KYSB20180020), supported by the International Partnership Program of the Chinese Academy of Sciences, is gratefully acknowledged.

Notes
The author declares no competing financial interest.