CO 2 Electrocatalyst Design Using Graph Theory

The electrochemical reduction of CO 2 (CO 2 RR) to higher-order hydrocarbons or oxygenates using low-carbon electricity offers a promising path to generate renewable fuels and chemicals; however, the low selectivity of present-day CO 2 RR catalysts toward more valuable C 3 products limits technoeconomically compelling avenues toward productization. Systematically enumerating possible intermediates and reactions to the set of possible {C 1 , C 2 , and C 3 } CO 2 RR products entails 3206 intermediates and 4506 reactions. Here we use graph theory to enumerate these possibilities comprehensively, treating intermediates as graph nodes and pathways as graph edges. C 3 products fall into two groups, group A (1-propanol, allyl alcohol, and propionaldehyde) and group B (acetone and hydroxyacetone); and we nd that an early branch reaction, three steps after C 1 -C 2 coupling, bifurcates these two groups of C 3 products. We develop a set of C 3 descriptors and screen catalysts for CO2RR to C 3 products: specically, we screen bimetallic (doped Cu) catalysts for their combination of CO binding, C 1 -C 1 coupling and C 1 -C 2 coupling. Cu and also Au- and Ag-doped Cu fulll the rst two requirements; but the set of promising C 1 C 2 coupling catalysts (Ni-doped Pb and Al-doped Pb) needed to get to C 3 is nonoverlapping with that for CO binding and C 1 -C 1 coupling. Our ndings agree with the experimental picture that Cu, while among the most productive to propanol, has been limited in F.E. to the 10-20% range; and that, to date, no single catalyst has achieved exceptional C 3 productivity. We discuss tandem catalyst designs where a rst catalyst promotes CO binding and C 1 -C 1 , and where these C 2 and C 1 intermediates can come together for coupling on a second distinct class of integrated catalysts.


Introduction
The electrochemical reduction of carbon dioxide (CO 2 RR) converts CO 2 to fuels and chemical feedstocks using renewable electricity 1 . CO 2 RR can in principle produce a vast array of possible hydrocarbons and oxygenates: the carbon, oxygen, and hydrogen atoms from carbon dioxide and water are building blocks that constitute an extensive combinatorial space. Experimentally, only a few products have been reported to have high selectivity: these are carbon monoxide 2 , formic acid 3 , methane 4 , ethylene 5 , and ethanol 6 .
Designing catalysts that go beyond these products is an important task in CO 2 RR.
Rational catalyst design requires an understanding of reaction pathways to the desired product, and also to the competing products. Jaramillo and co-workers observed 16 different products with numbers of carbon from one to three (C 1 -C 3 ) in CO 2 RR on metallic copper catalysts 7 ; yet the routes to these diverse products -such as, among C 3 products, hydroxyacetone, acetone, allyl alcohol, propionaldehyde, and 1propanol -have not been reported.
For C 2 products, different groups report distinct reaction mechanisms: for example, C 1 -C 1 coupling mechanisms have been proposed that range from OC-CO 8,9 to OC-COH 10 to OC-CHO 11,12 to OC-CHOH 12 and H 2 C-CH 2 . 13 The key intermediate species branching ethylene and ethanol pathways are also under debate: Graza et al. 11 proposed that this bifurcation occurs immediately after OC-CHO coupling forming OCCHO* (5etransferred from CO 2 ) on Cu(100). Goddard and co-workers 9 have suggested that HCCOH* (8etransferred) is the key surface species on Cu(100) before the bifurcation of ethylene and ethanol.
Calle-Vallejo and Koper 14 have reported that the key intermediate is H 2 CCHO (9etransferred).
The lack of authoritative mechanisms to C 3 products and debate over C 2 routes limits the rational design of new CO 2 RR catalysts toward desired products: only once reaction pathways are known does it become feasible to pursue strategies to decrease the energies of intermediates along the pathways toward the desired product while increasing the energies of the intermediates in the competing pathways. For example, the design of CO 2 RR catalysts with high ethylene selectivity requires knowledge of the C 1 -C 1 coupling steps and the branch reactions to other C 2 products. This can inform the engineering of active sites that favour C 1 -C 1 coupling and that limit branching to other products.
The number of possible reaction steps and intermediates in the CO 2 RR up to C 3 species is large ( Figure   1a): there are 20 elementary steps in the pathway to 1-propanol, including 18 hydrogenation steps and 2 coupling steps. There are many possible choices for each step: for example, in the simple case from CO to CH 4 , there are ten different pathways shown in Figure 1b. A comprehensive analysis of CO 2 RR must enumerate systematically the possible intermediates and pathways. Further, a quantitative approach will be required: rather than reject certain pathways as unlikely, we propose a quantitative analysis of the energetics of each possible pathway.
We turned to graph theory 15 to enumerate the possibilities systematically, density functional theory (DFT) to calculate reaction energies, and a steady-state solver to evaluate the most favorable pathways. In a reaction network (Figure 1b), species (reactant, intermediates, and products) are represented as nodes; and reactions (edges) convert one species into another. Here we investigate CO 2 RR to products up to and including C 3 on both Cu(111) and Cu(100).

