Theoretical Study on the Mechanism of CO* Electrochemical Reduction on Cu(111) under Constant Potential

: In order to understand the mechanism of the electrochemical reduction of CO* on Cu(111), which competitively generates two intermediates, CHO* and COH*, we performed a ﬁrst principles calculation on these two electrocatalytic reactions, including the solvent effect and imposing a constant potential. The transition states of the two reactions under constant potential conditions were located by the electrochemical nudged elastic band (eNEB) method and the charge effect in the two reactions was elucidated by charge correction in the constant potential model and Bader charge analysis. It was found that the two reactions have different potential dependencies and involve different degrees of partial charge transfer. In addition, we conﬁrmed that the COH* reaction pathway exists only in a certain potential range when using the implicit solvent model.


Introduction
It is a promising green chemistry strategy to recycle carbon dioxide (CO 2 ) in the environment and convert it into useful fuels and chemicals [1]. Electrocatalyzed CO 2 reduction, as one of the numerous converting methods, has become a research hotspot in recent years due to its advantages such as mild reaction conditions, cheap reactants, and easy control of the reaction process. Among multitudinous electrocatalysts, copper has experimentally been reported as the only pure metal catalyst that can further reduce CO molecules into hydrocarbons [2][3][4], possibly due to its negative adsorption energy for *CO and positive adsorption energy for surface hydrogen H* [5]. However, various reaction paths in the process of CO* reduction lead to poor reaction selectivity, thus the elucidation of the mechanism of these reactions are of great importance for the designation of highly active and selective catalysts.
Plenty of theoretical studies on the mechanism of these reactions based on experimental observations have been reported up to now. In 2010, Peterson et al. performed density functional theory (DFT) calculations on the electroreduction of CO 2 on a Cu surface and proposed a detailed reaction path including a CHO* intermediate that leads to methane production [6]. After that, in 2013, Nie et al. proposed an H-shuttling COH* path and showed that the activation barrier for the formation of COH* would be reduced by the water-assisted H-shuttling process [7]. Later in 2020, Zhao et al. re-evaluated two reaction paths involving CHO* and COH*, respectively, using embedded complete active space second-order perturbation theory (emb-CASPT2) [8]. These works showed that the first step in CO reduction is critical to the selectivity of the whole reduction process.
So far, the reported theoretical works mostly used models that contain constant charge, which may lead to an error that the surface potential of catalysts changes with the reaction process. To make the computational model closer to the reality that electrochemical reactions proceed under constant potential and undergo charge transfer from or to the electrodes, some groups developed constant potential methods [9][10][11]. In this work, we apply a constant-potential local minimize and saddle point search method including both We first calculated the activation barrier and reaction energy of the reaction path under each potential. It is worth noting that there exists a shallow local minimum (denoted as intermediate state, MS) close to the TS. We performed another eNEB calculation between IS and MS and found a TS, which is the same as the TS between IS and FS. Figure 2a Previous studies believed that the process of CO* forming CHO* through the direct transfer of surface hydrogen is a chemical process rather than an electrochemical process [8] since no electrons transfer. That is to say, this reaction can only be written as CO * + H * −→ CHO * instead of CO * + H + + e − −→ CHO * . Based on this, the activation energy and reaction energy of the reaction should be independent of the potential. However, our result shows that the activation energy and reaction energy of this reaction does vary with the potential, albeit to a small magnitude. The activation energies decrease from 0.84 eV at U = 0.0 V to 0.57 eV at U = −1.6 V. This is because the system charge of the IS, TS, and FS are changed to keep the potential constant ( Figure 2c  Bader charge analysis of each atom (listed in Table 1) shows that charge redistribution occurs during this process from IS to TS, and from TS to FS. Overall, with the formation of CHO*, partial charge transfer occurred from the surface species back to the copper electrode. Specifically, when the surface hydrogen H* leaves the copper surface, it returns partial charges to the copper atom, and the carbon atom absorbs partial charges from the surface and then forms a C-H bond with hydrogen. By the time CHO* is finally formed, the charge of the surface species again partially returned to the copper surface. Comparing the charge changes of carbon atoms at various potentials, we found that they all need to obtain a charge of about 0.25 e in the process from IS to TS, while the charge difference of the system between the TS and the IS increases from close to 0 to about 0.25 e as the potential becomes negative. Combined with the change in the activation barrier with the potential, we speculate that one reason for the decrease in activation energy as the potential becomes negative is that it is easier for the charged surface to donate electrons, which makes the formation of the TS easier.

The COH* Path
The mechanism of *CO forming COH* intermediates through water-assisted H-shuttling was proposed by Nie et al. in 2013 [7], the subsequent steps are the formation of an additional O-H bond and then dissociating a water molecule to form a surface carbon atoms C* and then gradually hydrogenate to form methane. This reaction path can better explain the phenomenon of methane production in the CO 2 electroreduction experiment, and the participation of a water molecule also solves the problem of the high potential barrier for the formation of COH*.
When exploring the CO* to COH* mechanism using implicit solvent model, we found that the hydrogen of the final COH* intermediate could not point to the oxygen of the water molecule and be stabilized by hydrogen bonds, it would transfer back to the water molecule to form H 3 O + species. The energy of this structure is lower than the initial state energy proposed by Nie et al. in 2013 [7] and can be regarded as the initial state of a different path. This phenomenon is similar to the situation that Nie et al. encountered in 2014 [13], they stabilized COH* by pointing the hydrogen of COH* in the opposite direction. We also tested this opposite direction COH* structure and found that it can exist in implicit solvent model. This shows that the reason why the hydrogen of COH* transferred back to the water molecule is not the instability of the O-H bond of COH* in implicit solvent model, but the energy difference between O w − H · · · O C and O w · · · H − O C hydrogen bond structure ("· · · " represent a hydrogen bond, O w represent the oxygen atom of water, and O c represent the oxygen atom of COH*), once there is not enough energy barrier between the two structures, one will change into another structure spontaneously.
The relative energy of these two hydrogen bond structures will vary with the potential, both O w · · · H − O C and O w − H · · · O C structures have their existing potential conditions; in our model, a potential more negative than about −0.6 V is needed (−0.6 V is not an exact threshold) to prevent COH* hydrogen returns to the water molecule; similarly, the H 3 O + structure is stabilized under a potential more positive than −1.2 V. That is to say, the COH* reaction path we described here only exists in the range of about −1.2 V to −0.6 V. Figure 3 shows the optimized structures of IS, TS, and FS of the COH* reaction path under the potential of −0.8 V.  Figure 4a shows energies of structures along the MEP of CO* reduction to COH* from −0.8 V to −1.0 V. The activation barrier for this reaction path is quite low and goes lower when a more negative potential is applied. The result for −1.0 V are close to a non-barrier process, which means the IS of the reaction path will no longer exist and the reaction will proceed spontaneously once a H 3 O + species approaches the CO* at a potential lower than −1.0 V. On the other hand, when the potential is higher, the COH* intermediate will never come to the formation, as we mentioned before. The barrier between the IS and FS is quite interesting. At first, we thought the barrier was caused by the transition of the proton from the water molecule to the CO molecule; however, the vibration mode of the imaginary frequency of the TS is not the H stretching vibration, but the orientation twisting of the water molecule (see Supplementary Materials for details), which means the barrier is actually caused by the orientation change of the water molecule (and the orientation change of the water molecule is due to the charged surface), thus the COH* path has no truly TS. This result is consistent with the study of Zhao et al. in 2021 [14], as they calculated the COH* formation via Proton-Coupled Electron Transfer (PCET) pathway and find that no transition state observed, which means the bond break between the water molecule and the proton (or between the CO molecule and the proton, from the other side of the reaction) needs no activation energy in this reaction. This phenomenon may be due to the presence of hydrogen bonds.
The activation barrier and reaction energy obviously varies with the potential due to the significant system charge difference along the reaction. Figure 4b shows that the charge of adsorbed species changes almost synchronously with the system charge, and about 0.5 e charge was transferred during the reaction.
The Bader charge analysis (listed in Table 2) shows that the reduction of CO* to COH* needs about 0.28 e of charge, and the charge of the transferred hydrogen atom does not change much and remains around 0.7 e, as we can see in the Table 2. Such a hydrogen atom could be regarded as a proton. Then, we can describe the whole process clearly as the formation of the O-H bond on CO* via proton transfer, accompanied by charge transfer from the electrode to the adsorbed species, which can be written as CO * + H x+ + xe − −→ COH * .

CO Poisoning
It is necessary to investigate CO poisoning at different potentials. In order to describe the binding strength between CO and copper atoms, we calculated a crystal orbital Hamilton population (COHP). The projected density of state (pDOS) and projected COHP (pCOHP) of a Cu(111) slab with a CO adsorbate at the potential of 0.0 V are shown in Figure 5.  Figure S2). It can be seen from Figure 6a that the binding strength of CO decreases when the CO coverage goes higher, which shows that Cu(111) surfaces are not easy to be poisoned by CO. Additionally, different CO coverage structures have different responses to the variation of potential, high CO coverage structures have stronger CO binding strength as the potential becomes lower, but structures with low CO coverage have a smaller tendency. We know that during the process of CO adsorption on the copper surface, the copper atom donates electrons to σ − π bond between Cu and C.
As the negative potential is applied, the copper surface becomes more inclined to provide electrons; thus, the bonding will become stronger. One speculation as to the cause of the behavior of low CO coverage structures may be due to the shift of the d-band center of copper. As we can see in Figure 6b, the d-band center of the copper surface has been lowered when the potential increases, and a lower d-band center caused weaker Cu-C bonding and low CO coverage structure may be more sensitive to the effect of d-band center shift.

