Mechanism and Selectivity of Electrochemical Reduction of CO2 on Metalloporphyrin Catalysts from DFT Studies

Electrochemical reduction of CO2 to value-added chemicals has been hindered by poor product selectivity and competition from hydrogen evolution reactions. This study aims to unravel the origin of the product selectivity and competitive hydrogen evolution reaction on [MP]0 catalysts (M = Fe, Co, Rh and Ir; P is porphyrin ligand) by analyzing the mechanism of CO2 reduction and H2 formation based on the results of density functional theory calculations. Reduction of CO2 to CO and HCOO− proceeds via the formation of carboxylate adduct ([MP-COOH]0 and ([MP-COOH]−) and metal-hydride [MP-H]−, respectively. Competing proton reduction to gaseous hydrogen shares the [MP-H]− intermediate. Our results show that the pKa of [MP-H]0 can be used as an indicator of the CO or HCOO−/H2 preference. Furthermore, an ergoneutral pH has been determined and used to determine the minimum pH at which selective CO2 reduction to HCOO− becomes favorable over the H2 production. These analyses allow us to understand the product selectivity of CO2 reduction on [FeP]0, [CoP]0, [RhP]0 and [IrP]0; [FeP]0 and [CoP]0 are selective for CO whereas [RhP]0 and [IrP]0 are selective for HCOO− while suppressing H2 formation. These descriptors should be applicable to other catalysts in an aqueous medium.


Introduction
Carbon dioxide emissions from fossil fuel combustion have resulted in an accelerated increase in atmospheric CO 2 levels, and consequently, serious climate impacts and ocean acidification [1,2]. Conversion of CO 2 to value-added chemicals and fuels with sustainable and renewable energy sources helps to mitigate increasing CO 2 levels in the atmosphere [3][4][5]. Among the options for CO 2 conversion and utilization, electrochemical reduction driven by electricity generated from renewable energy sources such as sunlight and wind is a promising approach [6][7][8]. However, poor product selectivity and competing hydrogen evolution reactions (HER) due to proton reduction are among the major obstacles to the deployment of electrochemical reduction of CO 2 [9][10][11].
Tuning of the active site of a molecular catalyst and optimization of the operating pH have been used to control selectivity while suppressing HER. The intrinsic activity of a molecular catalyst determines the reduction potential of the electron transfer steps, pK a of protonated intermediates, and the binding energy of the adsorbate. This activity can be tuned by changing the metal ion or modifying the ligands. Computational electrocatalysis can provide a mechanistic understanding of the reaction chemistry, predict the reduction potential of a particular catalyst, and determine the pK a value and adsorption strength of intermediates. This information will in turn help screening of the metal centers and ligands and guide the selection of experimental conditions, such as applied potential and pH.
Molecular complexes, as homogenous catalysts, can reduce CO 2 to CO or HCOO − and proton to gaseous hydrogen (two electron-reduction products). CO can be fed to the Fischer-Tropsch process [12] for the synthesis of alkanes, whereas HCOO − can be

Results
Electrochemical reduction of CO 2 to CO and HCOO − , including the competing hydrogen evolution reaction, was studied using the reaction mechanism shown in Scheme 1. In Scheme 1, the catalytic cycle starts from the reduction of the catalyst ([MP] 0 ) to the reduced form ([MP] − ) at the corresponding reduction potential (E • ) (Step 1) [35][36][37] . Desorption of CO from [MP-CO] 0 regenerates the catalyst and completes the CO cycle (Step 15). Scheme 1. Schematic reaction steps of the electrocatalytic reduction of CO2 to CO (blue background) and HCOO − (yellow background). The yellow regions also contain proton reduction to hydrogen (red part). The grey region shows common intermediates shared by both cycles. M = Fe, Co, Rh and Ir; P = porphyrin ligand.  [40]. Moreover, the second reduction potentials are more negative than the corresponding first reduction potentials, indicating that the formation of [MP] 2− is less favorable and Step 4 is unlikely to happen under a typically applied potential of CO2 reduction (−1 V to −1.5 V) [46,47]. Consequently, [MP] 2− is less likely to have a significant contribution to either the CO or HCOO − /H2 cycle.  (Step 3 and 9), branching the reaction to the HCOO − /H2 and CO cycles, respectively. As shown in Table 1, the second reduction potentials of [MP] − (Step 4) are more negative than the corresponding first reduction potentials, indicating that protonation of [MP] 2− (Step 6) and CO2 binding on [MP] 2− (Step 9) are unlikely to contribute to the overall reaction in either of the cycles. Consequently, Steps 2 and 3 will be the dominant pathways. If protonation of [MP] − is favorable, the reaction will proceed toward the HCOO − /H2 cycle via [MP-H] 0 formations. Otherwise, CO2 will be Scheme 1. Schematic reaction steps of the electrocatalytic reduction of CO 2 to CO (blue background) and HCOO − (yellow background). The yellow regions also contain proton reduction to hydrogen (red part). The grey region shows common intermediates shared by both cycles. M = Fe, Co, Rh and Ir; P = porphyrin ligand.  Table 1 [40]. Moreover, the second reduction potentials are more negative than the corresponding first reduction potentials, indicating that the formation of [MP] 2− is less favorable and Step 4 is unlikely to happen under a typically applied potential of CO 2 reduction (−1 V to −1.5 V) [46,47]. Consequently, [MP] 2− is less likely to have a significant contribution to either the CO or HCOO − /H 2 cycle. Steps 2 and 3 will be the dominant pathways. If protonation of [MP] − is favorable, the reaction will proceed toward the HCOO − /H 2 cycle via [MP-H] 0 formations. Otherwise, CO 2 will be activated upon binding on [MP] − , steering the reaction toward the CO cycle. Therefore, the reaction free energy of protonation of [MP] − to [MP-H] 0 (Step 2) will dictate the selectivity between the CO and HCOO − /H 2 cycles.
The proton transfer step (Step 2) is represented by Equation (1). The dependence of the reaction free energy, ∆G • rxn (PT), on the pH and pK a of [MP-H] 0 can be derived with Equation (2). Equation (2) clearly shows that ∆G • rxn (PT) can be regulated by adjusting pH relative to the pK a of [MP-H] 0 . Consequently, the pK a of [MP-H] 0 can be used as a descriptor for HCOO − /H 2 selectivity. Since ∆G • rxn (PT) is positive for pH > pK a of [MP-H] 0 , proton binding to the metal center will not be spontaneous, making [MP] − available for CO 2 binding. For pH < pK a of [MP-H] 0 , ∆G • rxn (PT) will become negative, indicating the spontaneity of protonation and initiation of the HCOO − /H 2 cycle.  (3) and (4)) [14,39]. There is a constant difference of 34.88 kJ/mol between ∆G • rxn (HCOO − ) and ∆G • rxn (H 2 ) (H2O) , ∆G • rxn (H 2 ) (H2O) = ∆G • rxn (HCOO − ) + 34.88 kJ/mol [20], as presented graphically in Figure 1. This result indicates that HCOO − formation is always more thermodynamically favorable than H 2 formation. On  (Table 2), respectively, indicating that a pH below those values will make HCOO − /H2 thermodynamically favorable.  Figure 1.
The reaction free energy of [MP-H] − with a proton to produce H2 (Step 8, Scheme 1, Equation (5)) depends on the pH.
The pH at which ΔGrxn(H2) is zero is referred to as ergoneutral pH. A quantitative relationship between ΔG°rxn(HCOO − ) and the ergoneutral pH has been derived: pH(ergoneutral) = (41.84 kJ/mol − ΔG°rxn(HCOO − ))/5.70 kJ/mol (6) Equation (6) determines the pH at which ΔGrxn(H2) becomes zero. At a pH less than the ergoneutral pH, hydrogen evolution prevails. Equation (6) was used to construct the product distribution diagram shown in Figure 2. The vertical line at pH = 3.75 separates the HCOOH and HCOO − zone and the horizontal line at ΔG°rxn(HCOO − ) = 0 divides the zone of gaseous CO2 and HCOOH/HCOO − . The green diagonal line divides the zone between H2 and H + . The green region in Figure 2 corresponds to the selective reduction of CO2 to HCOO − where HER is suppressed. Similar schemes have been used to determine the optimum pKa value of the acids for non-aqueous solvents [48,49].    Figure 1.
The reaction free energy of [MP-H] − with a proton to produce H 2 (Step 8, Scheme 1, Equation (5)) depends on the pH.
The pH at which ∆G rxn (H 2 ) is zero is referred to as ergoneutral pH. A quantitative relationship between ∆G • rxn (HCOO − ) and the ergoneutral pH has been derived: pH(ergoneutral) = (41.84 kJ/mol − ∆G • rxn(HCOO − ))/5.70 kJ/mol (6) Equation (6) determines the pH at which ∆G rxn (H 2 ) becomes zero. At a pH less than the ergoneutral pH, hydrogen evolution prevails. Equation (6) was used to construct the product distribution diagram shown in Figure 2. The vertical line at pH = 3.75 separates the HCOOH and HCOO − zone and the horizontal line at ∆G • rxn (HCOO − ) = 0 divides the zone of gaseous CO 2 and HCOOH/HCOO − . The green diagonal line divides the zone between H 2 and H + . The green region in Figure 2 corresponds to the selective reduction of CO 2 to HCOO − where HER is suppressed. Similar schemes have been used to determine the optimum pK a value of the acids for non-aqueous solvents [48,49].
The the ergoneutral pH, hydrogen evolution prevails. Equation (6) was used to construct the product distribution diagram shown in Figure 2. The vertical line at pH = 3.75 separates the HCOOH and HCOO − zone and the horizontal line at ΔG°rxn(HCOO − ) = 0 divides the zone of gaseous CO2 and HCOOH/HCOO − . The green diagonal line divides the zone between H2 and H + . The green region in Figure 2 corresponds to the selective reduction of CO2 to HCOO − where HER is suppressed. Similar schemes have been used to determine the optimum pKa value of the acids for non-aqueous solvents [48,49].