Results And Discussions
Graph theory framework We began from CO species since the activation of CO 2 to CO has been systematically investigated before 16,17 . We rst obtained the possible chemical formulas up to the maximum carbon, oxygen and hydrogen atoms subject to C≤3, H≤8 and O≤3 (Figure 2a We then calculated the formation energies of these species on Cu(111) and Cu(100). We chose the most unsaturated atom as the binding atom, and calculated the adsorption of the unsaturated species on the possible sites, including atop, bridge, hcp hollow and fcc hollow on Cu(111) and atop, bridge, and fourfold hollow on Cu(100) (Figure 2b), using an adsorption vector algorithm 18 . Following optimisation using DFT, we checked the connectivity matrix of the intermediates to ensure that there occurred neither dissociation nor reconstruction, and the dissociated or reconstructed species are removed from the species list. The effects of potential, solvent, and ions on the species formation energies are important, and there are several methods in literature to consider these effects, such as the ab initio molecular dynamics with biased sampling method 9 , the grand canonical quantum mechanics (GCQM) method 8,19 , charged water model using global optimization method 20 , and implicit solvent model 21 Figure S1. The most common reaction types are hydrogenation, C-C coupling, and water formation, in that order ( Figure 2c). We checked the species pairwise to enumerate the possible reactions. For example, if species X and XH had the same structures except for one extra hydrogen, we added the reaction X* + H + « XH* to the reaction list. Using this algorithm, we consider the possible hydrogenation, coupling, and dehydroxylation reactions in CO 2  With the formation energies of the possible species and the comprehensive reaction list, we then obtain the reaction energies ( Figure 2d). We used the reaction energies in one pathway to describe the favourability of this pathway (more details are found in the Supplymentary Information). The effect of applied potential is considered with the computational hydrogen electrode, 22 and pH value is taken into account via the concentration of the reactant (proton source) in the steady-state solver ( Figure 2d).
In sum, the reaction pathway search framework ( Figure 2) generates structures of the species, calculates the adsorption structures of the intermediates, obtains the reaction energies, and builds the reaction graph using Catkit 18 and an ab initio simulation program (VASP) [23][24][25][26] .
Reaction graph of CO 2 RR Figure 3a shows the reaction graph for CO 2 RR, including the possible intermediates and products up to three carbons and their corresponding elementary steps. The number of C 1 species is 13 (Figure 3a), 55 for C 2 species, 385 for C 3 species, and 5 species without carbon. The species react with protons, forming other species with the same number of carbon atoms (black edges for hydrogenation and green edges for dehydroxylation in Figure 3a), or couple with other species to produce species with a higher number of carbons (red edges in Figure 3a). Except for methane, the products in Figure 3a are linked with more than one edge, suggesting that each product can be produced via several possible routes.
Using the simple path algorithm 27 , we obtained the number of possible pathways for several products (Table S1): for example, there exist 10 distinct pathways from CO to methane as shown in Figure 1b. Ethylene has 4744 different pathways from CO that are at least one intermediate different, while for CO to 1-propanol the possible pathways are more than 27 million. These make it infeasible to enumerate pathways manually.
In graph theory, the degree of the nodes refers to the number of edges connected to the nodes. This corresponds to the number of possible reactions in which a given species is implicated. The average degree of all nodes is ~13.2, and the maximum degrees are CH 2 O and CH 3 O, with 84 edges connected, suggesting that these two species are the most important in the reaction network from the point of view of connectivity. The high-degree species are C 1 species, since these species participate in coupling to form C 2 and C 3 species.

Formation energies of intermediates
The energy distributions of the intermediates are shown as a function of the number of electrons transferred in Figure 3b. The formation energy is referenced to the free energies of CO 2 , H 2 O, and H 2 , which are the three reactants of CO 2 RR. Thus, the formation energy of each intermediate in Figure 3b indicates the reaction energies for the step from CO 2 and H + /H 2 O to this intermediate. The number of electrons transferred per carbon atom indicates the degree of hydrogenation. In general, the formation energies increase rst, and then decrease with further increase in the number of electrons transferred per carbon atom. Thus, the reaction energies of the rst few steps are endothermic at low applied potential, indicating that the rate-determining step of CO (2) RR is in the rst few steps. When a negative potential is applied, the formation energies of the intermediates shift by the number of electrons transferred × potential: thus the overall reaction energies of intermediates with an electron transferred become more favourable. These results agree with the high overpotential due to the endothermic steps in the initial stages of CO 2 RR. Favouring the rst few steps using active site engineering enables the design of CO 2 RR catalysts with low overpotential.
Figure 3b also shows that the formation energies of intermediates on Cu(100) are more stable than those on Cu(111), indicating that Cu(100) has a greater driving force for CO 2 reduction and is thus the more likely active site for CO 2 RR 28 . The most energetically favourable pathway is in Figure 3b: it is the pathway from one intermediate to the most stable intermediate with one more electron transferred reachable in the graph. However, if the energies of two intermediates are similar, it is di cult to determine the one most likely in energy space; instead, it is necessary to obtain the most representative pathway(s) in rate space.
The formation energy distribution for intermediates having with 0-3 carbons is provided in Figure 3c. The median and the range of C 2 intermediate formation energy is lower than in the case of C 1 intermediate, and this is true both on Cu(111) and Cu(100). This suggests that C 2 is energetically favoured and is the main product of CO 2 RR. Regarding C 2 to C 3 , the driving force is much lower on both Cu(111) and Cu(100) (Figure 3c), suggesting that it will be hard to make C 3 as the main product of CO 2 RR on pure Cu catalysts.
Experiments also show that Cu-based catalysts achieve high C 2 productivities in CO 2 RR 5 , and lower C 3 productivities 7 . C 2 Reaction mechanisms C 1 -C 1 coupling is the key step for the selectivity of the C 2 products 29 . Many different C 1 -C 1 coupling reaction mechanisms are proposed including OC-CO 8,9 , OC-COH 10 , OC-CHO 11,12 , OC-CHOH 12 , H 2 C-CH 2 13 .
However, no conclusion has been drawn regarding the dominant elementary step for C 1 -C 1 coupling. In our reaction graph, we considered all 55 C 1 -C 1 possible coupling reactions on both Cu(111) and Cu(100).
With the steady-state solver, we obtained the rates for the elementary steps up to C 2 including the C 1 -C 1 coupling steps at different applied potentials from 0 to -1.5 V vs. standard hydrogen electrode (SHE).
Despite the large number of possible coupling reactions, there are only a few dominant C 1 -C 1 coupling reactions. Figure 4a and 4b show the main C 1 -C 1 coupling reactions and their likelihood among the C 2 formation steps at neutral pH and applied potential from 0 to -1.5 V vs. SHE. The C 1 -C 1 coupling mechanisms change with applied potentials and facets due to changes in surface coverage. On Cu(100), OC-CO is the only possible pathway at 0 V, in agreement with the reaction step proposed in prior literature 8,9 . With an increase of applied potential, the hydrogenation steps become more favourable: thus, CO* is likely to be further hydrogenated. Therefore, the coverage of HCO* increases and the rate of OC-HCO coupling rises. In the applied potential region higher than -0.5 V, OC-HCO is the main reaction mechanism, a nding that has been reported by Nørskov and Head-Gordon 11,12 . The H 2 CO-H 2 CO and COH-COH coupling steps are also possible at high applied potential; however, these applied potential regions correspond to the formation of methane and hydrogen 30 . On Cu(111), CO-HCO is the main coupling step at potentials ranging from 0 V to -0.75 V, while CO-COH is also possible with a low likelihood ( Figure 4b). HCO-HCO coupling takes the place of CO-HCO at more negative applied potentials.
The ndings reported above span coupling mechanism in prior reports; and they conclude that there is not a single unique coupling mechanism. Instead, different coupling steps dominate at different conditions and facets.
We then went on to evaluate the set of paths to C 2 products. We show in Figure 4c the reaction mechanisms on Cu(100) to ethanol and ethylene, two main C 2 products for Cu-based catalysts of particular interest experimentally. 6,31 These are shown at -0.5 V and neutral pH. The reaction starts from CO-HCO coupling, forming OCCHO*, the most likely coupling mechanism at this condition (Figure 4a). The next steps are the hydrogenation of OCCHO* to OHCCHO* and then HOCHCHO*. The hydroxyl group in HOCHCHO* is then removed, giving HCCHO* and water. The removal of the second hydroxyl group happens at the intermediate HCCHOH*, and after the dehydroxylation reaction, HCCH* proceeds to ethylene with two more hydrogenation steps. This reaction pathway explains 49% of the ethylene formed.
As an alternative next step following the dehydroxylation step, HCCHOH* can be hydrogenated to H 2 CCHOH* (Figure 4c). There are two main hydrogenation steps of H 2 CCHOH* to H 2 CCH 2 OH* and H 3 CCHOH*, both leading to the formation of ethanol. Interestingly, we found that the H 2 CCH 2 OH* can also form ethylene, similar to the mechanism reported by Koper and co-worker 14 . Therefore, more than one branch reaction exists between ethylene and ethanol on Cu(100) at U=-0.5 V vs. SHE and pH=7, and the HCCHOH* and H 2 CCH 2 OH* are the two key intermediates for the selectivity of these two products.
The reaction pathways to all eight C 2 products on both Cu(100) and Cu(111) are listed in Table S3-S18.

C 3 Reaction mechanisms
As shown in Table S1, there exist 27 million different pathways from CO to 1-propanol. More carbon atoms mean more isomers due to the possible arrangements of carbon, hydrogen, and oxygen. We sought to enumerate the possible reaction pathways to C 3 products and explore the highest-rate paths.
The formation of C 3 species starts with C 1 -C 2 coupling 32-34 . There are 700 coupling reactions for C 1 -C 2 coupling due to the large numbers of possible C 2 species. We considered the coupling steps along with their ensuing reaction pathways to C 3 products. Figure 5a shows the reaction mechanism to ve products reported experimentally 7 -1-propanol, allyl alcohol, propionaldehyde, acetone, and hydroxyacetone -at U=-0.5 V and neutral pH; the reaction mechanisms to all ten C 3 products are found in Supplementary Information (Cu(100) in Table S19-S28 and Cu(111) in Table S29-S38). The most favourable C 1 -C 2 coupling mechanism at -0.5 V and neutral pH is OC-OCCHO coupling on Cu(100) (Figure 5a). Three hydrogenation steps follow the formation of OCC(O)CHO* (IM3_1 in Figure 5a Regarding the group A products, the rst hydroxyl group removed in the C 3 species is on the middle carbon ( Figure 5a). OHCCCHOH (IM3A_2) is then reduced to OHCCHCH 2 OH* (IM3A_4) before the removal of the second hydroxyl group (Figure 5), forming OHCCHCH 2 * (IM3A_5). This is followed by two possibilities: hydrogenation of CH 2 -to CH 3 -forming propionaldehyde, and hydrogenation of ketone group to the hydroxyl group. The branch between allyl alcohol and 1-propanol occurs at HOCHCHCH 2 * (IM3A_7); then the 1-propanol is formed with three hydrogenation steps. For group B products, the bifurcation of hydroxyacetone and acetone happens at HOCH 2 C(O)CH 2 * (IM3B_4) (Figure 5a).
Interestingly, previous experiments 7 have suggested that partial current densities of these ve C 3 products in CO 2 RR also fall into two groups at pH=6.8 (Figure 5b): 1-propanol has the highest current density, followed closely by allyl alcohol and then propionaldehyde. The current densities to these three products (group A) are generally one order of magnitude higher than to the other two products, acetone and hydroxyacetone (group B). Furthermore, the branch reactions of 1-propanol and other C 3 products occurs in the following orders in the reaction process: {acetone, hydroxyacetone}, propionaldehyde, and allyl alcohol. The current density ordered from lowest to highest proceeds acetone ≈ hydroxyacetone << propionaldehyde < allyl alcohol < 1-propanol.
In sum, there are two main branch reactions in the C 3 pathways: one at the initial stage after C 1 -C 2 coupling, OHCC(O)CHOH* (IM3_3), leading two groups of products; and the other one occurs at the very end of the reaction process, OHCCHCH 2 * (IM3A_5) and HOCHCHCH 2 * (IM3A_7), and these give different products.
We also obtained information on reaction pathways to other C 3 products, including ones that have not been detected experimentally. Obtaining such results offers inputs into catalyst design beyond traditional products already seen in CO 2 RR. For example, 1,2-propylene glycol is used in the production of polymers and food processing, and it has not been reported as a direct product of CO 2 RR. Our results (Figure 5a Table S20 and Cu(111) in Table S30), propane (Table S24 and S34), 2-propanol (Table S25 and S35), glycerol (Table S27 and S37).
Catalyst design for C 3 products The framework presented herein allows systematic screening of possible pathways; we acknowledge that the use of reaction energies along the pathway as a surrogate for the favourability of a pathway is an approximation, and we recognize that, in future works, the graph theoretic enumeration can bene cially be combined with (for example) a high-throughput transition state searching method, kinetic Monte Carlo, high-level quantum chemistry, or machine learning.
With that caveat, we evaluate the implications of Figure 6 -a plot that captures CO binding (a), C 1 -C 1 coupling (b), and C 1 -C 2 coupling (c).
Optimal CO binding is reached near -0.67 V 35 . Not surprisingly, Cu (both pure and doped) comes closest to ful lling this condition (Figure 6a). Au and Al include some doping options such as Pd-doped Au, Pddoped Al, Ir-doped Al, and Rh-doped Al that come within 0.2 eV of the optimum.
Favourable C 1 -C 1 coupling is the next prerequisite along the ultimate pathway towards C 3 products. Here the most strongly negative reaction energy in Figure 6b is the most desirable. Cu doped with Ag and with Au is promising here, and Al is particularly favourable, as Cu-doped Ag and Pb-doped Ag.
To go the rest of the way to C 3 products, a favourable reaction energy of C 1 -C 2 coupling is further required. Here again the most strongly negative reaction energy is desired. Ni-doped Pb performs exceptionally well. Pure Cu is a reasonable candidate but not strongly favourable. Ag doped with Al, Au, and Pb, shows promise, but the C 1 -C 2 coupling energy is not su ciently strong as to favour C 3 formation overwhelmingly compared to Ni-doped Pb and Al-doped Pb.
Overall these ndings are consistent with the observation that, in prior reports, Cu to propanol has been among the most productive, but its Faradaic e ciency has been limited to the 10-20% range; and that, to date, no single-material nor single-material-doped catalyst has achieved exceptional C 3 productivity.
Intriguingly, the nding that doping one catalyst with a single dopant does not simultaneously enable optimal CO binding; and highly reactive C 1 -C 1 coupling; and C 1 -C 2 coupling; suggests that a new strategy is needed to achieve productive C 3 generation.
One promising solution is tandem catalysts that combine two or more catalysts: a rst catalyst (such as Cu) is in charge of CO* intermediate formation and C 1 -C 1 coupling; and a second catalyst is speci cally designed for C 1 -C 2 coupling (in Figure 6, a candidate is Ni-doped Pb, though doped Au variants also suggest promise).

Conclusions
We enumerated herein the intermediates in CO 2 RR up to C 3 H 8 O 3 on both Cu(111) and Cu(100). The most favourable reaction pathways were studied as a function of pH and applied potential. We obtained the intermediates in CO 2 RR up to C 3 with the optimized structures. This offers candidates for experiment studies related to intermediates, such as in-situ Raman and infrared spectroscopy. We enumerated possible reaction pathways to 20 products, and 1144 hydrogenation reactions, 755 coupling reactions, and 354 dehydroxylation reactions were found. We found that the C We report systematically, for the rst time, the suite of reaction pathways available to C 3 products. We found there exists two groups of products bifurcate at the early stage after C 1 -C 2 coupling, namely group A (propionaldehyde, allyl alcohol, 1-propanol) and group B (acetone and hydroxyacetone). Previous experimental results showed that products in the same groups have similar partial current densities, and the partial current densities in the case of group A are an order of magnitude higher than those of group B.
We also obtained the reaction pathways to products not previously reported. These results offer the possibility to tune selectivity toward new chemicals and to design catalysts beyond those whose products are already known. The framework can be extended to more complicated products with higher numbers of carbon in CO 2 RR and beyond. Based on these understanding, we offered designs of catalysts toward C 3 products.

Computational methods
Density functional theory calculations. In this work, the density functional theory (DFT) calculations were carried out with a periodic slab model using the Vienna ab initio simulation program [23][24][25][26] (VASP) (https://www.vasp.at/). Detailed theoretical methods are found in the Supplementary Information.

Declarations
Data availability The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.   Reaction graph and formation energy distribution. a. Reaction graph from C1 to C3. Each node represents one reactant, intermediate or products; and each edge represents one reaction. Red, purple, and blue nodes stand for species with one, two, and three carbons, respectively. Coupling, hydrogenation, and dehydroxylation reactions are shown in red, black, and green, respectively. Formation energy distribution of all the species on Cu(100) and C(111) according to b. the number of electron transfers per carbon atom and c. the number of carbon atoms.   Screening of CO2RR catalysts of C3 products. a. The adsorption energy of CO, b. the reaction energy of C1-C1 coupling, c. the reaction energy of C1-C2 coupling on 100 screened alloy catalysts. All the energies are in eV.