Constant Potential Model
Here we describe the constant potential model developed by Duan et al. [12] in brief. In a constant potential system, the number of electrons were changed during atomic structural changes, such a system can be treated as a grand canonical system. The corrected energy of the grand canonical system (H) that considers both atomic positions (r) and number of electron s(n) as variables can be written as: where E is the internal potential energy, r is the Cartesian coordinates, n 0 is the total number of nuclear charges, and φ is the effective energy level of the counter electrode. The relation between φ and the applied voltage U when taking the standard hydrogen electrode (SHE) as the reference electrode is: we determined the absolute potential of the SHE to be 4.44 V according to experiments [15,16]. The electron chemical potential then can be written as: where E f is the Fermi level of the electrode which we are concerned with, and can be obtained from DFT calculations. The full gradient of H(r, n) is then (∂H/∂r, ∂H/∂n) = (−f, µ ). Geometry optimization and saddle point search under constant potential can then be performed using this gradient combined with optimization algorithms. In an eNEB run, in addition to the given initial guess of the geometric structure, the initial guess of the number of electrons in the system is also needed. The number of electrons, as with the geometric structure, is also obtained through self-consistent calculations. More details and the program package can be found in Ref. [12].

Electronic Structure Methods
A periodic slab model was used to describe catalyst surfaces. We constructed a supercell containing four layers of 3 × 3 Cu atoms with the bottom 2 layers fixed to simulate the bulk properties. A vacuum space of 15 Å thickness along the z-axis was set to avoid interactions between surface adsorbed species and the lower surface of another periodic structure above.
All electronic structure calculations in this paper were performed in the Vienna Abinitio Simulation Package (VASP) [17,18] within the framework of DFT. The plane wave pseudopotential basis set was used in calculations, and the cutoff energy of the plane wave basis set was set to 400 eV (The results of the convergence test of the cutoff energy are provided in the Supplementary Materials, Table S1). A 3 × 3 × 1 Gamma-centered k-point mesh was used to sample the Brillouin region. Electronic exchange and correlation were described using the Perdew-Burke-Ernzerhof (PBE) functional [19], and the threshold for electronic energy in SCF iteration was set to 1 × 10 −6 eV. Since the system we are concerned with is an open shell structure, spin polarization is considered. We describe the weak interactions such as hydrogen bonds by the Grimme DFT-D3 dispersion correction method [20]. The solvent environment of the system is described by an implicit solvent model implemented in the VASPsol program [21][22][23] by Hennig's group. In the implicit solvent model, the relative permittivity was set to 80 to simulate the aqueous environment. Geometry optimization calculations proceeded under constant potential, transition states were located using a climbing image-nudged elastic band (CI-NEB) method combined with the constant potential algorithm, as implemented in the eNEB code [12]. The local minimum and transition state structures were fully relaxed under constant potentials until the maximum atomic force was less than 0.05 eV/Å. All transition states were confirmed by frequency calculations, and the imaginary frequencies and vibration modes are provided in Supplementary Materials. The MEP curves are plotted using energies vs. reaction coordinates, each reaction coordinate corresponds to the relative position of each atom in a given reaction system, thus the specific form of reaction coordinates depends on the motion trajectory of atoms along the MEP. The energies we use here were corrected by adding −(n − n 0 )V sol to the corrected energy in Equation (1) to remove the interaction between the charged system and its compensating charge, which is actually nonexistent. Bader charge analysis was finished by the program package developed by Henkelman's group [24][25][26][27], and COHP calculation was done using Lobster program [28,29].