CO Cycle
The pK a values of [MP-H] 0 determine the lowest pH at which protonation of the metal center of the complexes is not spontaneous, freeing the metal center for CO 2 14) is spontaneous with a reaction-free energy of −132.7 kJ/mol. Desorption of CO is also spontaneous (Step 15). These results indicate that CO formation from electrochemical CO 2 reduction at an applied potential of −1.28 V occurs via a series of electron and proton transfer steps at neutral pH or slightly alkaline conditions. formation of [FeP-COOH] can happen through a number of electron transfer steps in different orders combined with CO2 binding and proton transfer. The reduction of [FeP] − requires a significantly more negative potential (−1.75 V) and is unlikely to play a role in CO formation. Starting from [FeP] − , the free energy of binding of CO2 on [FeP] − is 22.19 kJ/mol. The pKa value of [FeP-COOH] 0 is 2.23, indicating that protonation of [FeP-COO] − becomes spontaneous only at a pH < 2.23. Further reduction of [FeP-COOH] 0 to [FeP-COOH] − happens at a potential of −0.40 V. As we showed above, a low pH of 2.23 will make the HCOO − /H2 pathway competitive and, therefore, CO formation through this reduction less likely for CO formation. Generally, electrochemical CO2 reduction happens in a slightly acidic or neutral pH. At this pH, the pathway following [FeP]  is also exergonic with a desorption free energy of −31.75 kJ/mol (Step 15). Therefore, CO2 reduction to CO can happen at −1.24 V in a slightly acidic or neutral pH on [FeP] 0 . This conclusion is consistent with the experimental results reported by Costentin et al. [37].

Discussion
The initial competitive binding of proton and CO 2 to [MP] − determines the pathway that the reaction will follow. Protonation will lead to the HCOO − /H 2 cycle. A more nucleophilic metal center tends to bind a small proton more preferably than CO 2 , making the reaction follow the HCOO − /H 2 cycle. On the other hand, a less nucleophilic metal center will bind CO 2 and steer the reaction toward the CO cycle. The charge distribution and electrostatic potential on the metal center in the complex provide an understanding of the preference for either proton or CO 2 binding [51]. The results of Bader charge analysis for [MP] 0 and [MP] − in Table 3 show that the electron density is significantly delocalized on the porphyrin ligand in [FeP] − and [CoP] − following the one-electron reduction. In contrast, the electron density is shared between the metal center and ligands in [RhP]  The orbital overlap between the metal center and its ligands provides a channel for the electrons to be shared and will affect the redox activity of the transition metal-porphyrin complex. Lack of such orbital overlap will hinder the electron channeling between the metal center and ligands and results in the increased nucleophilicity of the metal centers, steering the reaction towards the HCOO − /H 2 cycle. ing these metal centers prone to proton binding. In contrast, the negative electrostatic potential is delocalized among the metal centers and the porphyrin ligands in  The first reduction potentials of the metal complex also offer some insights into the preference for proton/CO2 binding. The first reduction potentials of [CoP] 0 and [FeP] 0 are −1.24 V and −1.28 V, respectively, whereas those of [RhP] 0 and [IrP] 0 are only −0.02 V and −0.26 V, respectively. It has been shown that CO2 binding is favored at a potential more negative than −1 V [35,52]. The more negative reduction potentials for [CoP] 0 and [FeP] 0 can be attributed to the more effective charge delocalization between the metal center and ligands. Further reduction of the activated CO2 is driven by the applied potential for the electron transfer steps. In contrast, the proton transfer steps are influenced by the chemical potential of the protons.
In summary, the electronic structure analyses shed additional light on the mechanism of electrochemical reduction of CO2 on transition metal-porphyrin complexes: CO2 reduction to CO can happen on [FeP] 0 and [CoP] 0 , whereas HCOO − is the main product on [RhP] 0 and [IrP] 0 , consistent with experimental observation [27,43,44]. Also consistent with the prediction from this study, hydrogen formation on [FeP] 0 in a strongly acidic aqueous medium has been reported [27,53]. We note that other factors could affect the stability of reaction intermediates and the chemical potentials of the reactive species, thereby changing the selectivity of the overall reaction. For example, Costentine et al. reported that H2 production The first reduction potentials of the metal complex also offer some insights into the preference for proton/CO 2 binding. The first reduction potentials of [CoP] 0 and [FeP] 0 are −1.24 V and −1.28 V, respectively, whereas those of [RhP] 0 and [IrP] 0 are only −0.02 V and −0.26 V, respectively. It has been shown that CO 2 binding is favored at a potential more negative than −1 V [35,52]. The more negative reduction potentials for [CoP] 0 and [FeP] 0 can be attributed to the more effective charge delocalization between the metal center and ligands. Further reduction of the activated CO 2 is driven by the applied potential for the electron transfer steps. In contrast, the proton transfer steps are influenced by the chemical potential of the protons.
In summary, the electronic structure analyses shed additional light on the mechanism of electrochemical reduction of CO 2 on transition metal-porphyrin complexes: CO 2 reduction to CO can happen on [FeP] 0 and [CoP] 0 , whereas HCOO − is the main product on [RhP] 0 and [IrP] 0 , consistent with experimental observation [27,43,44]. Also consistent with the prediction from this study, hydrogen formation on [FeP] 0 in a strongly acidic aqueous medium has been reported [27,53]. We note that other factors could affect the stability of reaction intermediates and the chemical potentials of the reactive species, thereby changing the selectivity of the overall reaction. For example, Costentine et al. reported that H 2 production and CO 2 reduction in a formic acid buffer on [FeP] 0 is roughly 50:50 [43]. Costentine et al. also showed that 90% faradaic efficiency for CO can be achieved by maintaining pH = 6.7 via NaOH addition. Based on Equation (2), ∆G • rxn (PT) becomes 24.23 kJ/mol at pH 6.7, preventing [FeP-H] 0 formation while allowing CO 2 to bind and react with [FeP] − and facilitating CO formation. Although the present study focuses on the thermodynamic aspect of the electrochemical reduction of CO 2 , we note that kinetics plays an important role in determining operating potential and product selectivity [54,55].