Conclusions
We performed a detailed mechanistic study of the two competing reactions for the Cu(111) catalyzed electroreduction, the formation of *CHO intermediate and the formation of *COH intermediate, via a constant potential model, taking charge effects and solvent effects into account.
We found that the elementary step of CO* forming CHO* intermediate through surface hydrogen H* transfer, which was previously considered to be a non-electrochemical and potential-independent process, actually involved partial charge transfer and have potential dependence on the activation barrier and reaction energy.
We explored the elementary step of CO* reduction to COH* intermediate, confirming that this path exists only in a specific potential interval. This process involves significant charge transfer, thus the activation energy barrier and the reaction energy are greatly affected by potential changes.
Comparing the reaction paths of these two intermediates, we infer that the first step of reduction of CO* favors the formation of CHO* at high potentials and COH* at lower potentials. The instability of COH* intermediates at high potentials may be able to explain the relatively negative production potential of methane in Cu catalyzed CO 2 electroreduction experiments [2].
We investigated the CO poisoning on a Cu surface. It is found that the binding strength between CO and copper atoms decreases when the CO coverage goes higher. CO binds stronger on the Cu surface when the potential goes lower because of more electron donation to the Cu-C bond, but this trend is smaller for low CO coverage structures, possibly due to greater sensitivity to the effect of d-band center shift.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/catal13060960/s1, Table S1: The convergence test of the cutoff energy of the plane wave basis set; Table S2: The imaginary frequencies of the TSs; Figure S1: The vibration mode of the imaginary frequencies; Figure