Computational Details
DFT calculations with B3LYP hybrid functional [56] and Grimme's D3 dispersion corrections [57] have been performed using Gaussian 16 [58]. The Stuttgart-Dresden effective core potentials [59] were adopted for transition metals, and a 6-31g(d,p) basis set was used for the main group elements. The choice of the basis set has been extensively tested against 6-311++G(d,p) in our previous study [20]. Our results showed that the computed reduction potentials using 6-311++G(d,p) and 6-31G(d,p) basis sets differed by less than 50 mV. The same level of theory was used for geometry optimization and frequency analysis. Optimized geometries were confirmed by the absence of any imaginary frequency. The SMD implicit solvation model for water was used [60]. The stability of wavefunctions for all species was also checked. Wavefunction optimization has been performed to include possible unrestricted open-shell diradicals. In this case, single-point energy of a specific spin state, including singlet (doublet), triplet (quartet), quintet (sextet), and septet (octet) was calculated first. These states were then compared, and the lowest-energy spin state was selected for wavefunction optimization. The optimized wavefunctions were used for subsequent geometry optimization and frequency calculations. This strategy allowed us to include possible open-shell radicals' structures. The spin multiplicities of the complexes and intermediates are presented in Table S1 in the Supplementary Materials.
The methods of determining the reduction potential of the electron transfer steps and pK a value of the carboxylate adducts ([MP-COOH] 0 and ([MP-COOH] − ) have been presented in our previous study [20]. The reduction potentials were reported by referencing the standard hydrogen electrode (SHE). The pK a value of [MP-H] 0 was calculated using the reaction free energy for [MP-H] 0 → [MP] − + H + . The chemical potential of H + , G • solv (H + ) = −1125.96 kJ/mol was determined using the solvated proton with six water molecules. The Multiwfn program [61] was used to perform Bader's atom-in-molecule (AIM) charge [62] analysis.

Conclusions
A mechanistic analysis of CO 2 reduction and H 2 formation on transition metalporphyrin complex catalysts has been performed based on the results of density functional theory calculations. Reduction of CO 2 to CO and HCOO − proceeds via the formation of carboxylate adduct ([MP-COOH] 0 and ([MP-COOH] − ) and metal-hydride [MP-H] − , respectively. The pK a of [MP-H] 0 determined the pH or buffer acid condition for selectively producing CO or HCOO − /H 2 . Furthermore, an ergoneutral pH determined from ∆G • rxn (HCOO − ) defines the lowest pH above which HER will be suppressed. Based on the pK a of [MP-H] 0 . These findings are in agreement with experimentally observed product selectivities. Furthermore, our study helps to rationalize the experimental observations and provides guidelines to design selective catalysts and choose optimal experimental parameters for CO 2 reduction.