Theoretical Insights into Heterogeneous (Photo)electrochemical CO2

Review Cite This: Chem. Rev. XXXX, XXX, XXX−XXX

pubs.acs.org/CR

Theoretical Insights into Heterogeneous (Photo)electrochemical CO2 Reduction Shenzhen Xu† and Emily A. Carter*,‡ Department of Mechanical and Aerospace Engineering and ‡School of Engineering and Applied Science, Princeton University, Princeton, New Jersey 08544-5263, United States

Downloaded via TULANE UNIV on December 19, 2018 at 00:09:44 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

ABSTRACT: Electrochemical and photoelectrochemical CO2 reduction technologies offer the promise of zero-carbon-emission renewable fuels needed for heavy-duty transportation. However, the inert nature of the CO2 molecule poses a fundamental challenge that must be overcome before efficient (photo)electrochemical CO2 reduction at scale will be achieved. Optimal catalysts exhibit enduring stability, fast kinetics, high selectivity, and low manufacturing cost. Identifying catalytic mechanisms of CO2 reduction in (photo)electrochemical systems could accelerate design of efficient catalysts. In recent decades, numerous theoretical studies have contributed to our understanding of CO2 reduction pathways and identifying rate-limiting steps. Although a significant body of work exists regarding homogeneous electrocatalysis for CO2 reduction, this review focuses specifically on the theory of heterogeneous (photo)electrochemical reduction. We first give an overview of the relevant thermodynamics and semiconductor physics. We then introduce important, widely used theoretical techniques and modeling approaches to catalysis. Recent progress in elucidating mechanisms of heterogeneous (photo)electrochemical CO2 reduction is discussed through the lens of two experimental systems: pyridine (Py)-catalyzed CO2 (photo)electrochemical reduction at p-GaP photoelectrodes and electrochemical CO2 reduction at Cu electrodes. We close by proposing strategies and principles for the future design of (photo)electrochemical catalysts to improve the selectivity and reaction kinetics of CO2 reduction. 2.1.3. Debate about PyH+ as Merely a Proton Donor 2.1.4. Proposed Dihydropyridine (DHPaq vs DHP*) Intermediate 2.1.5. Adsorbed 2-Pyridinide (2-PyH−*) as a Possible Intermediate 2.2. CO2 Reduction on Metal-Oxide and MetalChalcogenide Photoelectrodes 2.2.1. Titanium Oxide (TiO2) 2.2.2. Indium Oxide (In2O3) 2.2.3. Co3O4 Hexagonal Platelets 2.2.4. Cadmium Telluride (CdTe) 2.2.5. Copper Indium Sulfide (CuInS2) 2.3. Non-Metal-Based Photoelectrodes 3. Heterogeneous CO2 Electrochemical Reduction 3.1. CO2 Reduction Mechanisms on Cu Surfaces 3.1.1. Mechanisms Predicted by the CHE Model 3.1.2. Linear-Scaling Relationship and Overpotential Volcano Plot 3.1.3. Mechanisms Predicted by Water-Solvated and H-Shuttling Model

CONTENTS 1. Introduction 1.1. The Big Picture: Why CO2 Reduction 1.2. Thermodynamics of (Photo)electrochemical CO2 Reduction and Relevant Semiconductor Physics 1.3. Theoretical Techniques for Examining Heterogeneous CO2 Reduction on Electrode Surfaces 1.3.1. Density Functional Theory (DFT) Calculations and Beyond 1.3.2. Solvation Models 1.3.3. Computational Hydrogen Electrode (CHE) 1.3.4. Water Solvated and H-Shuttling Model 1.3.5. Potential-Dependent Activation Energy 1.3.6. Ab Initio Molecular Dynamics (AIMD)Based Metadynamics and Constrained Molecular Dynamics (CMD) Simulations 1.3.7. Experimental Observables for Validating Theoretical Models 2. Heterogeneous CO2 Photoelectrochemical Reduction 2.1. CO2 Reduction on p-GaP Photoelectrodes with Pyridine (Py) as a Cocatalyst 2.1.1. Proposed Catalytic Mechanisms Arising from Experiment and theory 2.1.2. Heterogeneous vs Homogeneous © XXXX American Chemical Society

B B

B

E E H I I J

K K

N N N O O Q Q R R R S S S S U

L L

Special Issue: Computational Design of Catalysts from Molecules to Materials

L N

Received: July 29, 2018

A

DOI: 10.1021/acs.chemrev.8b00481 Chem. Rev. XXXX, XXX, XXX−XXX

Chemical Reviews 3.1.4. Mechanisms Predicted by AIMD-Based Metadynamics and CMD Simulations 3.1.5. Effects of Solvated Cations and Anions on CO2 Electrochemical Reduction 3.2. CO2 Reduction Mechanisms on Other Metal (non-Cu) Electrode Surfaces 3.3. CO2 Reduction Mechanisms on Metal Alloy Electrode Surfaces 3.4. Mechanistic Understanding of C2 Product Formation 4. Discussion and Outlook 4.1. Possible Intermediates for Py-Catalyzed CO2 Photoelectrochemical Reduction on Semiconductor Surfaces 4.2. Theoretical Insights for Cocatalyst Design for Heterogeneous CO2 Photoelectrochemical Reduction 4.3. Theoretical Insights for the Design of CO2 Electrochemical Reduction Catalytic Electrodes 4.4. Limitations of the Current Electrochemical Models and Outlook 5. Summary Author Information Corresponding Author ORCID Notes Biographies Acknowledgments References

Review

progress in the field of CO2 electrochemical reduction. Photoelectrochemical CO2 reduction as “artificial photosynthesis” has been studied since the late 1970s.20−32 Photoelectrochemical reduction systems fall into two main categories based on the system’s configuration: photoelectrochemical cells (PECs) and photocatalytic particles. Readers can refer to the reviews from the Bocarsly group19 and the Peng group33 for detailed discussions of the experimental progress in the field of CO2 photoelectrochemical reduction. Besides (photo)electrochemical CO2 reduction, several other CO2 heterogeneous conversion processes deserve mention, e.g., thermocatalytic CO2 reduction to methanol, CO2 conversion to CO by reverse water−gas shift (RWGS) and by solar thermochemical reduction. For the thermocatalytic CO2 reduction to methanol at gas−solid interfaces, the reaction is CO2 + 3H2 → CH3OH + H2O with a standard formation enthalpy of −49.3 kJ/mol at 298.15 K. This CO2 conversion process proceeds on a Cu/ZnO/Al2O3 catalyst at moderate pressure (50−100 bar) and temperature (473−573 K) conditions on an industrial scale.34,35 References 36 and 37 review theoretical insights combined with experimental evidence into mechanisms of thermocatalytic CO2 reduction. For the RWGS process, the reaction is CO2 + H2 ⇔ CO + H2O with a standard formation enthalpy of 42.1 kJ/mol at 298.15 K. This endothermic reaction requires a high temperature to drive the equilibrium conversion to CO. References 37 and 38 review the main type of catalyst, i.e., metals (Cu, Pt, Rh, Au, etc.) supported on oxides (Al2O3, ZnO, ZrO2, SiO2, CeO2, etc.), and their corresponding catalytic mechanisms. For solar thermochemical reduction, the main idea is to utilize oxygen vacancies, generated in metal oxides (e.g., CeO2) at high temperatures by concentrated solar fluxes, to extract oxygen from H2O to make H2 and reduce CO2 to CO by filling the oxygen vacancies in the solid material at a relatively low temperature. References 36 and 37 also review progress in solar thermochemical reduction. A detailed discussion of thermocatalytic CO2 reduction, RWGS, and solar thermochemical reduction is beyond our scope; here we focus only on heterogeneous (photo)electrochemical CO2 reduction.

W X Y Z AA AB

AB

AB

AC AD AE AF AF AF AF AF AF AF

1. INTRODUCTION 1.1. The Big Picture: Why CO2 Reduction

The world’s total energy consumption is projected to increase from 30% to 50% in the next 20 years.1,2 Burning fossil fuels, the traditional source of energy for humans, causes excessive emissions of greenhouse gases (e.g., CO2) and accelerates climate change.3 Moreover, burning fossil fuels is unsustainable due to finite resources. Using captured CO2 as a feedstock to generate organic fuels could help mitigate the above problems by reducing fossil fuel extraction and use. Rather than simply spilling the electrical energy from excess renewable sources (e.g., solar mid-day and wind at night in the middle of the U.S.), it could be stored as chemicals and fuels, by converting renewable energy to electronic energy and then on to chemical energy. Two promising future technologies, corresponding to the above path of energy conversion, are electrochemical and photoelectrochemical CO2 reduction, which could produce, e.g., alcohols or hydrocarbons as valuable fuels or fuel precursors from CO2 recycling. At the same time, a common side reaction, the hydrogen evolution reaction (HER), must be suppressed. Although H2 gas is also a fuel, it has the drawback of low energy density and inherent infrastructure problems relevant to H2 transport, storage, and safety concerns.4 That is why the goal here is to inhibit production of H2 while maximizing production of energy-dense carbon-based fuels. Electrochemical CO2 reduction on various metal electrodes with different cell designs has been studied extensively by Hori, Bard, Lvov, Fujishima, Sakata, etc., over the past ∼30 years.5−18 We refer the readers to Hori’s work9,10 and the review19 from the Bocarsly group for a comprehensive review of the experimental

1.2. Thermodynamics of (Photo)electrochemical CO2 Reduction and Relevant Semiconductor Physics

Thermodynamically, reactants can be driven to form products only if the system’s free energy decreases after the reaction. For example, in the case of CO2 being reduced to hydrocarbons, xCO2 + y(H+ + e−) → CxHy‑2nO2x‑n + nH2O, the electrons must have a chemical potential (i.e., free energy) higher than the free energy required to drive the reaction, i.e., G(reactants) ≥ G(products) where G means the Gibbs free energy. In electrochemical reduction, the energy of electrons can be tuned by an applied potential. The potential at which G(reactants) = G(products) is defined as the reduction potential. Reduction potentials for CO2 being reduced to CO, HCOOH, CH3OH, CH4, and CH2O, and H+ being reduced to H2 are shown in Table 1. In photoelectrochemical systems, the energy of photoexcited electrons is determined generally by the conduction band edge (CBE) (introduced in the next paragraph) of the semiconductor photoelectrode. Thermodynamically, CO2 photoreduction is favorable if the energy of the CBE is more negative than that of the reduction potential. For CO2 photoelectrochemical reduction, the electrode materials must be comprised of semiconductors. Here, we briefly introduce some basic concepts about them to understand B


Chemical Reviews

Review

Table 1. Formal Electrochemical Reduction Potential (E0) Relative to the Standard Hydrogen Electrode (SHE) at pH = 7 for Relevant Reactions of CO2 Reductiona Reaction −

•−

CO2 + e → CO2

E0 vs SHE (V) −1.85

CO2(g) + 2H+ + 2e− → HCOOH(aq) +

−0.665

−

CO2(g) + 2H + 2e → CO(g) + H 2O

−0.521

−

CO2(g) + 4H + 4e → CH2O(aq) + H2O CO2(g) + 6H+ + 6e− → CH3OH(aq) + H2O

−0.485 −0.399

CO2(g) + 8H+ + 8e− → CH4(g) + 2H 2O

−0.246

2H+ + 2e− → H 2(g)

−0.414

+

semiconductor acts as the photocathode with a EF lower than that of the electrolyte solution. Once the p-type photocathode is immersed in the electrolyte, electrons will flow from the electrolyte to the semiconductor, building up an electric field in the semiconductor near the interface. This electric field causes a band bending of the semiconductor’s electronic structure. In the case of a p-type semiconductor−electrolyte interface, the semiconductor’s CB and VB will bend downward (direction: semiconductor → electrolyte, as shown in Scheme 1). This band Scheme 1. Schematic of Band Bending at the p-Type Semiconductor/Electrolyte Interfacea

a

Data taken from refs 39, 40.

why. Most materials we categorize as metals, semiconductors, or insulators, depending on their distribution of electronic energy levels. Because of the sheer number of atoms in a macroscopic crystal, the electronic wave functions on each atom combine to form states characterized by continuous bands of energy levels. The energy levels of interest are a set of lower-energy bands comprised of filled electronic states called the valence band (VB), and the VB maximum (VBM)the highest energy filled electronic state in the darkis the VB edge. A set of higherenergy bands of unoccupied electronic states (in the dark) is called the conduction band (CB), and the lowest energy CB state is called the CB minimum (CBM) or the CB edge. For metals, the CB and VB overlap, forming a continuum of electronic energy states. The electrical conductivity of metals is high because electrons have access to delocalized, empty electronic states accessible with negligible thermal activation energy. For semiconductors, the two bands are separated by a small or modest band gap (Eg), which is a quantum mechanically forbidden energy zone. A substantial energy is required to excite electrons from the VB to the CB in a semiconductor, leading to much lower electrical conductivity compared to metals. However, even in the dark, the conductivity of a semiconductor is nonzero because of the Boltzmann thermal tail of electronic energies that allows for a nonzero occupation of the CB. By contrast, for insulators, Eg is so large that the VB is completely filled and the CB is completely empty under typical experimental conditions (e.g., at room temperature), making it impossible to redistribute electrons with an applied electric field. Thus, insulators exhibit an infinitesimal electrical conductivity. Semiconductors generally achieve higher conductivity instead via deliberate photoexcitation of electrons from the VB to CB. Indeed, it is their abilitybecause of their band gapsto absorb light in an energy range that overlaps well the solar spectrum that enables their use as photovoltaics, photoelectrodes, and photocatalysts.41 Next, we discuss the variation of the semiconductor band structure at the semiconductor−electrolyte interface. We first need to define the concept of the Fermi level (EF). The most relevant definition for our purposes is that it defines the chemical potential of electrons in a material, a thermodynamic quantity that we can use to analyze electron transfer at material interfaces. If the EF of a semiconductor electrode placed in an electrolyte solution is not equal to the EF of the electrolyte, then electron transfer will occur at the semiconductor−electrolyte interface following the direction from high EF to low EF. Let us use a ptype semiconductor as an example (recall that “p-type” refers to a material doped to be electron-deficient). In a PEC, a p-type

a

EF is the Fermi level of the semiconductor. ER is the redox potential, which can be considered as the Fermi level of the electrolyte.

bending facilitates electron transport to the photocathode− electrolyte interface, explained as follows: when a light source illuminates a semiconductor, some fraction of electrons in the VB will be excited to the CB by the incident photons, leaving holes (positive charges due to the absence of electrons at those energies) in the VB. The photocathode’s CB/VB bending near the interface applies a force on electrons/holes in opposite directions, driving electrons/holes toward/away from the interface. Thus, band bending enhances electron−hole pair separation, which is critical to photoconversion efficiency that relies on irreversible separation of charge carriers (vide infra). The electrons drifting to the photocathode−electrolyte interface are available for reduction reactions, provided the energy of the CBM at the interface lies above the reaction free energy required (otherwise violating thermodynamics!). A similar mechanism for holes occurs at the photoanode−electrolyte interface. The main difference is that a photoanode is an n-type semiconductor (doped to have excess electrons), leading to an opposite CB/ VB-bending direction compared to the photocathode. This opposite bending (upward toward the interface) facilitates hole transport to the photoanode-electrolyte interface, allowing oxidation of electrolyte molecules (e.g., H2O). The above description of a p-type photocathode and an n-type photoanode corresponds to the two possible photoelectrodes in a PEC; one can have both types simultaneously in a tandem PEC or one type only with a metal counter-electrode in a standard PEC. An alternative configuration of a photocatalytic system is the suspension of semiconductor particles in an electrolyte C


Chemical Reviews

Review

ions are preferred elements because they have no exchangeenergy driving force to trap electrons,41 consistent with experimental observations that d0 and d10 transition-metal oxides (TMOs) exhibit promising photocatalytic capability.43 Besides exchange-energy driving forces, other descriptors may help design materials to avoid charge traps. One example is the Carter group’s earlier work studying electron transport in pure and doped hematite using electrostatically embedded clusters and the small polaron model.44 Here, the motivation for doping hematite was to increase the carrier concentration. However, we need to choose doping elements carefully to avoid introducing trap states. Because the electron carriers in hematite are localized Fe 3d electrons, the electron transport mechanism involves small polaron hopping, in which the electron motion is coupled to lattice distortions. Four n-type substitutional dopant elements were studied (Ti, Zr, Si, and Ge), and it was found that only Ti4+ strongly traps electrons. Ti3+ has a much larger ionization potential (IP) than Zr3+; thus, Ti4+ has a stronger driving force to trap an electron compared to Zr4+. For Si and Ge dopants, they form polar covalent bonds to oxygen in contrast to the ionic bonding between Ti/Zr and oxygen. The photoexcited electron must occupy a Si/Ge−O antibonding state, which is energetically unfavorable. Si and Ge dopants therefore do not trap electrons. The findings in this work44 suggest that doping the system with ions having a low IP or forming more covalent bonds to oxygen may increase the charge-carrier concentration without inducing trapping. Following the discussion of electron−hole pair lifetime and charge transport in semiconductors, we briefly discuss band gap optimization. For the photocathode used in the CO2 reduction, we want to maximally utilize the solar spectrum and have a sufficiently high CBM to drive the desired reduction reaction. The optimal band gap for photoconductivity is ∼1.5 eV45 but that energy generally does not provide sufficient driving force for the reaction, with most estimates suggesting band gaps around 2 eV being optimal for photocatalysis.43 Let us first consider conventional photocatalysis, namely TMOs. For many, the VBM consists mostly of O 2p character, which is +3 V or even more positive vs the normal hydrogen electrode (NHE). Considering the band gap has to be less than 3 eV to capture a significant fraction of the solar spectrum, the corresponding CBM is likely to be lower than the reduction potentials of CO2 conversion to many fuels (e.g., HCOOH, CO, CH2O, and CH3OH) and hence not viable. One possible solution is to substitute elements that shift up the VBM of the semiconductor photocatalyst. For TMOs, some anion dopants (C, N, S, and P) have lower electronegativity than O. The p-orbitals of these dopants are higher in energy compared to the O 2p state. Therefore, the VBM of anion-doped TMOs might be shifted up.41 Use of metal chalcogenides, i.e., full substitution of all oxide ions, is of course the limiting case that indeed is explored further in this review (vide infra). Another strategy for tuning the band gap and band edge positions is alloying. For example, a 1:1 solid solution of MnO:ZnO was predicted to be a promising material with a band gap of 2.6 eV and favorable CBM/VBM positions.46 Solid solutions of Fe2+ oxides with other M2+ oxides (Mg, Mn, Ni, and Zn) also showed promise.47 In terms of a driving force for CO2 reduction, all of the Fe1−x(Mg/Mn/Ni/ Zn)xO, solid solutions are predicted to have a sufficiently high CBM to reduce CO2. Another descriptor was discovered related to band gap center differences: alloys exhibit favorably reduced band gaps (compared to relatively large band gaps of pure endmembers) if the band gap centers of the pure materials are

solution, where both the reduction and oxidation reactions occur on the same particle. The photocatalytic particle can be considered as a combination of the photocathode and photoanode. A schematic explaining the mechanism of photocatalysts for CO2 reduction is shown in Figure 1.

Figure 1. Simplified scheme of the mechanism of photocatalytic CO2 reduction over a semiconducting photocatalyst.

To achieve efficient photon-induced reduction reactions at a photocathode surface, it is essential to have sufficient lifetime of the photoexcited electrons. Electron−hole pair recombination is the limiting factor for the lifetime of photoexcited electrons. We need materials in which the photoexcited electrons can reach the photocathode−electrolyte interface before recombining with holes. One possible strategy is to exploit the spin-forbidden nature of triplet to singlet transitions (inspired by a design strategy for organic light emitting diodes, which are essentially inverse photovoltaics42). The idea41 is to convert the initially formed singlet state of the electron−hole pair to a nearly degenerate but slightly lower-energy triplet state (the latter to provide driving force but the former for maximum coupling) so that subsequent electron−hole recombination will be significantly suppressed due to the forbidden transition. An only slightly lower-energy triplet state avoids a substantial energy loss while allowing the electron−hole pair to live long enough to eventually separate rather than recombining. In compound semiconductors, such as a metal oxide or chalcogenide, ligand to metal charge transfer (LMCT) transitions induced by incident photons might be optimal excitation types because the ligand hole and the metal electron are located on different atoms, leading to a small exchange splitting between the singlet and triplet states. The exchange splitting also might be tunable via doping, as long as such dopants do not produce midgap states that could act as charge carrier traps. In general, if we could devise a material with ligand states (e.g., O 2p states) dominant at the VBM and metal states dominant at the CBM, then photoexcited electrons and holes will reside on different atoms. The electron−hole pair then should have a sufficient lifetime to separate and diffuse to the photoelectrode surfaces. Beyond the lifetime of the electron−hole pair, fast transport of photoexcited electrons and holes in opposite directions without encountering trap sites is essential to photoconversion efficiency. In the case of transition metal (TM) dopants, the large intra-atomic exchange energies especially associated with first-row TMs could drive or inhibit charge trapping, depending on their particular electronic structure.41 For example, if a photoexcited electron encounters a high-spin d4 Mn3+ ion, that electron very likely will get trapped to form a high-spin d5 Mn2+ ion, because of the net four d−d exchange terms gained (the origin of the stability of half-filled shells). Based on this same principle, one can readily show that d0, d10, and high-spin d5 TM D


Chemical Reviews

Review

failure of the electronic structure prediction (e.g., significant underestimations of charge transfer barriers and Eg). Further, the lack of a derivative discontinuity in pure XC density functionals adds to the artifacts observed.62 One of the more famous examples is predicting Mott insulators to be metals, such as the TMOs FeO and NiO.64,65 A charge-transfer example is the failure of DFT to predict any barrier at all for dissociative adsorption of dioxygen on Al (111), which was shown to involve abrupt charge transfer from higher-level theory, and from which a significant barrier consistent with experiment naturally emerges.66 Although not relevant for the materials discussed here, for completeness we also mention that the mean-field nature of DFT inaccurately describes strongly correlated systems, in which dynamical electron correlation is critical to obtaining even a qualitatively correct description.67,68 Another limitation of conventional DFT with LDA and GGA XC functionals is the lack of nonlocal correlation such as induced dipole−induced dipole interactions, i.e., the dispersion or van der Waals (vdW) interaction. The latter is important for describing weak binding such as solvent-electrode interactions and adsorbate−adsorbate lateral interactions. Thus, a dispersion correction (e.g., the Grimme D269 and D370 corrections) is typically required, as pointed out below in section 2.2.1. Unfortunately, a major class of materials used as photoelectrodes and photocatalysts are TMO semiconductors, which contain localized d-electrons on the transition metal cations that suffer severe SIE. We therefore need approaches beyond conventional DFT, which can correct the SIE and overdelocalization of these strongly correlated electrons. The two most common theories are DFT+U and hybrid DFT. DFT+U was originally proposed by Anisimov et al.71 The key idea is to specifically treat the XC interactions of the strongly correlated delectrons within a HF-like theory, while treating the rest of the system with conventional DFT. The intra-atomic HF treatment is employed to remove the SIE of the localized d-electrons. The computational cost of DFT+U is not significantly higher than that of standard DFT, which makes it a practical and attractive method for (photo)electrochemical systems involving TMOs. The main weakness of the DFT+U method is the need to choose a U value, actually U-J as proposed by Dudarev et al.72 Here, U/J corresponds to the intra-atomic Coulomb/exchange interaction. Frequently, U is chosen empirically, by comparing with experimental observables of known materials; in this case such DFT+U calculations must be considered semiempirical rather than fully ab initio. Several first-principles schemes have been proposed to determine U, such as constrained DFT calculations to determine U,73,74 a self-consistent evaluation of U at the level of DFT+U,75 and an approach proposed by the Carter group that determines U-J based on ab initio unrestricted HF (UHF) calculations on electrostatically embedded clusters.76,77 Unlike other schemes for deriving U that still depend on approximate XC and its attendant error, the UHF approach provides values uniquely defined based on exact exchange calculations. The other approach to SIE correction is hybrid DFT. The key idea is to include a fraction of the exact HF exchange into the XC functional applied to all electrons in the system. Some commonly used hybrid functionals are the B3LYP78 functional, which contains 20% HF exchange and is popular for studying molecules; the PBE079,80 functional, which incorporates 25% HF exchange into the Perdew−Burke−Ernzerhof (PBE) exchange;61 and the Heyd−Scuseria−Ernzerhof (HSE) functional,81 which contains 25% screened (short-range) HF exchange with the remaining XC functional described by PBE.

significantly different. Moreover, the CB and VB edges can have different elemental character in the alloy, which additionally could reduce electron−hole recombination as the photoexcited electrons and holes will be located on different atoms. Finally, a direct band gap vs an indirect band gap is also a tunable factor in band gap engineering, which equivalently corresponds to a balance of high light absorption vs slow electron−hole recombination. Because this review does not focus on the semiconductor physics of photoelectrochemical systems, we refer readers to the reviews from the Carter group,41 the Bocarsly group,19 and the Lewis group48−51 for a comprehensive review of basic semiconductor concepts and engineering approaches. 1.3. Theoretical Techniques for Examining Heterogeneous CO2 Reduction on Electrode Surfaces

1.3.1. Density Functional Theory (DFT) Calculations and Beyond. Approximations to formally exact DFT are the most widely used techniques for the quantum simulation of chemistry and materials science.52−54 They provide a reasonable balance between computational efficiency and accuracy. For (photo)electrochemistry, DFT can be a powerful tool for exploring electronic, thermodynamic, and kinetic properties of such phenomena. We describe briefly some ab initio methods that build upon standard DFT, as well. The last part of this subsection then discusses modeling approaches implemented with DFT for simulation of material surfaces and their attendant properties that can be obtained. Two fundamental theorems underlying DFT were proposed and proved by Hohenberg and Kohn (HK) in 1964.55 They first proved a one-to-one mapping between the ground-state density and the external potential (HK theorem I). The second theorem (HK theorem II) states that the density minimizing the total energy (a functional of electron density) is the system’s exact ground-state density. Although DFT is an exact theory, it does not tell us anything about the specific expression of the energy density functionals. Thus, we have to approximate the kinetic and exchange-correlation (XC) functionals to put DFT into practice. The standard Kohn−Sham (KS) approach56 to DFT maps an interacting, many-body problem onto a noninteracting, single-particle problem, in which a set of noninteracting orbitals (the KS orbitals) are solved for self-consistently, and the total density is calculated by summing the occupied KS-orbital densities. The kinetic energy is approximated by applying the single-particle kinetic energy operators on the noninteracting KS orbitals. As a consequence, the entire accuracy of KS-DFT hinges on the formulation of the XC functional, which has to account for the difference in kinetic energy between interacting and noninteracting electrons as well as capturing accurately the quantum effects of exchange and correlation. Two of the first and still commonly used forms of the XC functional are the local density approximation (LDA)57 and the generalized gradient approximation (GGA).58−61 With the computational power available today, KS-DFT calculations (using LDA and GGA) can be applied to surface electrochemical systems containing several hundred atoms within a reasonable computing time. The physical limitations of the LDA and GGA XC functionals originate from the so-called self-interaction error (SIE)62 such that their approximate exchange energies do not exactly cancel the spurious Coulomb repulsion of each electron with itself, as happens naturally in Hartree−Fock (HF) theory.63 This problem in pure DFT approximations often produces inappropriate electron overdelocalization and a qualitative E


Chemical Reviews

Review

the quasiparticle wave functions and eigenvalues should not be thought of as single-particle quantities. Instead, they reflect properties of the correlated many-electron system. The quasiparticle eigenvalues from the GWA give the IP and EA values and thus can be compared directly to PES and IPES measurements. If, however, we want to compare to measured optical band gaps that correspond not to IP-EA but to a neutral single excitation, then the electron−hole interaction must be taken into account. The Bethe-Salpeter equation (BSE)93 subsequent to a GWA calculation can be applied to calculate the material’s optical spectrum and band gap. Note that the importance of the electron−hole interaction varies strongly with the dielectric properties of the material as well as the dimensionality of the system. In terms of computational cost, both GWA and BSE are much more expensive than conventional DFT methods. Here, we briefly summarize the level of accuracy and appropriateness of the various first-principles techniques mentioned above for various materials properties. Standard DFT methods with LDA or GGA as the XC functional are the most computationally efficient. They generally are reliable for structural optimization; the relaxed structures can be used as input for higher accuracy, more costly calculations of other properties. For metallic systems and non-TMOs, standard DFT also are generally trustworthy for calculating electronic structure (e.g., densities of states (DOS), atomic charges, etc.) and various thermodynamic properties (energy, equation of states, etc.). For materials with TM ions (e.g.., strongly correlated electron materials), standard DFT methods tend to fail. We therefore need the DFT+U method at least, or even the more expensive hybrid XC functional method. DFT+U is only slightly more expensive than standard DFT, but the nonunique choice of U is a drawback. In principle, the U value depends on the TM and its oxidation state, as well as the surrounding chemical environment. Thus, DFT+U is inappropriate for phenomena involving TM ions undergoing multiple spin/valence-state transitions,94,95 where multiple empirical U values might have to be chosen, which makes the predictions less convincing. The alternative strategy to reduce SIE in standard DFT is to use hybrid XC functionals, mixing overly ionic HF theory with overly metallic DFT, which often yields to a fortuitous cancellation of errors in predictions of band gaps, charge transfer barriers, and other electronic and thermodynamic properties. However, in its periodic implementation, hybrid DFT is much more expensive than conventional DFT and DFT +U, limiting the system size that can be examined. Moving on from ground-state properties, direct comparisons with measurements of excited-state properties require the much more expensive TDDFT, GWA, or BSE methods. GWA yields the IP-EA quasiparticle band gap while TDDFT and BSE furnish absorption spectra and optical band gaps. The quasiparticle gap is a very good approximation to the optical gap whenever the exciton binding energy is small (as it is in most inorganic solids). However, these excited-state GWA and TDDFT methods are not appropriate for calculations of thermodynamic properties or structural optimization as they do not deliver total energies from which such properties can be derived. Construction of appropriate physical models for simulation of a (photo)electrochemical environment also impacts the level of theory that can and should be employed. Because heterogeneous CO2 reduction reactions typically occur at photoelectrode surfaces, we next introduce two common approaches for the physical modeling of semiconductor surfaces.

The main differences between DFT+U and the hybrid functional method are (1) Instead of selecting U and J values, one selects the fraction of HF exchange, frequently adjusted to best reproduce measured properties, although adiabatic connection theory justifies the PBE0 choice.79 (2) The hybrid functional method applies the same XC functional to all electrons, essentially a compromise between an overly ionic HF treatment and the overdelocalized DFT solution. By contrast, DFT+U assumes the source of SIE to originate entirely from the localized d-/f-electrons. This oversimplification is problematic for mixed ionic/covalent materials (e.g., Cu2O82). (3) The hybrid functional method requires a full calculation of nonlocal HF exchange, which is computationally expensive, especially so for periodic systems. Thus, application of the hybrid functional method is more limited than the DFT+U method. Although DFT is a ground-state theory, the CBM/VBM eigenvalues often are used to estimate the materials’ band gap. However, the eigenvalues difference obtained from DFT does not have a direct physical correspondence to a band gap measured in experiments. Even though DFT+U and hybrid DFT account for the SIE and correct the underestimation of band gaps to some extent, the eigenvalue gaps based on ground-state DFT methods are not rigorously comparable to measured band gaps. To get the right band gap for the right reason, we therefore need methods that truly account for the physics of excited electronic states. The time-dependent DFT (TDDFT) method proposed by Runge and Gross83 is able to probe a system’s excited states. Excitation energies can be evaluated by the positions of poles in the density response function. TDDFT also can be used to calculate and directly compare to absorption spectra measured experimentally, although the accuracy depends on the XC functionals and typically exhibits an error of ∼0.5 eV.84 The primary applications of TDDFT are for molecules and nanostructures. Some recent examples of TDDFT applied in periodic surfaces/solids (relevant to (photo)electrochemical systems) include the absorption spectrum of dye-sensitized TiO2,85 photogenerated hole dynamics,86 and methoxy photo-oxidation dynamics87 at rutile TiO2 (110) surfaces (based on real-time TDDFT). Besides TDDFT, a many-body perturbation theory based on the Green function method provides a scheme to calculate the materials’ IP (measured by photoemission spectroscopy (PES)) and electron affinity (EA, measured by inverse photoemission spectroscopy (IPES)). In the GW approximation (GWA),88 the inserted electron/hole in the system together with the surrounding Coulomb hole/electron induced by the XC effect is treated as a single entity (the so-called quasiparticle). Due to the screening effect of the Coulomb hole/electron (corresponding to the inserted electron/hole), quasiparticles interact via a weak and screened interaction W instead of the strong Coulomb interaction. Hedin’s equations89 for quasiparticles are solved by an approximation of the self-energy operator in terms of G0 (the single-particle Green function of the noninteracting system, obtained from, e.g., DFT-based methods such as KS-DFT, DFT +U, hybrid DFT, etc.) and W. As the DFT-based wave functions and eigenvalues can serve as initial inputs for the GWA, a judicious choice of the KS-DFT setup is important for obtaining meaningful results by GWA. Several versions of GWA with different levels of self-consistency provide users the freedom to balance accuracy and computational efficiency, e.g., nonselfconsistent GW (G0W0), updating only eigenvalues in G (GW0), or in both G and W (GW),90 updating wave functions as well as eigenvalues (quasiparticle self-consistent GW).91,92 Note that F


Chemical Reviews

Review

It is also a challenge to model electric fields across the solid− liquid interface that we encounter in electrochemistry. Nørskov and co-workers thus proposed an approach to simulate the interfacial electric field without using a charged system.109−111 In their model, water molecules are added on top of a metal surface, while hydrogen atoms are further placed on top of the water molecules. The hydrogen atoms ionize to produce hydronium ions, and electrons are transferred to the metal surface. Thus, an electric field builds up across the metal−water interface. The capacitance C of this modeled, electrochemical double layer determines the integrated free energy stored in the double layer as a function of the electrode potential, which is an important physical parameter to compare with experiments. In their test calculations on Pt(111), they found that C = 22.7 μF/ cm2 for one water layer, C = 23.6 μF/cm2 for two water layers, and C = 24.5 μF/cm2 for three water layers, all of which are consistent with the measured C = 20 μF/cm2 for Pt(111).112 This study indicates that even a single water layer is sufficient to describe the electrochemical double layer at the metal/water interface.110 The other approach to modeling electrode surfaces is the embedded cluster model. The crystal surface can be studied by a cluster consisting of a limited number of atoms located in the periodic surface’s lattice sites. Boundary effects have to be taken into account, typically by some form of embedding that mimics the extended system. In the case of covalent solids with localized bonds, capping atoms are used to saturate the dangling bonds at the boundary, with care taken to preserve the correct local geometry.97,113,114 In the case of metals, more complicated embedding techniques generally are required.115,116 In the case of ionic materials, point-charge embedding methods are sufficient to represent the long-range electrostatic interactions in the infinite crystal.116−118 It is also important to test the convergence of relevant properties with respect to cluster size and validate the cluster model using the results of the periodic slab model as the benchmark.113 Advantages of the cluster model vs the periodic surface slab model101 are: (1) a reliable, uncomplicated simulation of systems with net charge, (2) hybrid DFT and post-HF correlated wave function methods are readily employed at a fraction of the cost compared to the periodic case, and (3) continuum solvation models are well validated and computationally efficient. Our own tests show that inclusion of implicit solvent in cluster model calculations (in the quantum chemistry package ORCA version 3.0.3) increases the computational cost by a factor of 1.2 (vs the vacuum case). By contrast, inclusion of implicit solvent in the periodic slab model (in VASP version 5.4.4) increases the computational cost by a factor of 2.3 (vs the vacuum case). To close this subsection, we outline some important surface physical/chemical properties (relevant to heterogeneous (photo)electrochemical reactions) obtainable from electronic structure calculations using one of these two physical models:

The first approach is the periodic slab model. DFT calculations with periodic boundary conditions (PBCs) properly account for the periodicity of crystals. The general approach is to build up a slab with multiple atomic layers in a periodic cell, where the surface slab typically is cleaved from a structurally relaxed bulk material that defines the surface lattice vectors of the slab. A vacuum region is introduced in the surface normal direction, thick enough to avoid interactions between the periodic images of the slab. Spurious long-range, periodically repeated interactions of the dipole moments of the surface and/ or adsorbates have to be eliminated. One solution is to symmetrize the two surfaces (top and bottom) of the slab including the adsorbates, where the bottom surface is a mirror image of the top surface.96 This approach is relatively expensive computationally, not only because we have to double the number of adsorbates in the supercell, but also because we need to incorporate more layers in the slab to isolate the effect of the mirrored adsorbates from each other. An alternative in the case of covalent solids is to study phenomena on one surface while saturating the dangling bonds with (pseudo)hydrogen atoms on the reverse side.97 To cancel any spurious net dipole−dipole contribution to the system’s total energy and local potential, selfconsistent corrections can be added, e.g., by turning on IDIPOL and LDIPOL in the software package VASP.98−100 In the case of metals, the artificial interaction between periodic images of the adsorbates’ dipole moments is usually minimal due to screening by the metal slab. The main limitations of the periodic slab model101 are: (1) for systems with net charge as might be studied in electrochemical phenomena, an unphysical uniform compensating charge background is required to avoid divergence of the Ewald summation and (2) hybrid DFT and post-HF methods are exceedingly expensive. To elaborate on the first limitation, a charged periodic slab is problematic unless large supercells are used to damp out the error. The origin of this error lies in the unphysical electrostatic interaction of the charge with its periodic images and the background. The error converges slowly with respect to the supercell lattice constant L. The leading term scales with q2/L, which can be calculated from the Madelung energy of an array of point charges q with a neutralizing background.102 Makov and Payne103 later proved that, for isolated ions, the quadrupole moment of the charge distribution further contributes a term scaling as Q/L3, where Q is the quadrupole moment. These scaling terms (L−1 + L−3) require large supercells for convergence. This empirical extrapolation is also questionable,104 as demonstrated by comparing the L−1 and (L−1 + L−3) extrapolation schemes for the formation energy of the +2 vacancy in diamond (plotted vs L−1). Other approaches of truncating or compensating the longrange tail of the bare Coulomb interaction during computation of the electrostatic potential had been proposed to remove the unwanted interaction.105−107 However, these methods neglect the polarization outside the supercell. A relatively new and robust approach was proposed by Van de Walle and co-workers in 2009,104 in which ab initio finite-size corrections for this unphysical electrostatic interaction could be calculated with fast convergence with respect to the supercell size L. This scheme works well for isotropic three-dimensional solids.104,108 For a periodic slab in contact with an electrolyte, however, this approach cannot ameliorate the unphysical behavior of a neutralizing background charge penetrating into the electrolyte environment (vide infra). Large supercells therefore are still required to damp out the error induced by PBCs.

(1) Electrode surface reconstructions and electrode−electrolyte interface structure, with several examples relevant to the heterogeneous CO2 reduction.114,119−123 (2) Adsorption energies, used to determine if a chemical species can adsorb spontaneously on a material’s surface, with several examples relevant to heterogeneous CO2 (photo)electrochemical reduction.113,114,124−128 (3) Band edge positions. One approach is to reference them to redox potentials, e.g., the hydrogen electrode. Van de Walle and Neugebauer129 proposed a semiempirical G


Chemical Reviews

Review

The first is the “explicit solvation” model that directly includes solvent molecules in the electronic structure calculation. The interactions between solvent, solutes, adsorbates, and electrode surfaces can be described accurately by this model. However, an explicit solvation model is computationally expensive. A specific atomic configuration of the electrode/solvent system lacks the thermodynamic averaging over the locations of atoms comprising the solvent environment, as well. Although one can sample over the configurational phase space using various approaches to approximate the thermodynamic ensemble of the system, e.g., molecular dynamics (MD), the Monte Carlo method, or even creating multiple configurations manually, such sampling makes the computation much more expensive and usually impractical. A quantum mechanical/molecular mechanical (QM/MM) partitioning (first introduced in the 1976 paper of Warshel and Levitt140) can help reduce the computational cost while sampling. The key idea is to treat the reacting portions of the system with relatively expensive QM, while the rest of the environment is treated with much cheaper MM classical force fields. Please refer to Field et al.141 for a detailed introduction of QM/MM potential development and implementation. For heterogeneous electrochemical reactions, the electrode surface and the adsorbates should be treated quantum mechanically while the solvent region could be modeled with explicit solvent molecules governed by a classical force field. An example can be found in a theoretical study of CO2 adsorption on rutile TiO2 (110)142 where QM/MM simulations were carried out to preequilibrate a water/TiO2 (110) interface. Another example used the QM/MM method to calculate the homogeneous reduction potential of aqueous TM complexes,143 where the TM complexes were treated quantum mechanically and the surrounding solvent environment was modeled by the classical MM method. Despite QM/MM’s success in reducing computational cost, one major issue lies in how to treat the boundary between the QM and MM regions, especially when the boundary cuts through a covalent bond, leading potentially to an unphysical charge state of the QM atom.144 The second is the implicit solvation model that represents the solvent as a polarizable medium characterized by a dielectric constant ε. The charge distribution in the solvent is represented as an electric field that gets polarized by and responds to the presence of the solute. Detailed discussion on this topic can be found in various review articles.145,146 Two commonly used implicit solvation models are the COnductor-like Solution MOdel (COSMO)147 and the Solvation Model based on solute electron Density (SMD).148 The COSMO model treats the continuum medium as a conductor, in which the calculation of the electrostatic stabilizing interaction arising from the solute polarizing the continuum medium is greatly simplified. A variant of the COSMO model is called the Conductor-like Polarizable Continuum Model (CPCM). There are multiple options for the cavity definitions used by CPCM, e.g., cavities defined by Bondi vdW radii,149 United Atom HF (UAHF) radii, and S(Simplified) UAHF radii.150 The SMD model is based on the generalized Born approximation, which approximates the Poisson and Poisson−Boltzmann equations so that they can be solved analytically. As a relevant example, the SMD, CPCMvdW, and CPCM-SUAHF solvation models were compared for treating solvation effects on acidity constants of substituted pyridinium (PyH+) ions and pyridinyl (PyH•) radicals.136 The drawback is that these continuum medium models cannot describe specific interactions (e.g., hydrogen bonds) between the solute and solvent molecules.

approach to align band edges of different materials by calculating interstitial hydrogen levels in the bulk. The hydrogen level is defined by the Fermi-level position, where the formation energy of the positive (interstitial H+) charged state equals to the formation energy of the negative (interstitial H−) charged state. They showed that this hydrogen level is universal (∼−4.5 eV vs vacuum) and almost independent of the host material. The absolute energy scale of a certain material’s band edge then is assigned by aligning its hydrogen level to that of Si and setting the CBM of Si to the measured EA of Si. A second approach later proposed by Toroker et al. calculates absolute band edge positions by using the band gap center (BGC) to reference the vacuum level.130 Perdew and Levy proved that the BGC position predicted by DFT is formally exact.131 Its absolute position is referenced to the vacuum level and is equal to the work function for an intrinsic semiconductor (i.e., no doping), calculated from the electrostatic potential for a free electron in vacuum and the BGC inside the slab. The absolute positions of the VBM and CBM then are determined by subtracting or adding one-half of the band gap to the BGC position. The band gap is obtained from a separate, higher-level calculation on the bulk material (e.g., hybrid DFT for an estimate or the GWA for a rigorous value). The major differences between these two methods is that the earlier one is semiempirical while the later one is not, and secondly the later method calculates the band edges of a specific material surface (and could therefore include the effects of adsorbates and explicit electrolyte), while the earlier method determines the band edge positions within the bulk material with no dependence on the interface structure or environment. (4) Homogeneous/heterogeneous reduction potentials132 of the chemical species formed by electrochemical or photoelectrochemical reduction reactions, which indicates either the least-required applied potential on an electrode or the least negative CBM of a photoelectrode for producing a targeted chemical species. Calculated standard reduction potentials also can help exclude possible intermediates by comparing their reduction potentials with either experimental applied potentials or CBM positions of the photoelectrodes. (5) Acidities/hydricities of chemical species in solution or adsorbed on surfaces.125,133−139 These two quantities can help us decide whether a certain molecule (or H*; “*” refers to a surface adsorption site and X* refers to adsorbed X) can easily donate a proton/hydride in a proton-coupled electron transfer (PCET)/hydride transfer (HT) step in the CO2 reduction process. The acidity of a molecule in solution also tells us if the molecule is more stable in a protonated or deprotonated state under an experimental (pH) condition. (6) Pathways of surface reactions and the corresponding activation and reaction free energies, which are essential to understanding the mechanisms of the reactions and for estimating reaction rates. 1.3.2. Solvation Models. Since (photo)electrodes are in direct contact with electrolyte solution, employment of an appropriate solvation model is critical to properly simulating heterogeneous CO2 reduction at (photo)electrode surfaces. In this section, we introduce several solvation models that treat solvent effects on surface reactions. H


Chemical Reviews

Review

experimental system under study), then the corresponding electron energy level is defined to be the RHE potential. The RHE scale therefore differs from the SHE scale by the proton free energy difference between pH = 0 and pH equal to the value in the system under study. The convenience of using the RHE scale is that the pH variance of the electrolyte solution is incorporated into the definition. If a certain potential U (vs RHE) is applied to the electrode, then the free energy of a (H+ + e−) pair can be represented by 1/2 μ(H2) − eU. Thus, a complete thermodynamic energy landscape of a series of PCET steps with respect to the applied potential (vs RHE) can be generated by the CHE model. Finally, the SHE reference defined above allows one to explicitly incorporate a pH dependence. If a certain potential U (vs SHE) is applied to the electrode at a certain pH, then the free energy of a (H+ + e−) pair can be represented by 1/2 μ(H2) − eU − 2.3 × kBT × pH. One of the most important applications of the CHE model is the prediction of a reaction’s minimum overpotential, consisting of multiple PCET steps from a purely thermodynamic perspective. For example, the minimum overpotential for CO2 electrochemical reduction to CH4 can be obtained by taking the difference between the potential required to make the complete thermodynamic energy landscape downhill and the formal electrochemical reduction potential of CO2 to CH4 (0.17 V vs RHE).159 The advantage of this CHE model is that it can generate the free energy landscape (without explicit kinetic barriers) of a reaction consisting of multiple PCET steps and then identify the potential-determining step (PDS). Here we first define the concept of the limiting potential (UL) of a PCET step. Consider the generic PCET reaction, A* + H+ + e− → AH*, where A is a CO2 reduction intermediate. Then the free energy of the (H+ + e−) pair is conveniently equal to (1/ 2)μ(H2(g)), where μ(H2(g)) is the free energy of H2 gas at standard state. If the free energy change of the generic PCET reaction is ΔG, then the UL (vs RHE) is simply −ΔG/e. Now consider a reduction reaction consisting of multiple PCET steps, such as the example shown in Scheme 2: A* + 3H+ + 3e− → AH* + 2H+ + 2e− → AH2* + H+ + e− → AH3. The PCET step with the most negative UL (vs RHE) is the PDS. The thermodynamic overpotential then can be calculated from |UL − U0|, where U0 is the standard reduction potential. Scheme 2 illustrates how the PDS is determined. The central idea of the CHE model, namely that of formulating the energies of H+ and e− as a pair, based on the equilibrium of H+ + e− ⇔ 1/2H2, conveniently enables calculation of the free energy changes of the PCET steps at a certain applied potential. The main drawback of the CHE model is that it neglects the activation energies of all PCET steps and therefore is not reliable for accurate kinetics. The step with the most positive free energy difference is considered to be the ratelimiting step based on the assumption that the PCET steps occurring at electrode/solvent interfaces have small kinetic barriers, surmountable at room temperature.159,160 For example, barriers for O2 reduction to OOH on Pt161 and OH reduction to H2O on Pt160 were calculated to be in the range of 0.15−0.25 eV at the potential needed to make the elementary step exergonic. However, this assumption might be valid only for limited cases, so its generality is still in question. 1.3.4. Water Solvated and H-Shuttling Model. The water solvated and H-shuttling model considers detailed kinetics,162,163 in contrast to the CHE model. It is a waterassisting, hydrogen-shuttling model used for calculating the activation energy of a hydrogenation process, A* + H* →AH*.

By combining the implicit and explicit solvation models, the mixed implicit−explicit solvation model emerges, which balances accuracy and computational efficiency. Mixed solvation models treat the solvent region close to the solute/electrode with atomic-level representation of solvent molecules, while the longer-range remaining solvent polarization is treated by a continuum model, allowing, e.g., hydrogen bonding to solutes to be accounted for explicitly while still accounting for longer range electrostatic response. For instance, in Lessio and co-workers investigations of CO2 photoelectroreduction on solvated GaP surfaces,120,124,133,151 a monolayer of half-dissociated water moleculesthe most favorable configuration for water found by ab initio thermodynamicswas explicitly treated at the GaP surface, while the remaining region was treated by the SMD implicit solvation model. Another approach to handling solvent effects at electrode surfaces is joint DFT (JDFT), which is an ab initio description of an electronic system in equilibrium with a liquid environment, e.g., solid−liquid interfaces. The free energy functional of the system is divided into three contributions: the electronic subsystem (solute), liquid subsystem (solvent), and coupling interactions between the two subsystems.152−154 Each of these contributions can be described with different levels of approximations. For example, KS-DFT can be applied to the electronic system, while a classical density functional, which is developed on the basis of either an implicit solvent model (e.g., the linearized Poisson−Boltzmann model)154 or an explicit solvent model155,156 (taking into account the nonlinear dielectric response of the solution and the solvent structure), can be applied to the solvent environment.154,157 Advantages of JDFT are as follows. (1) It can treat properly charged systems with PBCs, electrochemical double layer structures, electrode surfaces under applied potential, etc.154 (2) Users have more freedom to balance accuracy and computational efficiency in terms of choosing functional levels for the electronic system, solvent environment, and coupling term. (3) The physics of thermodynamic averaging of the solvent environment in equilibrium with the electronic system is incorporated into JDFT. One major drawback (an open area of JDFT research) is the less-developed functionals for the solvent environment and coupling term.154 With the above discussions of the approaches for building up electrode surfaces and the surrounding solvent environment in hand, we now introduce in the next subsections 1.3.3−1.3.6) models for calculating free energy changes and activation energies of elementary electrochemical reaction steps, needed to gain insights into reaction mechanisms. 1.3.3. Computational Hydrogen Electrode (CHE). The CHE model generates the thermodynamic energy landscape of a series of PCET steps under an applied potential (referenced to the reversible hydrogen electrode (RHE) potential).158 The CHE model provides an elegant way of avoiding the explicit treatment of solvated protons and electrons in a PCET reaction. The free energy of a (H+ + e−) pair in a reaction A* + H+ + e− → AH* at zero applied potential (i.e., at the RHE potential) is defined by the energy of a (H+ + e−) pair in equilibrium with gaseous H2 at 1 atm, i.e., H+ + e− ⇔ 1/2H2.158,159 Note that the absolute energy level of the electron in the above reaction depends on the pH of the electrolyte solution. If the pH is zero, then the absolute energy level of the electron (making the reaction in equilibrium) is defined as the standard hydrogen electrode (SHE) potential. If the pH is an arbitrary value (i.e., typically the pH value of the electrolyte solution in an I


Chemical Reviews

Review

electrode to adsorbates, due to the limited surface area in a constant-charge simulation (a finite-size error). A capacitor-like model was developed by Chan and Nørskov to obtain constantpotential activation energies based on work function calculations and Bader charge analysis.166,167 The activation energy of a surface reaction with charge Δq transferring from the electrode to adsorbates at a constant work function (or potential) Φ can be obtained from ΔE(Φ) = ΔE(Φref) − Δq(Φ − Φref),167 where Φref is the electrode work function of the initial state before the charge transfer occurs and ΔE is the activation energy; ΔE(Φref) can be calculated by an approach shown in ref 166. Examples applied to CO2 (photo)electrochemical reduction are CO2 to CO reduction kinetics on Cu(100),123,168,169 and formation kinetics of 2-PyH−* on GaP(110).170 In summary, the advantage of this correction scheme is that it takes into account the finite-size error caused by the unphysical change of the electrode work function during the simulation of electron transfer from the electrode to adsorbates. The required inputs of this approach are straightforward, which are the work function change and the amount of transferred charge across the electrode surface. The calculation of this capacitor-like correction is based on a constant-charge simulation and does not involve any charged cell calculations. One potential drawback of this method is its fundamental assumption that the electrostatic energy change caused by electron transfer is purely capacitive. This assumption works well for simple PCET steps where no large dipoles are associated with the adsorbates and the solvent-reorganization energy is negligible.166 However, if the PCET steps involve adsorbed molecules with a redox center far from the surface or with a very large dipole, then the capacitor-like assumption may be invalid. Another technical issue of this correction technique is the estimation of the amount of charge transferred across the electrode surface. The corresponding accuracy may be sensitive to the choice of an algorithm for attributing charge density to electrode surfaces. For example, recent embedded correlated wave function calculations171 verified that because of the overlap of tails of the electron density between the electrode and the solvent, the actual charge state of a proton in the electrochemical double layer is not unity. Fractionally charged protons will affect the energetics and kinetics of, e.g., a Volmer reaction (H+ + e− + * → H*). Even though this capacitor-like model accounts for the finitesize error caused by the constant charge simulation, it does not physically simulate a constant-potential system that should be represented by a grand canonical ensemble where the amount of charge in the electrode can vary. Explicit charging schemes can control the potential of the electrode but require an approach to deal with a net charge in the periodic slab supercell. Proposed strategies include: (1) Applying a homogeneous background charge to compensate for the extra charge in the electrode slab and referencing the electrode’s absolute potential via the socalled double-reference method.172 (2) Constructing a charged plane of continuum opposite charge in the supercell cell, performing as a counter electrode.173 (3) Introducing a conductor-like continuum medium (with an infinite dielectric constant, performing as a metal electrode) in the vacuum region, isolated from the surface slab. The counter charge can be represented by the continuum image charge density in the region of the conductor-like medium.174 (4) Introducing implicit solvent approaches where the distribution of continuum counter charge is determined by the Debye screening length of the electrochemical double layer,154,175 as employed in the

Scheme 2. Schematic Illustrating How to Determine the Potential-Limiting Step and Limiting Potentiala

“A” is a representative chemical species being reduced. Here, the reaction AH + H+ + e− → AH2 is the potential-limiting step because it has the largest free energy increase (+1.5 eV) at U = 0 V (vs RHE). At U = −1.5 V (vs. RHE), the energy landscape becomes favorable for all steps, and −1.5 V defines the limiting potential. The length of each arrow represents -U (vs. RHE) eV. The number of arrows is equal to the number of electrons in each step. For example, the energy of the step A + 3H+ + 3e− increases by 3 × (−(−1.5)) = 4.5 eV when U changes from 0 to −1.5 V (vs. RHE), because there are three electrons in this step. Three arrows represent this energy change, as shown in the plot. a

H* is used to represent a proton−electron couple under an “equilibrium” potential U0, where H+ + e− ⇔ H*. The calculated activation energy is labeled as Eact0(U0). At an arbitrary potential U, Eact(U) = Eact0(U0) + β(U−U0) based on Butler−Volmer theory,164 where β is set to be 0.5. The water-solvated and Hshuttling model consists of two mechanisms, a water-solvated mechanism (direct H* transfer to A*) and an H-shuttling mechanism (water molecule shuttling H* from the surface to A*). The advantage of this model is that the kinetic barriers of the PCET steps are calculated explicitly. Moreover, the dependence of the barriers on the applied potential is considered. However, the generality of its major assumption that the kinetics of a PCET step can be represented by an adsorbed hydride transferis still in question. Also, inclusion of only one or two explicit H2O molecules for the H-shuttling mechanism (assisted by H2O) might be insufficient for certain systems, as well. For example, our recent work165 on the protonation kinetics of 2-pyridinide (2-PyH−*) on GaP(110) included a half-dissociated monolayer of explicit H2O molecules adsorbed on the surface to properly represent the semiconductor−water interface, and four explicit H2O molecules to solvate a proton in solution (an Eigen cation). Three explicit H2O molecules were involved in the PT process;165 thus, inclusion of only one or two explicit H2O molecules would have been insufficient. Details of the water-solvated and H-shuttling model can be found in the original paper163 and a review.36 1.3.5. Potential-Dependent Activation Energy. In this subsection, we discuss the issue of constant electrode potential simulation vs constant electrode charge simulation. In electrochemical experiments, PCET steps occur under the condition of a constant electrode potential. However, most computational studies are done with a a fixed number of electrons, i.e., constant charge, leading to the problem of electrode work function (potential) changes during a charge-transfer process from the J


Chemical Reviews

Review

reproduces measured electrode potentials. For example, a recent theoretical study185 of CO2 electrochemical reduction used experimental PZC values for Cu, Ag, Au, and Ni electrodes as benchmarks to validate use of the Poisson−Boltzmann continuum solvation model; in this case, the average error was small, only 0.09 V. Such PZC comparisons have been used in other electrochemical model benchmarking studies.154,157,186 A third metric is the electrostatic potential profile from the electrode surface into the bulk electrolyte. Although it is very difficult for experiments to directly measure this potential profile, qualitative comparisons with the well-established Gouy−Chapman−Stern model187−189 are useful for validating electrochemical simulations. For example, a recently developed constant-potential electrochemical model yielded a simulated electrostatic potential profile that reproduced well the constantfield region of the internal Stern layer and the exponential-decay region of potential contributed by the Gouy−Chapman diffuse layer.176 Beyond the above metrics for benchmarking electrochemical models of the surface double-layer, another validation scheme is to predict the current density vs voltage (C−V) curves via microkinetic models to make a direct comparison with a macroscopic measurement. (A note of caution: As there are multiple inputs to such models, agreement of theory and experiment does not necessarily guarantee that the model is correct, as multiple solutions could give rise to the same C−V curve. Therefore, one can only conclude that the model is consistent with the data.) Microkinetic models are constructed as follows. First, a complete set of catalytic reaction steps is identified. Then, the rate constant for intermediate A going to intermediate B is calculated typically from transition state theory (kA→B = (kBT/h) × exp(−ΔG⧧/ kBT), where kA→B is the rate constant of the reaction step A→ B and ΔG⧧ is the free energy barrier). The reaction rate of this step is rA→B = kA→B × θA, where θA is the surface coverage/concentration/partial pressure of intermediate species A. Then, the reaction rates of all the catalytic steps are coupled and one solves a steady-state condition, in order to predict the total reaction rate J for the final product formation. For CO2 electrochemical reduction, the free energy barrier ΔG⧧ of each catalytic step depends on the applied potential U. The predicted total reaction rate J therefore is a function of U. This J(U) can be compared with the experimental C−V curve to verify as a function of U the ratelimiting step and the corresponding activation energies predicted, e.g., by first-principles-based electrochemical models. For instance, simulations of C−C bond formation on Cu(100) predicted that at low overpotentials, CO* dimerization is the central pathway whereas at high overpotentials, the reaction of CO* with CHO* dominates. Microkinetics simulations assuming these pathways matched the measured C−V curve very well.185 Another study simulated CO2 electrochemical reduction to CH4 on Cu(100), Cu(111), and Cu(211);190 microkinetic modeling predicted that Cu(211) is the most active facet for CO2 reduction. Their predicted slope of the log(J)−V curve matched experiment well, indicating a probably correct identification of the rate-limiting step and the corresponding barrier as a function of applied potential. Beyond the nonuniqueness of its solutions mentioned above, another limitation of the microkinetic model is that it only represents the intrinsic rate of a multistep reaction, leaving out other factors such as those related to spatial variations in concentration. For example, local pH variations at electrode surfaces compared to the bulk pH and local CO2 concentration

JDFT method. (5) Building up a broad region of uniformly distributed continuum counter charge in the vacuum region above electrode surface, defined as a jellium slab of counter charge. This jellium slab is immersed in an implicit solvent environment that occupies a fraction of the vacuum region above the electrode surface, the rest is just vacuum connecting to the back side of the surface slab. This strategy is called the solvated jellium method.176 It is beyond the scope of this review to dig further into the details and the pros/cons of each chargecompensating strategy. Please refer to recently published work176 for a more comprehensive introduction/analysis of the above schemes and benchmarking calculations for comparing the solvated jellium method with the previously proposed homogeneous-background-charge technique172 and the counter electrode method.173 1.3.6. Ab Initio Molecular Dynamics (AIMD)-Based Metadynamics and Constrained Molecular Dynamics (CMD) Simulations. Metadynamics is an accelerated sampling method used to explore the system’s free energy landscape (as a function of collective variables (CVs)).177 A process can be accelerated by adding positive Gaussian potentials (bias potentials) to the system’s real energy landscape during the MD simulation. A multidimensional (multiple CVs) free-energy landscape can be obtained from the summation of all added bias potentials if the simulation time is sufficiently long enough. The drawback here is that one has to choose CVs of interest that might bias the outcomes and insights gleaned. The choice of Gaussian potential functions may affect the results as well, although the influence is expected to be minor. Examples relevant to heterogeneous CO2 reduction include refs 123, 168, and 169 discussed further in subsection 3.1.4. CMD is an enhanced sampling method,178−182 where CVs are dragged along a certain reaction pathway. Multiple windows then are inserted into the reaction trajectory, where MD steps are carried out with the CV fixed at each window to produce a potential of mean force (PMF).181−184 The free-energy profile can be obtained by integrating the PMF with respect to the CV along the trajectory. The simulation of a multidimensional (multiple CVs) reaction pathway (e.g., a concerted proton transfer involving multiple protons) can be computationally impractical for the CMD. A definition of an appropriate single CV describing the multidimensional reaction therefore is needed; an example of defining a single CV for a concerted proton-transfer process can be found in ref 123, which studied CO reduction on Cu(100). Just as for metadynamics, different assumptions/definitions of CVs may affect results (e.g., the freeenergy profile), which is a drawback of both methods. 1.3.7. Experimental Observables for Validating Theoretical Models. Several theoretical electrochemical models were introduced in subsections 1.3.3−1.3.6. Here, we will discuss experimental metrics for validating these models. Because most electrocatalytic reactions occur within the electrochemical double layer at electrode surfaces, we first introduce three experimental observables relevant to the simulation of the double layer’s properties. The first is the capacitance of the electrichemical double layer, which represents the variation of the electrode’s potential with respect to the surface charge. Early work111 from the Nørskov group simulating electrified solid−liquid interfaces (discussed in section 1.3.1) used the capacitance as a benchmark. The Arias group used the capacitance as a metric to validate their JDFT model as well.154 The second is the potential of zero charge (PZC), which is used to verify whether a given solvation model K


Chemical Reviews

Review

Scheme 3. Proposed Paths of Py-Catalyzed CO2 Reduction. Reprinted figure with permission from Senftle, T.; Lessio, M.; Carter, E.A. Chemistry of Materials, 2016, 28, 5799−5810.114 Copyright 2016 from the American Chemical Society

overpotentials were required to drive the reaction,192,193 suggesting that Py or a Py-derived species acts as a cocatalyst to reduce the overpotential. Based on an experiment on a Pt electrode, Bocarsly and coworkers proposed that the initial reaction step involves a oneelectron homogeneous reduction of pyridinium in solution (PyH+sol) to solvated pyridinyl radical (1-PyH•sol) at −0.58 V vs the saturated calomel electrode (SCE), followed by reaction with CO2 to form a carbamate complex (Scheme 3A).194 Subsequent experiments by the Bocarsly group showed cyclic voltammetry evidence demonstrating a one-electron reduction occurring on the Pt electrode.195 However, theoretical work of Tossell,196 Keith and Carter,134 Musgrave and co-workers,137 and Batista and co-workers197 excluded homogeneous mechanisms involving 1-PyH•sol, based on the very large discrepancy between computed reduction potentials to form solvated 1PyH• (−1.31 V to −1.58 V vs SCE) and the required applied potential on Pt in the experiment (−0.58 V vs SCE). The KeithCarter error estimate for the theoretical results is less than 0.2 V, when comparing experimental reduction potentials of a series of substituted PyH+ species.134 Therefore, a discrepancy of at least 0.7 V between the computed and experimental reduction potentials of PyH+ to 1-PyH• rules out the hypothesis of this homogeneous reduction mechanism at metal electrodes. Further investigations comparing semiconductor photoelectrode CBMs with a homogeneous PyH+ reduction potential were carried out by the Carter group.114,151 Their calculations

variations under transport-limited conditions are not considered. This deficiency affects the comparison between theoretical predictions and measurements. A recent multiscale simulation191 coupled a microkinetic model with a continuum model describing mass transport in an electrochemical cell. The coupled model significantly improved the qualitative comparison with experimental C−V curves after taking into account the local CO2 concentration and pH near the cathode surface. Thus, this work shows the promise of fully coupled multiscale simulations of the electrochemical interface, in which atomicscale, quantum-evaluated energy landscapes of catalytic steps couple to microkinetic-evaluated intrinsic overall reaction rates couple to continuum models of mass transport in the electrolyte.

2. HETEROGENEOUS CO2 PHOTOELECTROCHEMICAL REDUCTION 2.1. CO2 Reduction on p-GaP Photoelectrodes with Pyridine (Py) as a Cocatalyst

2.1.1. Proposed Catalytic Mechanisms Arising from Experiment and theory. Experimental observation of pyridine(Py)-catalyzed CO2 photoelectrochemical reduction at p-GaP electrodes was originally reported by Bocarsly and coworkers.24 They measured almost 100% faradaic efficiency of CO2 reduction to methanol at an underpotential of ∼0.3 V. Earlier experiments also showed good faradaic efficiency for CO2 reduction at the p-GaP photoelectrode. However, large L


Chemical Reviews

Review

Scheme 4. Paths of Py-Catalyzed CO2 Reduction via Surface-Bound (a) 2-PyH•* and (b) 2-PyH−*. Reprinted figure with permission from Senftle, T.; Lessio, M.; Carter, E.A. ACS Central Science, 2017, 3, 968−974.202 Copyright 2017 from the American Chemical Society

configuration at the GaP(110)/water interface,120 where a hydride-like species can be generated on the surface.121 Lessio et al.125 showed that protons adsorbed on GaP(110) are very stable, and while they may be converted to adsorbed hydride (H−*) under a negative bias and illumination, they will not be a proton source for DHP* formation. Such protons instead would have to come from the electrolyte solution. Lessio et al.’s subsequent work133 on hydricities of relevant possible intermediates in the Py*→DHP* reaction suggest that HT steps leading to DHP* and solvated formate (HCOO−sol) formation (Scheme 3D) are thermodynamically favorable. However, H−* transfer to any of the relevant species is predicted to be kinetically hindered (H−* to Py*/PyH+/CO2); therefore, H−* is unlikely to play a role in CO2 reduction at the p-GaP electrode. This suggested that either DHP* is formed via a more kinetically favorably pathway or another Py-derived intermediate might be the active catalytic species on the surface, such as 2-PyH•* and 2-PyH−*. The Carter group’s recent calculations170,202 showed a high HT barrier from DHP* to CO2 on the GaP surface. DHP* was therefore excluded as the active intermediate. After DHP* was proposed by Carter and co-workers to be a candidate catalytic intermediate for CO2 reduction, another mechanism involving DHP as the catalytic species was proposed by Musgrave, Hynes, and co-workers, whereby PyH+sol is homogeneously reduced to DHPsol (Scheme 3B).138 Reasons arguing against this mechanism are that (a) the first reduction step of PyH+ to 1-PyH• was already ruled out, even after considering the presence of the GaP photoelectrode;114,124,151 and (b) a strictly homogeneous mechanism cannot explain readily why CO2 reduction only occurs readily on certain photoelectrode surfaces and why its activity and selectivity is also electrode-dependent.26,29,30,205 As mentioned above, the Carter group’s work excluding DHP* inspired them to consider the possibility of 2-PyH•* and 2-PyH−* as catalytic intermediates. One heterogeneous mechanism proposes Py* to convert into 2-PyH•* (Scheme 4a).114,151 Their work on the stability and reduction potential of 2-PyH•* formation101,151,202 showed that 2-PyH•* is less stable and has a more negative reduction potential than 2-PyH−* due to the radical nature of 2-PyH• vs the closed-shell nature of 2PyH−. An alternative heterogeneous mechanism was proposed in which Py* is converted into adsorbed 2-PyH−* (Scheme 4b).101,133,170,202 The Carter group’s recent work shows that 2PyH−* is stable on the GaP surface, and that the reduction potential for 2-PyH−* formation via Py* + H+ + 2e− → 2-PyH−* is below the GaP CBM, independent of crystal facet.101,170,202

predict almost zero thermodynamic driving force for the PyH+ → 1-PyH• homogeneous reduction under the potential defined by the GaP CBM114,151 and also no PyH+ adsorption on the GaP surface.113,114,124 Their results indicate a much stronger driving force for the formation of Py* + H*, a possible intermediate state perhaps leading to formation of adsorbed dihydropyridine (DHP*).114,151 Py* is also a starting point for the formation of other Py-derived adsorbed species that might catalyze CO2 reduction (e.g., adsorbed 2-pyridinyl (2-PyH•*) and 2-PyH−* from Py*).114,151 A heterogeneous mechanism was proposed subsequently by Batista and co-workers (Scheme 3C)197 for Pt electrodes, and is supported by a number of measurements on Pt electrodes,195,198−201 showing that on Pt, PyH+ simply behaves as a weak acid, and instead of accepting an electron to form 1PyH•sol, it instead just donates its proton to the metal surface to form adsorbed hydrogen atoms (H*). Although this mechanism fails to explain why PyH+ in particular is needed in the semiconductor photoelectrochemistry, since its role is postulated to simply be that of a Brønsted acid, it inspired us to consider the role of adsorbed hydride on semiconductor photoelectrode surfaces.133,202 Keith and Carter proposed a different heterogeneous mechanism for Py-catalyzed CO2 reduction at p-GaP photoelectrodes. They predicted that the Py→DHP homogeneous two-electron reduction has a moderate potential (−0.72 V vs SCE),203 offering the first proposal of DHPsol as another possible intermediate that might form under electrochemical conditions (−0.58 V vs SCE on Pt). Bocarsly and co-workers suggested that DHP is not a reduction product because of the observed oneelectron reduction feature in cyclic voltammograms;195 however, their experiments were performed on a Pt electrode rather than a semiconductor photoelectrode. Despite the debate between theory and experiment, the earlier work of Keith and Carter on the Py→DHP homogeneous reduction potential inspired them to consider the possibility of DHP* as an intermediate. Carter and co-workers’ later calculations of Py and DHP adsorption energies and Py*→DHP* two-electron reduction potentials113,135,204 suggested that DHP* is a thermodynamically feasible intermediate on p-GaP that hence potentially could play a role as a catalytic intermediate. Keith and Carter proposed that Py* could be converted into DHP* (Scheme 3D),135,204 perhaps via a surface hydride. The Carter group’s earlier theoretical work121 showed that a full monolayer of half-dissociated water molecules is the most stable water M


Chemical Reviews

Review

both Py and PyH+ to be present, as would be the case at these intermediate pHs. In the mechanism proposed by Batista and co-workers (for the Pt electrode case),197 PyH+ first donates a proton to form H*, which then acts as an adsorbed hydride species (H−*) in the following step in which coupled HT/PT produces formic acid and eventually methanol. However, the most recent experiments206 on Pt suggest that hydrogen gas evolution and not CO2 reduction occurs, which is consistent with Batista and coworkers’ first step but not subsequent ones. Based on the Carter group’s predictions of hydricities and HT kinetics,133 transfers of H−* to Py*/PyH+sol/CO2sol are thermodynamically neutral/ favorable/favorable, respectively but also are kinetically hindered. Thus, H−* is unlikely to participate in catalytic CO2 reduction on p-GaP. This is further evidence excluding Batista’s mechanism for the p-GaP semiconductor photoelectrode. 2.1.4. Proposed Dihydropyridine (DHPaq vs DHP*) Intermediate. The Carter group first and then the Musgrave group proposed DHP as a candidate catalyst for CO2 reduction. Carter’s group considered DHP*, while Musgrave’s group suggested DHPaq in solution as the relevant intermediate, since homogeneous two-electron reduction of Pyaq to DHPaq has a moderate predicted potential (−0.72 V vs SCE).203 Since homogeneous mechanisms can be excluded, however, it inspired the theorists to spend more effort to test DHP* as an active catalytic intermediate. The Carter group’s studies demonstrated that DHP should be a stable adsorbate on p-GaP (110) and (111).114,124,151 Furthermore, hydricity calculations suggest that thermodynamically DHP* might be able to provide the hydride species to reduce CO2.133 However, a high kinetic barrier to form DHP* (the H−* transfer step is rate-limiting)133 and a high kinetic barrier for HT from DHP* to CO2170,202 indicate that formation of DHP* and reaction of DHP* with CO2 are both kinetically hindered. We therefore conclude that DHP* cannot be the active intermediate for CO2 reduction at p-GaP photoelectrodes. 2.1.5. Adsorbed 2-Pyridinide (2-PyH−*) as a Possible Intermediate. 2-PyH−*, the catalytic intermediate most recently proposed by the Carter group, is predicted to adsorb favorably on the p-GaP (110) and (111) surfaces, so it will be present at the electrode/electrolyte interface and not in solution.101,170,202 Moreover, a surmountable kinetic barrier (0.58 eV on GaP(110) and 0.19 eV on GaP(111)) for its formation via Py* + H+ + 2e− → 2-PyH−* is predicted.170 The predicted kinetic barrier for the subsequent HT from 2-PyH−* to CO2 is low (0.33 eV on GaP(110) and 0.59 eV on GaP(111)),101,170,202 and follows the measured trend of faster kinetics on the (110) facet205 compared to the (111) facet. Therefore, it can form and react as hoped. A subsequent question is whether it can survive long enough as a transient species to do so. Xu and Carter recently estimated the lifetime of 2-PyH−* against protonation to DHP or H2 by comparing the corresponding barriers to those of its formation and its transfer of hydride to CO2. They found that 2-PyH−* situated next to an adjacent OH−* should have sufficient lifetime to enable HT to CO2, and a high activation energy of HER via either 2-PyH−* or H−*.165 The high HER barrier is consistent with the observed selectivity. Until now, it was unclear whether 2-PyH−* could catalyze the subsequent CO2 reduction steps leading all the way to CH3OH, a critical last piece of the puzzle. Carter and co-workers very recently published207 a viable pathway, involving sequential HT steps to very specific CO2 reduction intermediates from 2-

Moreover, they demonstrated that the activation barrier for 2PyH−* formation (0.58 eV on GaP(110), 0.19 eV on GaP(111)) and HT from 2-PyH−* to CO2 (0.33 eV on GaP(110), and 0.59 eV on GaP(111)) are surmountable.170,202 The lower kinetic barrier for HT to CO2 from 2-PyH−* on GaP(110) compared to GaP(111) is consistent with the higher rates (current densities) measured for etched GaP(111) crystals exposing GaP(110) facets,205 providing further support for the involvement of this intermediate. Another factor to be taken into account is the lifetime of 2-PyH−*, since this negatively charged species is likely to be protonated to form DHP*, which as we have discussed has been determined to be an inactive adsorbate for CO2 reduction. The Carter group’s recent work165 demonstrated that 2-PyH−* situated next to an adjacent OH−* should have a sufficient lifetime for HT to CO2 to occur. Further, a high barrier to HER via either 2-PyH−* or H−* is predicted, consistent with the observed high selectivity toward CO2 reduction.24,192 Based on the above overview of the proposed mechanisms and the relevant experimental evidence, we next discuss the validity of these mechanisms and summarize a few conclusions in subsections 2.1.2 to 2.1.5. 2.1.2. Heterogeneous vs Homogeneous. Among the above six mechanisms, two of them are homogeneous catalytic reactions and the rest are heterogeneous. To narrow down the scope of our discussion, we discuss the possibility of heterogeneous vs homogeneous mechanisms in this subsection. The dependence of the CO2 photoelectroreduction activity and selectivity on choice of photoelectrode material, such as GaP/ CdTe/CuInS2,24,26,29,30 as well as on different facets,205 strongly suggests that the CO2 reduction goes through heterogeneous catalytic reactions. But does the electrode material choice necessarily mean that the reaction is happening on the surface? It could also mean that the CBM/metal Fermi level position is all that matters and that only certain electrodes have the appropriate CBM/metal Fermi level. However, examination of the CBM/metal Fermi level shows essentially no thermodynamic driving force for PyH+sol being reduced to 1-PyH•sol, both in the case of the Pt electrode under an applied potential of −0.58 V vs SCE134,137,196,197,203 and in the cases of semiconductor electrodes GaP and CdTe under illumination.114,151 Moreover, the significantly differing reaction kinetics for two different surface facets of GaP,205 each of which have essentially the same CBM,114,151 strongly suggests adsorbed species must be controlling the reaction kinetics. We therefore can safely conclude that the Py-catalyzed CO2 reduction at the p-GaP semiconductor surface is a heterogeneous process. 2.1.3. Debate about PyH+ as Merely a Proton Donor. Based on the predicted low acidity and high stability of adsorbed protons in a variety of environments on GaP(110),125 Lessio et al. concluded that they are not involved in CO2 reduction. The proton source therefore must come from a solvated protonated species in solution, either solvated, protonated water (H+sol) or PyH+sol. Experiments conducted on Pd electrodes18 measured no reduction phenomena associated with Py at pH > 7, which indirectly supports a role of PyH+ as proton donor (pKa = 5.3). However, the role of PyH+ cannot be merely to serve as proton donor because that would imply that any acidic environment in the absence of PyH+ should accelerate CO2 reduction, which is clearly false. Moreover, reduction activity was observed to decrease when the pH < 5.2.32 This observation, combined with the evidence of no reduction at pH > 7,18 suggests the need for N


Chemical Reviews

Review

Scheme 5. Formation of Surface Ti3+ with a Photogenerated Electron, Leading to the Bending and Binding of CO2 and the Formation of CO2•−. Reprinted figure with permission from White, J.L., et al. Chemical Reviews 2015, 115, 12888−12935.19 Copyright 2015 from the American Chemical Society. Note that the resonance structure on the right is incorrect (the C−O bonds should be equivalent)

PyH−*. Predicted barriers suggest that HT to HCOOH requires it to be adsorbed (HCOOH*) in order to react with 2-PyH−*, a new role for the electrode surface. HT to HCOOH* readily produces CH2(OH)2 upon protonation but subsequent HT to CH2(OH)2 to form CH3OH upon protonation is predicted to be hindered. However, CH2O, dehydrated CH2(OH)2, is found to react easily with 2-PyH−* to produce CH3OH upon protonation. Further reduction of CH3OH to CH4 via HT from 2-PyH−* encounters a high barrier, consistent with the observed absence of methane in experiments. The finding that the GaP surface plays a chemical role in enabling HT to HCOOH* further offers insight into why the primary CO2 reduction product at CdTe photoelectrodes is HCOOH rather than methanol: HCOOH is predicted to not adsorb on CdTe and so the reaction stops. This latest study thus suggests that 2PyH−* can catalyze the complete CO2 reduction pathway to CH3OH on GaP electrodes. Finally, the experimental evidence showing that CO 2 reduction activity drops when the pH < 5.232 is consistent with the postulate of 2-PyH−* as the catalytic intermediate, as lower pH would accelerate formation of the relatively inactive DHP*. Having ruled out essentially all other possibilities (1PyH•sol, Py*, H*, PyH+sol, DHPsol, DHP*, 2-PyH•*), we therefore conclude that the most likely catalytic intermediate by far is 2-PyH−*, at least on p-GaP surfaces. As we shall see in ections 2.2.4 and 2.2.5, this intermediate is predicted to be the only one thermodynamically accessible for other relevant semiconductors as well. The high probability of its protonation to form DHP is likely responsible for the overall low overall activity. In section 4.2, we discuss possible ways to overcome this challenge.

(LUMO) decreases in energy as the O−C−O bond angle decreases.33,220 A fairly recent theoretical study provides new insights into CO2 adsorption on the anatase-TiO2 (001) and (101) surfaces.126 CO2 adsorption energies were calculated with DFT using the hybrid XC functional PBE0. Experiments indicate the existence of a less stable (001) surface.221,222 This work126 demonstrated a stronger Lewis basicity of the surface oxygen site on the (001) surface compared to the (101) surface. The most favorable adsorption configuration is schematically shown in Figure 2. A strong coupling of the surface oxygen’s 2p

Figure 2. CO2 adsorption geometry of a monodentate carbonate (MC) configuration on the anatase (001) surface. Oxygen, titanium, and carbon atoms are represented by yellow, red, and blue circles, respectively. Reprinted with permission from Mino, L.; Spoto, G.; Ferrari, A.M. The Journal of Physical Chemistry C 2014, 118, 25016− 25026.126 Copyright 2014 from the American Chemical Society.

states with bent CO2 leads to the electron transfer from TiO2 to absorbed CO2. This study126 found weak CO2 adsorption on the (101) surface, as expected given the inert character of CO2. The comparison of CO2 adsorption activity on TiO2 (001) vs (101) could benefit from further consideration of the interaction between the surface Ti ions with CO2. For example, surface Ti3+ (formed by incident photons or oxygen vacancies) interacting with CO2 (Scheme 5) may lead to CO2 adsorption in the form of CO2•−, which was not considered in this study. Moreover, some adsorption configurations considered in this work involve physisorption of CO2. Thus, dispersion corrections should be included to accurately capture vdW interactions between CO2 and the surface. Others recently studied CO2 adsorption on anataseTiO2(101) in the presence of subnano-Ag/Pt clusters.127 The addition of subnanometer-scale metal clusters enhances catalytic activity of various reactions, e.g., CO oxidation223 and oxidative dehydrogenation of propane,224 due to structural fluxionality and a larger fraction of under-coordinated surface atoms.223,224 In this work,127 DFT calculations with the GGA-PBE XC functional were performed to obtain CO2 adsorption energies and charge distributions of CO2 adsorbed on the TiO2 surface. This work demonstrated that Ag/Pt octamer clusters on the

2.2. CO2 Reduction on Metal-Oxide and Metal-Chalcogenide Photoelectrodes

2.2.1. Titanium Oxide (TiO2). Several reviews comprehensively discuss CO2 reduction mechanisms on TiO2 photocatalysts.19,33,208,209 Anatase is the most utilized and studied phase due to its abundance and ease of synthesis compared to brookite, and its enhanced photocatalytic activity compared to rutile.19 CO and CH4 are the main products when TiO2 is in contact with steam, while methanol dominates if the reaction occurs in liquid water.208 An early proposed mechanism assumed that CO2 gains electrons from the TiO2 conduction band.210 Photoexcitation of TiO2 is thought to induce LMCT and produces Ti3+ (from O2−) on the TiO2 surface. CO2 then can be reduced by Ti3+, forming bent CO2•− (activated radical) adsorbed on the surface (Scheme 5).19,209 Experimental evidence exists for Ti3+ at the surface.211−213 Ti3+ can be formed by oxygen vacancies (n-type), dopants, or incident photons.211,214−219 CO2•− formation correlates with the extent of CO2 bending on the surface. The more bent the adsorbed CO2 is, the easier it would be for photoexcited electrons to transfer to CO2, because the CO2 lowest-unoccupied molecular orbital O


Chemical Reviews

Review

Figure 3. Proposed pathway for CO2 dissociation to CO on the surface of a reduced anatase-TiO2 (101)-supported Pt octamer. Reprinted with permission from Yang, C.T.; Wood, B.C.; Bhethanabotla, V.R.; Joseph, B. The Journal of Physical Chemistry C 2014, 118, 26236−26248.127 Copyright 2014 from the American Chemical Society.

Figure 4. Free energy profiles of CO formation from COOH* with and without a surface Ti3+ polaron, as computed by AIMD-based metadynamics within the GGA+U approach. Reprinted with permission from Klyukin, K.; Alexandrov, V. The Journal of Physical Chemistry C 2017, 121, 10476− 10483.142 Copyright 2017 from the American Chemical Society.

the presence of Ti3+ (Figure 4). This result is supported by an experimental observation of this reaction’s absence on a defectfree TiO2 surface and reactivity enhancement on a reduced surface with oxygen vacancies in the dark.227 Consistent with an earlier proposed mechanism for CO2 adsorption being activated by surface Ti3+,19,209 this new theoretical work demonstrated again the important effect of surface Ti3+ (induced by photons and/or surface oxygen vacancies) on heterogeneous CO2 reduction. Despite successfully explaining the enhanced catalytic performance of TiO2 in the presence of a Ti3+ polaron, the Ti3+ polaron model in this work is questionable. The authors added one extra electron to the system with a homogeneous compensating background to generate the Ti3+ polaron. The unphysical interaction of the extra electron localized on the Ti cation with an artificial homogeneous background charge is unavoidable using this approach and must be carefully vetted and the errors compensated for.102−104,108 A better model would be to introduce a surface oxygen vacancy, which in reality acts as an electron donor to generate a Ti3+ polaron, as was done in an earlier study.228 Moreover, although the activation energy for COOH* to CO* is predicted to decrease in the presence of a polaron (Figure 4), the barrier is still ∼1.2 eV, which is quite

TiO2 surface have a nonlocal effect on CO2 adsorption because of electron donation from the metal octamers to the TiO2 surface. The Pt cluster provides extra adsorption sites for bent CO2 due to hybridization between Pt d-states and CO2 molecular orbitals. A possible mechanism (Figure 3) was proposed to explain CO2 dissociation to CO* aided by the enhanced fluxionality of the Pt cluster due to oxygen vacancies at the TiO2 surface (i.e., a reduced surface). As above, it would have been better to include dispersion corrections in this study, as some of the stable configurations correspond to weak physisorption of CO2 on the surface. Besides theoretical work on anatase-TiO2, a recent AIMDbased metadynamics simulation studied CO2 adsorption and reduction to CO* on the rutile-TiO2 (110) surface.142 DFT Car−Parrinello molecular dynamics (CPMD) calculations were carried out using the GGA-PBE XC functional under a NoseHoover thermostat225,226 at T = 300 K. For the reactions involving a surface Ti3+ polaron, a GGA+U formalism71 was employed with U = 5 eV. It was reported that CO2 prefers to be activated and adsorbed at a bridging oxygen atom (Ob), forming a bent CO2− on the surface, followed by spontaneous protonation to form COOH*. This simulation suggested that further reduction from COOH* to CO* is greatly enhanced by P


Chemical Reviews

Review

(2000 μmol h−1 g−1) is about twice the rate at the (111) facet (1200 μmol h−1 g−1), and the faradaic efficiency is enhanced from 68% ((111) facet) to 77% ((112) facet). DFT-PBE calculations were carried out to explore the free energy of CO2* to CO* via COOH* on the (112) and (111) facets (Figure 6).232 The (112) facet was demonstrated to be much more

high for the reaction to proceed at room temperature. Measurements of the activation energy for this COOH* to CO* step are needed to validate the quantitative prediction of this work. 2.2.2. Indium Oxide (In2O3). Enhanced activity of the photocatalytic gas-phase RWGS reaction (CO2 + H2 + hν → CO + H2O) had been experimentally observed on a defective indium oxide surface (In2O3‑x(OH)y).229 A theoretical exploration of the mechanism was carried out by a cluster-model simulation,230 in which linear-response time-dependent DFT (LR-TDDFT) with the hybrid B3LYP XC functional were used for excited-state calculations. The unique atomic configuration of surface hydroxyls In−OH coupled with an adjacent O vacancy produces low-energy midgap states (1.5−2.1 eV above VBM) that could absorb visible light. The photon-induced electronic transition into these midgap states causes charge redistribution, rendering the under-coordinated surface In atom (adjacent to an oxygen vacancy, Vo) electron-deficient and the OH species electron-rich. The surface-frustrated Lewis pair (OHsurfaceδ‑ − [Vo] − Insurfaceδ+) gets more basic (referring to the OHsurfaceδ‑) and more acidic (referring to the Insurfaceδ+) upon photon activation, hence becoming more active for CO2 reduction (Figure 5). This study provides a unique example of

Figure 6. DFT calculations of adsorption and reduction of CO2 on Co3O4. Free energy diagram describing CO2* to CO* through the intermediate COOH* on the (111) and (112) surfaces of Co3O4 hexagonal platelets. Reprinted with permission from Gao, C., et al. Advanced Materials 2016, 28, 6485−6490.232 Copyright 2016 from WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.

active than the (111) facet. Because of the presence of undercoordinated Co ions at its surface, the (112) facet is relatively unstable compared to the (111) facet. It therefore makes physical sense that much stronger adsorption of CO2*, COOH*, and CO* are predicted for the (112) facet. From the adsorption geometries shown in Figure 6, we can see that multiple bonds are formed between the adsorbates and the (112) surface to saturate the under-coordinated Co ions (twofold-coordinated), while only one bond is formed between the adsorbates and the (111) surface, indicating a weaker interaction. The key argument in this work for a higher catalytic activity of the (112) surface is that the endoergicity going from CO2* to COOH* is 0.16 eV smaller on the (112) facet than compared to the (111) facet. This is based on the CHE model of the elementary PCET step CO2* + H+ + e− → COOH*, where the energy of (H+ + e−) is represented by the energy of H2 under standard conditions. However, this comparison is only valid if the catalytic reactions on the two different facets share the same free energy for the electrons involved. For the electrochemical CO2 reduction on different metal electrode surfaces, the applied potential determines the electron’s energy level, and thus the CHE model for an elementary PCET step is well-defined. For the photocatalytic reduction of CO2 on semiconductor photoelectrodes, the electron source is the CBM corresponding to a specific facet of the electrode. A more proper approach, i.e., to compare the thermodynamic driving force of a PCET step on different facets, is to calculate the reduction potentials (subsection 1.3.1) and compare the computed reduction potentials with the CBM levels of the corresponding facets. Because the CBM band-edge positions of the (112) and (111)

Figure 5. Schematic illustration of the origin of the difference in the activation energy for the reaction CO2 + H2 → CO + H2O involving the ground-state surface-frustrated Lewis pairs in the dark and the excitedstate surface-frustrated Lewis pairs under illumination. Reprinted with permission from Ghuman, K.K., et al. Journal of the American Chemical Society 2016, 138, 1206−1214.230 Copyright 2016 from the American Chemical Society.

the midgap state facilitating the photoelectrochemical reaction. The main reason is that the charge redistribution induced by the photoexcitation into this midgap state is localized on two separate surface species (Lewis acid and Lewis base), where the electron−hole recombination is hindered to some extent. On the other hand, relying on defect midgap states to drive chemistry is unlikely to be very promising, due to the small number of states available into which to excite electrons. Better to consider ways to alloy or dope materials to form bands of midgap states, as in intermediate band semiconductors that will provide a significant enough density of states to allow sufficient light absorption.231 2.2.3. Co3O4 Hexagonal Platelets. Recent experimental− computational joint work reported a high activity of CO2 photocatalytic reduction to CO on Co 3 O 4 hexagonal platelets.232 The rate of CO production at the (112) facet Q


Chemical Reviews

Review

electrons from the CBM and a proton from solution) on CdTe (111). This was suggested as a possible explanation as to why the main CO2 reduction product over CdTe is HCOOH (lessreduced product compared to CH3OH), given the smaller driving force compared to GaP or CuInS2 (vide infra). However, recent work207 in the Carter group points to yet another explanation, namely, that HCOOH does not adsorb on CdTe, rendering the subsequent HT to HCOOH kinetically forbidden (discussed above in section 2.1.5). 2.2.5. Copper Indium Sulfide (CuInS2). Very promising Py-catalyzed CO2 photoelectroreduction with a high selectivity (97% to methanol) and a low overpotential (20 mV) was observed on the CuInS2 photoelectrode.26 The Carter group used similar theoretical schemes to those done for GaP and CdTe to calculate the CBM of the reconstructed CuInS2(111)/ water interface, the reduction potential associated with 2-PyH−* formation, and the activation energy for HT from 2PyH−* to CO 2 on CuInS2(112).202 Unlike CdTe(111), a strong thermodynamic driving force was predicted for 2-PyH−* formation on CuInS2(112). Moreover, activation energies for HT from 2-PyH−* to CO2 on the three relevant photoelectrode surfaces (GaP(111), CdTe(111), and CuInS2(112)) are all surmountable ( (111). 190,236,237 In the following subsections 3.1.1−3.1.4), we discuss CO2 reduction mechanisms on Cu surfaces as predicted by various electrochemical models. 3.1.1. Mechanisms Predicted by the CHE Model. The CHE model for studying CO2 electrochemical reduction on Cu(211) was conducted originally by DFT-GGA calculations using the revised-PBE (RPBE)238 XC functional.159 The predicted PCET steps on Cu(211) are CO2 → COOH* → CO* → CHO* → CH2O* → CH3O* → CH4↑ + O*. The PDS is CO*→CHO* (Figure 11). For technical details of the CHE model, please refer to the previous subsection, 1.3.3. 3.1.2. Linear-Scaling Relationship and Overpotential Volcano Plot. Intermediates that bind to surface Cu atoms through the C atom have binding energies (here obtained from DFT with the RPBE238 XC functional) that linearly scale with each other, e.g., CO*, COOH*, CHO*, and CH2O* (Figure 12).128 This phenomenon can be explained by the d-band center theory, discussed in the next paragraph. Based on the linearscaling relationship, using the CO* binding energy EB(CO*) as a descriptor, a volcano plot emerges for the overpotential of the CO2 reduction to CH4 with respect to EB(CO*) on different metal electrodes (Figure 13).128 To understand Figure 13, please refer to section 1.3.3 where we introduced the CHE model and how the UL and the overpotential are calculated. With the linear scaling relationships for the binding energies of CO2 reduction intermediates (Figure 12), the UL (or ΔG) of a PCET reaction on various metal electrodes can be described by a linear function with respect to the binding energy of a single intermediate, which here is the binding energy of CO, EB(CO*), as shown in ref 128. Since each PCET step has a linear relation UL(EB(CO*)), if we draw the UL(EB(CO*)) lines of all the involved PCET steps (CO2 reduction to CH4) on a single plot, the most negative envelope gives us the thermodynamically least-required potential to drive the reaction. The equilibrium potential of CO2 reduction to CH4 (horizontal dashed line in Figure 13) is derived based on the reaction equation CO2(g) + 8H+ + 8e− → CH4(g) + 2H2O. Similar to the above definition of UL, the thermodynamically required potential for driving this reaction is +0.17 V (vs RHE). The overpotential of CO2 reduction to CH4 on a certain metal electrode (represented by a certain EB(CO*)) can be obtained by 0.17 − UL(EB(CO*)) (V). The final results are shown in Figure 13; the overpotential volcano is depicted by the CO*→CHO* red line and the CO2→ COOH* blue line with the intersection at about EB(CO*) =

Figure 9. Proposed photoelectrocatalytic cycle for CO2 to CO reduction on pyridinic N (blue atom)-doped coronene (C24H12). Reprinted with permission from Yang, K.D., et al. Advanced Functional Materials 2016, 26, 233−242.234 Copyright 2015 from WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim.

reduction to CO was proposed based on their calculations of CO2 adsorption at different N-defect sites in N-GQS.234 Their proposed catalytic pathway of CO2 to CO is CO2 → CO2* → COOH* → CO + H 2 O, is the same as that for CO 2 photoelectroreduction to CO on Co3O4(112)232 and for CO2 electrochemical reduction to CO*/CO on a Cu electrode159 (which has been extensively studied and reached consensus; see subsection 3.1.4).

3. HETEROGENEOUS CO2 ELECTROCHEMICAL REDUCTION 3.1. CO2 Reduction Mechanisms on Cu Surfaces

The high electrical conductivity and intrinsic catalytic activity of metal cathodes makes them attractive targets for investigation. Among all metal cathodes, Cu is the only one that exhibits high selectivity for CO2 electrochemical reduction to hydrocarbons (CH4 and C2H4), as observed long ago in the experiments of Hori and co-workers.5,6,8−10 The corresponding data are graphically summarized in Figure 10. Consequently, the mechanism of CO2 reduction at Cu cathodes is the most intensively studied topic in the field of CO2 heterogeneous electrochemical reduction.9−12,14,17,36 A review article in 201536 comprehensively discussed and summarized the progress in understanding heterogeneous CO2 electrochemical reduction mechanisms from theory. We briefly summarize the major points from that review and then discuss progress since 2015. CO2 electrochemical reduction occurs on three different facets of Cu: (211), (111), and (100). Under standard electrochemical conditions, polycrystalline Cu transforms to Cu(111) in 30 min and Cu(111) transforms to Cu(100) in 60 min without any further changes,235 indicating a higher stability of the (100) facet. Cu(211) is a stepped facet with undercoordinated Cu atoms exposed at the surface, which is less stable than Cu (100) and (111). Therefore, the corresponding stability order is (100) > (111) > (211) under standard electrochemical conditions. However, the activity order of the facets is (211) > S


Chemical Reviews

Review

Figure 11. Free energy diagrams for the lowest energy pathways to (a) H2, (b) HCOOH, (c) CO, and (d) CH4. In each diagram, the black (higher) pathway represents the free energy at 0 V vs RHE, while the red (lower) pathway represents the free energy at the indicated potential. Reprinted with permission from Peterson, A.A.; Abild-Pederson, F.; Studt, F.; Rossmeisl, J.; Nørskov, J. K. Energy & Environmental Science 2010, 3, 1311−1315.159 Copyright 2010 from The Royal Society of Chemistry.

−0.6 eV, which defines the peak of the volcano (corresponding to the lowest overpotential). Cu is the closest to the volcano’s peak among all metal electrode materials. The d-band center theory239,240 clarifies the fundamental reason for the above linear-scaling relationship between the adsorption energies of different adsorbates (sharing the same binding element, e.g., C atom among CO*, COOH*, CHO*, and CH2O*). The main idea is that the closer to EFermi the metal d-band center is, the stronger the adsorbate-metal interaction is due to a lower occupation of antibonding states.239 The d-band center theory provides a simple and often useful descriptor for the binding energies of various adsorbates on metal electrodes. Here, we expound on the limits of the d-band center theory. The adsorption energy of an adsorbate can be divided into two parts, ΔEads = ΔEsp + ΔEd, where ΔEsp is the binding energy contributed from the interaction between the adsorbate’s valence level and the metal sp-states, and ΔEd is the extra interaction with the TM d-states.241 The d-band center model assumes that ΔEsp is independent of the metal. This approximation may fail when the metal particle gets small enough (typically diameter −1.0 V vs RHE), making sure that variations in gas flow rates did not affect product distributions. Periodic DFT calculations with the Bayesian error estimation functional (BEEF)−vdW251 were carried out to evaluate the cations’ effect on the energetics of CO2 reduction intermediates. The cations were found to act as promoters via the electric field induced by the cations in the OHP (consistent with an earlier theoretical study252). The electric field significantly stabilizes CO2* and surface species containing C−C bonds, thereby enhancing the selectivity toward C2 products. By contrast, the intrinsic electric field induced by a cation nearer to the surface was found to vary unsystematically with cation size. They proposed that an increased concentration of larger cations in the OHP is the critical reason for enhanced C2 product formation, supported further by DFT predictions that larger cations prefer to locate near the Cu surface rather than in the bulk electrolyte. Several months after this work was published, another joint experimental-theoretical study253 by Koper and co-workers revealed a similar cation promoter effect, also elucidating how

Figure 18. Most favorable kinetic pathways for the electrochemical reduction of CO2 to CO* and formate (HCOO−) on the Cu(100) surface in water at pH = 7. Reprinted with permission from Cheng, T.; Xiao, H.; Goddard, W. A., III Journal of the American Chemical Society 2016, 138, 13802−13805.168 Copyright 2016 from the American Chemical Society.

With the same methodology, the Goddard group also studied the variation of the CO* reduction product distribution with respect to an applied potential on Cu(100).169 They reported that at an applied potential U less negative than −0.6 V (vs the RHE), the coupling of two CO* forming a C−C bond has a lower barrier than that of CHO* formation; thus, C2H4 is the dominant product at U > −0.6 V. When −0.8 V < U < −0.6 V (vs the RHE), H* displaces CO* as H* is more strongly bound at this potential. The predominant H* species block further CO* reduction, so C2H4 production decreases and H2 production increases in the −0.8 V < U < −0.6 V range. When U < −0.8 V (vs the RHE), the blocking effect of H* on the hydrocarbon production is mitigated and CHO* can form from nonadsorbed CO directly reacting with H*, leading to enhanced CH4 production (Figure 19). C2H4 formation also increases by

Figure 19. Reactive trajectories of *CHO formation for a high coverage of surface H (1 ML of H*) from AIMD simulations at 298 K for (A) nonadsorbed CO, (B) TSs, and (C) *CHO. Reprinted with permission from Cheng, T.; Xiao, H.; Goddard, W. A., III Proceedings of the National Academy of Sciences 2017, 114, 1795−1800.169 Copyright 2017 from the National Academy of Sciences of the United States of America.

nonadsorbed CO reacting with CHO*. This study demonstrated the importance of H* for the CO reduction pathway with respect to applied potential. As CO* is an intermediate in the CO2 reduction pathway, this study has important implications for CO2 reduction on Cu as well. It was demonstrated in this work that at a moderately negative potential range −0.8 V < U < −0.6 V (vs the RHE), H* inhibits CO* reduction to X


Chemical Reviews

Review

the effect changes with applied potential and crystal facet. Cation size was shown to have a stronger impact on the onset potential for C2H4 formation at low overpotentials. Potentials more negative than −0.65 V (vs RHE) favored CH4 production and the corresponding onset potential was found to be independent of cation size. Cation promotion of C2H4 production was more significant on Cu(100) than on Cu(111) or polycrystalline Cu, consistent with an earlier observation254 that Cu(100) favors conversion of CO2 to C2H4. Their DFT-PBE calculations revealed that cations more strongly stabilize adsorbates containing a C−C bond (OCCO* and OCCOH*) compared to C1 intermediates (CHO*), supporting the proposed intrinsic cation promoter effect for CO2 reduction to C2 products.250,253 Electrolyte anions also influence the activity and selectivity of CO2 reduction. Progress in mechanistic understanding of the anion’s role was reviewed recently,255 where buffering,8 stabilization of intermediates,256 and surface roughness257 were highlighted. Recent Raman spectroscopy experiments255 found that coadsorbed halide ions (Cl−, Br−, I−) enhanced the amount of CO adsorbed on Cu, consistent with earlier theoretical predictions that the CO* binding energy can increase slightly in the presence of coadsorbed anions.258 The resultant higher CO* coverage may assist C−C coupling, a critical step in C2 production. Furthermore, the authors proposed that coadsorbed anions altered the local charge state of the CO*. They claimed that the coadsorbed anions create a more positively charged carbon atom in CO*. If a nearby surface site is anion-free, they claimed that the carbon adsorbed on the anion-free site would be more negatively charged due to the πback-donation from Cu. Two oppositely charged CO* species then could form more readily the adsorbed CO dimer (OCCO*), enabling C2 product formation, which is consistent with Goddard and co-workers’ earlier prediction.259 However, we would argue that coadsorbed anions should induce stronger π-back-donation to CO*, building up more, not less, negative charge on the carbon atom in CO*, and weakening the C−O bond. When an anion adsorbs on a metal surface, bringing its own negative charge from solution, the adsorbed anion will donate some negative charge to the metal, rendering the latter even more electron-rich, enabling the increased π-backbonding mentioned above. Actually, their in-situ Raman spectra255 supports our argument. The authors report that the frequency of the C−O stretch (in CO*) changes from 2087 to 2060 cm−1 (indicating bond weakening) as the electrolyte changes from ClO4− → Cl− → Br− → I−. They also showed that the binding abilities of these anions on Cu increase in the order of ClO4− → Cl− → Br− → I−. Thus, the experimental evidence instead indicates that more coadsorbed anions nearby induce a stronger π-back-donation, which is likely the source of the greater reactivity due to the weakening of the C−O bond. Despite the controversy about the exact charge-state variation of CO* with respect to the anion adsorption, this study255 demonstrated that anions also play an important role in catalyzing C2 production, by tuning the CO* coordination environment and hence their charge, C−O bond strength, and binding affinity on electrode surfaces.

Figure 20. Reaction pathways for the electroreduction of CO2 to different products (top) and the competing HER (bottom). Reprinted with permission from Yoo, J.S.; Christensen, R.; Vegge, T.; Nørskov, J. K.; Studt, F. ChemSusChem 2016, 9, 358−363.260 Copyright 2016 from Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim.

suggested that CO2 reduction via HCOO* (rather than COOH*) is more promising for selectivity toward HCOOH (vs HER), because it may break the linear-scaling relationship between H* and COOH*. To be clear, HCOO* is bound to the surface via both oxygen atoms whereas COOH* is bound to the surface via the carbon atom (Figure 20). The general rule proposed is to find a catalytic electrode binding HCOO*/H* strongly/weakly. Ag and Pb were found to be the best candidates for CO 2 reduction to HCOOH through the HCOO* intermediate while binding H* weakly to suppress HER. Based on Goddard’s work discussed in subsection 3.1.4, the mechanism of nonadsorbed CO2 directly reacting with H* (along with a one-electron transfer) and forming formate should be considered here, as well. Again, Goddard and co-workers suggested168 from their kinetics analysis that formate is produced mainly through physisorbed CO2 reacting with H* on Cu(100). Although the reduction potential for the HER is less negative than the reduction potential for HCOOH formation, indicating that thermodynamically H* prefers the HER rather than reducing CO2 to HCOOH, the difference in the reduction potentials is only 0.17 V, so the conclusion might change if kinetic barriers are considered. It is possible that a moderately bound H* might enhance HCOOH production through physisorbed CO2 reacting with H*, if the corresponding kinetic barrier is lower than that for the HER (H* + H+ + e− → H2). Recent experimental-theoretical joint work studied Bi dendrites as catalytic electrodes for CO2 reduction to formic acid in an aqueous electrolyte.261 The Bi dendrite structure has the advantages of a large surface area and a large number of under-coordinated surface sites. In the experiment, the Bi dendrite exhibited a current density of 2.7 mA/cm2 with a faradaic efficiency of 89% for HCOOH at −0.74 V (vs the RHE), while a pristine Bi foil exhibited a current density of only 0.4 mA/cm2 with a faradaic efficiency of 56% at the same potential. The CHE model within DFT-PBE calculations on various Bi planes, suggested that the CO2 → OCOH* → HCOOH* path gives the lowest overpotential among all possible PCET pathways (Figure 21).261 Note that the intermediate referred to as OCOH* here is bound to the surface via the two oxygen atoms, exactly the same as the “HCOO*” intermediate discussed above.260 The high-index Bi dendrite facets (e.g., the (012) and (104) facets) studied satisfy the metal electrode screening rule for CO2 reduction to HCOOH proposed by Yoo et al.,260 where one needs a catalyst

3.2. CO2 Reduction Mechanisms on Other Metal (non-Cu) Electrode Surfaces

CO2 electrochemical reduction to HCOOH vs to CO vs the HER on a variety of metal surfaces was studied via the CHE,260 with periodic DFT employing a BEEF−vdW XC functional.251 Figure 20 displays the mechanisms considered. This work Y


Chemical Reviews

Review

Figure 21. DFT predictions of the CO2-to-formate reduction reaction pathways for various Bi planes: (003), (012), (110), and (104). (a) Schematic diagram of three possible reaction pathways: Path 1 via formation of *COOH, Path 2 via formation of *OCOH, and Path 3 via formation of *H. Reaction free energy diagrams for (b) Path 1, (c) Path 2, and (d) Path 3, when a bias potential U = −0.21 V (vs SHE) is applied. Note that the (012), (110), and (104) surfaces grow preferentially in the Bi dendrite compared with the close-packed (003) surface dominant on the Bi foil. Reprinted with permission from Koh, J.H., et al. ACS Catalysis 2017, 7, 5071−5077.261 Copyright 2017 from the American Chemical Society.

surface. The authors proposed three criteria for NSA A/B pair screening: (1) an optimal CO* binding energy; (2) a weak H* binding energy; and (3) weak OH* binding for selective production of CH3OH. Periodic DFT-RPBE calculations suggested that the W/Au subsurface alloy is the most promising NSA pair. A free-energy diagram based on the CHE model compared the energetics of CO2 reduction to CH3OH/CH4 on Cu vs W/Au (Figure 23). W/Au subsurface NSA has a lower overpotential for CO* to CHO*, a high H* binding energy (suppressing HER), and weak OCH3 binding, so CH3OH can be formed. However, the CO2 reduction mechanism on Cu(211) was later modified190 to be CHO* → CHOH* → CH* → CH2* → CH3* → CH4. This change negates the third criterion above because it is based on an outdated mechanism in which OCH3* is the precursor to both CH4 and CH3OH production. In contrast to extensive research into CO 2 electrochemical reduction to methane, further exploration of possible mechanisms for CO2 electrochemical reduction to methanol on metal electrodes is warranted. Finally, a combined experimental-theoretical study showed that a Au surface tuned by Cu enrichment could give a continuous product distribution varying from CO (produced by CO2 electrochemical reduction) to H2 (produced by HER).263 Periodic DFT-PBE calculations were conducted to explore the variations of CO adsorption energies and C−O vibrational frequencies with respect to the fraction of Cu substituted for Au in the Au (111) surface layer. Calculated vibrational frequencies matched well with Raman spectra, and suggested that higher Cu substitution fraction on the Au surface binds CO more strongly and increases the competition between HER and CO desorption. The stronger CO binding to the Cu-enriched Au surface also can be explained by the d-band center theory

binding HCOO* strongly while binding H* weakly (Figure 21c,d). The Bi electrode was not considered by Yoo et al.260 This study supplements that work260 and supports the idea of searching for catalytic materials stabilizing HCOO* as an intermediate species for the CO2 electrochemical reduction to HCOOH. 3.3. CO2 Reduction Mechanisms on Metal Alloy Electrode Surfaces

A theoretical investigation of CO2 reduction to methanol was conducted on a near-surface alloy (NSA) electrode in search of promising NSA pairs A/B (where both the host A and solute B are metal elements).262 The geometries of the overlayer alloy (A*/B) and subsurface alloy (A/B) defined in this work are shown in Figure 22. This study focused on the (211) stepped

Figure 22. Schematics for (211) stepped overlayers (A*/B) and subsurface (A/B) alloys. Red and black atoms denote the solute (A) and host (B) metals, respectively. Reprinted with permission from Back, S.; Kim, H.; Jung, Y. ACS Catalysis 2015, 5, 965−971.262 Copyright 2014 from the American Chemical Society. Z


Chemical Reviews

Review

Figure 23. Free-energy diagram for (A) CO2 electroreductions to CH4 or CH3OH (red), and (B) H2 evolution reactions at zero applied potential (U = 0 V vs RHE) for Cu (black) and W/Au (blue). Each reaction step involves a proton−electron pair (H+ + e−) transferred to the adsorbate on the electrode, where the proton/electron comes from the solution/electrode. Reprinted with permission from Back, S.; Kim, H.; Jung, Y. ACS Catalysis 2015, 5, 965−971.262 Copyright 2014 from the American Chemical Society.

(subsection 3.1.2). The d-band center of pure Cu(111) is closer to EFermi than that of pure Au (111). Likewise, this work263 showed that the calculated d-band center of the Cu-substituted Au (111) approaches EFermi as the Cu fraction increases. A metal surface with a d-band center closer to EFermi leads to a stronger interaction with adsorbates.239,240 That is why a stronger binding of CO was observed on the Au surface with a higher Cu substitutional fraction in the experiment.263

effects. A higher mass-transport rate may overconsume protons near the cathode, yielding a higher pH, while certain buffering cations (Cs+, Rb+)249 and anions (HCO3−, HPO42−) may alter the local pH and mitigate the pH increase. The second factor is the potential applied to the electrode. Since CO2 electrochemical reduction involves electron transfer from the electrode to adsorbed intermediates, variations in applied potential may affect the activation energies of the C1 vs the C2 pathways, and hence the product distributions may be potential-dependent. A double-peak feature appeared in the measured C−V curve for CO (the common intermediate in the C1 and C2 pathways) reduction to C2H4 on Cu(100),236 with selectivity toward C2H4 (CH4 and H2 are the two other reduction products) enhanced at low (−0.4 V to −0.6 V vs RHE) and high overpotentials (more negative than −0.9 V vs RHE). Two recent theoretical studies169,185 provided insights into the mechanism. The common finding is that a two-pathway mechanism is necessary to explain the observed potential dependence. Both studies predicted that C2H4 forms via CO* dimerization at a low overpotentials while CO reacting with CHO* contributes to C2H4 formation at high overpotentials. However, the predicted mechanistic details differ. Bell, HeadGordon and co-workers185 suggest that C−C bond formation at high overpotentials proceeds with adsorbed CO and CHO whereas Goddard and co-workers argue169 that C−C bond formation occurs via a CO molecule from solution reacting with an adsorbed CHO* on the electrode surface. Despite the controversy about which form of CO is involved, the experimental evidence combined with theoretical insights demonstrate that varying the applied potential can affect the selectivity toward the C2 product. The third factor is the exposed surface facet. The energy landscape of a catalytic pathway obviously depends on the surface atomic structure. The Nørskov group predicted that surface-bound CO dimerization can occur on Cu(100)266 but not Cu(211),267 consistent with the experimental observation that Cu(100) exhibits high selectivity toward C2H4 prodcution.254 The fourth factor is the effect of dissolved cations and anions near the electrode surface. Section 3.1.5 reviews recent progress in understanding electrolyte ion promoter and buffering effects. The fifth factor is the influence of surface metal atom oxidation state. Extensive experimental evidence exists that the

3.4. Mechanistic Understanding of C2 Product Formation

Until now, we have emphasized insights gleaned from theory regarding C1 product (CH4, CO, and HCOOH) formation. However, developing ways to enhance selectivity toward multicarbon products is of great importance because of their higher energy density. We therefore discuss next the formation of C2 products (primarily C 2 H4 and C2 H5 OH); little mechanistic information about C3+ production is available to review. A recent article summarized experimental observations of C2 product formation,253 including early studies by Hori, and Koper, among others. Another article recently reviewed theoretical insights, particularly regarding the proposed ratelimiting step of C−C bond formation,185 including those of Koper, Nørskov, Bell, Head-Gordon, and Goddard. Here, we present a high-level discussion of the physical and chemical factors affecting C2 selectivity during CO2 electrochemical reduction at metal electrodes. The first factor to consider is the local pH value. Early experiments264 showed that the onset potential for CH4 formation depends on pH but C2H4 production does not. The suggested origin of this difference is that CH4 formation involves a series of PCET steps whose rates have an obvious dependence on the local pH near the cathode, whereas the critical step for C2H4 evolution is C−C bond formation, which does not involve proton transfer.169,185,265 A recent theoretical study247 from the Goddard group discussed the pH dependence of CO2 reduction pathways and product distributions. CO2 reduction was predicted to proceed via a C1 pathway (CO* → COH* → CHOH*) at low pH on Cu(111). At neutral pH, C−C bond formation via CO* reacting with COH* was predicted to dominate. As the pH further increases to pH = 12, CO* dimerization makes a major contribution to C−C bond formation. Note that the local pH near the cathode is not only determined by the initial experimental setup but can be affected by other factors, e.g., the mass-transport rate and ion buffering AA


Chemical Reviews

Review

4. DISCUSSION AND OUTLOOK

activity and selectivity toward C2 products of Cu electrodes are enhanced if Cu is first oxidized and then reduced.268−272 Thus, the oxidation state of the surface Cu atoms in these oxidederived Cu electrodes may play an important role in catalyzing CO2 reduction. A recent theoretical study259 from the Goddard group predicted that CO2 activation and CO* dimerization is enhanced by coupling Cu+ and Cu0 sites in a metal embedded in an oxidized matrix (MEOM). This conceptual model was used to mimic a case in which both Cu+ and Cu0 sites are present on the surface of oxide-derived Cu electrodes under relatively reducing conditions. The authors found that CO2* on the Cu0 site can be stabilized by H2O* on an adjacent Cu+ site, leading to enhanced CO2 activation. For CO* dimerization, the carbon atoms in CO* adsorbed on Cu0 and Cu+ sites were predicted to be oppositely charged, with Cu0/Cu+ inducing a negatively/ positively charged carbon atom. The electrostatic attraction of two oppositely charged CO* adsorbates facilitates CO dimerization. A subsequent joint theoretical−experimental study273 explored the potential effect of subsurface oxygen in Cu(111) electrodes. They claimed that a moderate amount of subsurface oxygen (identified by ambient pressure X-ray photoelectron spectroscopy) could stabilize CO2 adsorption via adsorbed H2O* on the Cu+ induced by the subsurface oxygen. The mechanism of this CO2 activation effect is analogous to that reported in the Goddard group’s earlier study259 on the mixture of Cu+ and Cu0 sites in the MEOM system. However, new experiments274 raised a critical question as to whether the oxygen identified in the earlier study273 was oxygen retained from the preoxidation of Cu, or the oxygen was introduced during the rinsing process (for removal of reduced Cu) or during the steps to prepare the sample for spectroscopic analysis. A subsequent theoretical study275 from Bell, HeadGordon, and co-workers predicted that oxygen atoms are not stable in the subsurface layer. The calculated kinetic barrier for a subsurface oxygen migrating to the surface is sufficiently low (∼0.35 eV) that subsurface oxygen atoms will not persist during CO2 reduction. The retained oxygen atoms in the oxide-derived Cu electrodes instead prefer to locate on the Cu surface. The free energy of an oxygen atom at a hollow site on the Cu surface was predicted to be ∼2 eV lower than that of an oxygen atom at a subsurface site. The authors concluded that the subsurface oxygen does not play a role in catalyzing CO2 reduction due to its instability. These authors instead attributed the observed enhancement in CO2 reduction on the oxide-derived Cu electrodes to new surface defects (such as vacancies and under-coordinated surface metal atoms), supported by recent investigations of CO2 reduction on Cu and Ag.276,277 Despite the controversy about the stability of subsurface oxygen, the key idea of Cu+ and Cu0 sites at Cu electrode surfaces working in concert to enhance CO2 adsorption and CO* dimerization259 inspired researchers to design novel strategies for tuning the oxidation states of surface atoms in Cu electrodes. A very recent study278 proposed a boron-doping approach to achieve a system with both Cuδ+ and Cu0 regions present at the Cu electrode surface, exhibiting an analogous pattern to the MEOM system investigated in the Goddard group work.259 A high selectivity (∼80%) toward C2 products via CO2 reduction was observed on the boron-doped Cu(111) electrode, with an average Cu valence state of +0.35. This finding demonstrated the important effect of surface metal atom oxidation states on CO2 reduction pathways.

4.1. Possible Intermediates for Py-Catalyzed CO2 Photoelectrochemical Reduction on Semiconductor Surfaces

Based on the previous detailed mechanistic discussion of Pycatalyzed CO2 photoelectrochemical reduction on p-GaP electrodes (section 2.1), we can conclude that the catalytic reaction is heterogeneous and the active catalytic intermediate has to be an adsorbed Py-derived species. The closed-shell species DHP* and 2-PyH−* are predicted to be stable on GaP surfaces.114,124,135,202 The radical species 2-PyH•* is much less stable than 2-PyH−*, so 2-PyH•* can be excluded.151,170 Comparing DHP* and 2-PyH−*, DHP* is much less likely to be the active intermediate due to a high kinetic barrier for HT to CO2,170,202 while for 2-PyH−*, a low kinetic barrier for HT from 2-PyH−* to CO2 is predicted for both GaP (110) and (111) surfaces, with the trend in relative barriers consistent with experiment.170,202 Moreover, the predicted activation barrier for 2-PyH−* formation on GaP surfaces is surmountable.170 All of the above evidence strongly suggests that 2-PyH−* is the active catalytic intermediate. However, as the negatively charged species 2-PyH−* could be protonated easily (by protons in the aqueous electrolyte) to form the inactive DHP*, it is critical to establish at least qualitatively the lifetime of 2-PyH−*. Another relevant issue is whether 2-PyH−* can undergo the HER. Since 2-PyH−* transfers hydride to CO2 facilely, it is necessary to understand whether HT from 2-PyH−* to a solvated proton (forming H2) occurs easily or not. Xu and Carter recently predicted165 that 2-PyH−* situated next to OH−* should have a sufficient lifetime for HT to CO2 to occur. Moreover, a high activation energy for the HER via either 2-PyH−* or H−* was predicted. With all other possible Py-related species having been excluded, all pieces of evidence point to 2-PyH−* as the active catalytic species, at least for the first reaction step of CO2 reduction to formate. The question then remains whether the same holds true for the subsequent reduction steps taking formate to methanol. Excitingly, the Carter group’s most recent work207 shows that 2-PyH−* is a feasible catalyst for reducing CO2 all the way to CH3OH via a series of HT steps via specific intermediates (discussed in detail in section 2.1.5). 4.2. Theoretical Insights for Cocatalyst Design for Heterogeneous CO2 Photoelectrochemical Reduction

HT to CO2 is the key step in CO2 reduction. Without adsorbed cocatalysts (e.g., Py-derived species adsorbed on the electrode surface), either overbinding or underbinding of hydride on an electrode surface impedes HT to CO2.133,202 The main role of a cocatalyst is to optimize HT by forming an active intermediate with the following properties: (1) The catalytic species is stable on the electrode surface. (2) The required potential for catalytic intermediate formation is less negative than the potential corresponding to the photoelectrode CBM. (3) The activation energy to form the catalytic intermediate is low. (4) The activation energy of HT from the catalytic intermediate to CO2 is low. (5) The catalytic intermediate has a sufficiently long lifetime against certain side reactions with electrolyte molecules or dissolved ions in the electrolyte. Based on the above principles, we propose some ideas for improving the Py-catalyzed CO2 reduction on chalcogenide semiconductor surfaces. One possible approach is surface substitutional (cation) doping of the surface, where the electronegativity of the dopant metal atom can be a tunable knob. A doping cation with large electronegativity tends to AB


Chemical Reviews

Review

oxygen bound at the metal site), where the reduction of OCHO* leads to HCOOH production. The lack of a scaling relationship between H* and OCHO*presumably originating in the dissimilar character of metal-H and metal-O bonds provides opportunities to select metals that alter the energetics of OCHO* and H* independently; for example, Ni and Cr were found to bind OCHO* more strongly than H*.281 Furthermore, OCHO* was predicted to be preferentially stabilized over H* in this metal-porphyrin-like system when the single-metal active site is coordinated with an axial ligand, leading to enhanced HCOOH production and suppressed H2 evolution. The authors explicitly calculated the overpotentials of forming H2, CO, and HCOOH using the CHE model. The predicted overpotential for forming HCOOH is significantly lower than for forming H2 when the CO2 reduction intermediate OCHO* becomes more stable than H*. This indicates that the critical step for CO2 reduction in this porphyrin-like system is the first reduction step (forming the intermediates H*, COOH*, or OCHO*) rather than subsequent reactions. Aside from breaking the linear-scaling relationship, facilitating the initial capture of CO2 from the electrolyte solution onto the electrode surfaces could enhance the performance of CO2 electrochemical reduction, as well. Facet engineering may expose certain facets with more low-coordinated metal atoms, which are active toward CO2 adsorption and activation. Introducing defects (e.g., vacancies and dopant atoms) into the electrode surface may enhance CO2 adsorption and activation. It may also break the linear-scaling relationship of multiple intermediates binding on the same surface site. Further discussion of the role of H* in the CO2 reduction process is warranted. As proposed by Batista (subsection 2.1.1, Scheme 3C) and Goddard (subsection 3.1.4, Figure 18), H* may act as an active intermediate in CO2 reduction to formate on metal electrodes. If the targeted CO2 reduction product is formic acid, then a promising electrode material candidate should have a moderate (neither too weak nor too strong) H* binding strength, as well as a low activation energy for CO2 + H* + e− → HCOO− (at least, lower than that of H* + H+ + e− → H2). The important role of H* in the process of CO2 reduction to formate alters the traditional view of inactivating H* (by either extremely weakening or strengthening the H* binding) in order to achieve high selectivity toward CO2 reduction products. Maintaining an active H* species on the surface could lead to enhanced HCOOH production, as long as the kinetics of CO2 + H* + e− → HCOO− are faster than that of H* + H+ + e− → H2. Besides H*, the catalytic pathway through the intermediate HCOO* (with two oxygen atoms binding to the surface)260,261 may facilitate CO2 reduction to HCOOH, as well (section 3.2). Because of the weak scaling relationship between HCOO* and H*, we may find an electrode material binding HCOO*/H* moderately/unfavorably so that high selectivity toward HCOOH can be achieved260 (and HER suppressed). Finally, the local environment in the electrolyte solution near the cathode surface can influence the CO2 reduction pathways and product distributions. The local pH, which depends on mass transport rates in the electrolyte and buffering by electrolyte cations and anions, is an important factor. Under a high rate of CO2 electrochemical reduction, protons at the cathode are overconsumed, leading to a local pH increase. Selectivity toward C2 products increases at higher pH because the critical step of C−C bond formation does not involve PT. However, buffering cations, e.g., Cs+ and Rb+ which may hydrolyze easily249 ([Cs(H2O)n]+ + H2O → [Cs(H2O)n‑1OH] + H3O+), and

attract more electron density toward the surface. Thus, the catalytic intermediate 2-PyH−* would become less negatively charged and the DHP* formation (via 2-PyH−* protonation) could be less favorable, which is beneficial for catalyzing CO2 reduction. However, if the doping cation’s electronegativity is so strong that the electron density moves toward the surface too much, then the ease of HT from 2-PyH−* to CO2 might be affected, leading to a higher kinetic barrier for HT. In contrast, if the doping cation’s electronegativity is too weak, even though HT from 2-PyH−* to CO2 could be kinetically favorable, then 2PyH−* might be easily protonated to DHP* by reacting with solvated protons and losing its functionality as a hydride donor. We therefore need to make a judicious selection of the surface doping cation with an optimized electronegativity, where both the capability for HT and the stability of the 2-PyH−* (against protonation) can be maintained. Another possible approach for optimizing the catalytic intermediate in the CO2 reduction process is Py functionalization. The idea is to add electrondonating/extracting R groups to Py, e.g., −NH2 amine or −OR ethers,202 so that the electron density on the functionalized Py can be tuned. If the functionalized Py is more negatively charged, then the kinetics of HT to CO2 is expected to be faster but the stability against protonation could be lower. Similar to the idea proposed above of surface doping, we need to optimize carefully the functionalization group so that the stability of the catalytic intermediate could be improved while maintaining a good kinetics for HT to CO2. 4.3. Theoretical Insights for the Design of CO2 Electrochemical Reduction Catalytic Electrodes

Based on our discussion of CO2 electrochemical reduction mechanisms on various metal electrodes (Cu and beyond) in section 3, one major bottleneck for improving electrodes’ catalytic activity is the correlation of adsorption energies of CO2 reduction intermediates (e.g., COOH*, CO*, CHO*, etc.). Similar limits caused by the linear-scaling relationship of intermediates’ binding energies impede catalyst development for oxygen reduction/evolution reactions.279,280 Therefore, breaking the correlation of the reaction intermediates’ adsorption energies is essential for improving the performance of catalytic electrodes. Electrode-materials design beyond pure metal elements may provide us with possibilities (e.g., metal alloys, TMOs, and metal chalcogenides) to break the linearscaling relationship. A recent review article280 from the Nørskov group summarized some surface-engineering strategies to break the linear-scaling relationship, including surface multifunctionality, electrolyte ions as promoters for a specific surface intermediate, tethering/functionalization of the electrode surface, optimization of electrolyte molecules, utilization of adsorption sites at phase boundaries, and confinement by a pore or channel in electrode materials. Recently, one means to break the linear-scaling relationship on a single metal site was proposed using a metal-porphyrin-like system.281 In this theoretical study, the authors reported a linear scaling relationship between H* and COOH* (carbon bound at the single metal site), where the reduction of H*/COOH* leads to H2/CO production. They discussed that this scaling relationship between H* and COOH*presumably originating in the similarity of their single, covalent σ-bondsmakes it less likely to suppress the HER by CO production in this metal-porphyrinlike system because H* was always bound on the metal site rather than COOH*. However, no linear scaling relationship was found between the free energies of H* and OCHO* (one AC


Chemical Reviews

Review

buffering anions (HCO3− and HPO42−) in the electrolyte may act as proton sources to mitigate local pH variations. These factors interact with each other during CO2 reduction, influence the local chemical environment near the cathode surface, and thus must be tailored to optimize the selectivity of CO2 reduction products.

heterogeneous reaction step involves an electron transfer from the electrode surface slab to the adsorbate (e.g., a PCET step). Because of the surface’s limited area in the supercell, an unphysical work-function change of the electrode will be induced after the electron transfer. Under the experimental conditions, the work function (or the potential) of the electrode is fixed at a constant value during the redox reactions occurring at the electrode surface. The origin of this finite-size error is the difference between a constant charge simulation vs a constant voltage simulation. This topic and a corresponding correction strategy were discussed in subsection 1.3.5. Note that the correction scheme based on a capacitor-like model is limited to the cases where the adsorbates have small dipoles. The generality of this correction for modeling a constant potential reaction step is still in question. As for the constant potential simulation, an alternative strategy is modeling a grand canonical ensemble where the charge in the electrode slab can be varied. This approach will induce a net charge in the cell, thus various methods can be used to compensate the net charge. We summarized these methods in section 1.3.5. In principle, an ideal approach (for compensating the net charge) should perform well in the following aspects: (1) Low computational cost. We want to avoid too many explicit solvent molecules in the cell above the surface slab. This point relates to the deficiency of the double-reference method172 where many H2O molecules have to be included in the cell to represent the bulk environment of water. Another critical factor relevant to the computational cost is the size of the supercell. For the recently proposed solvated jellium model,176 four distinct layers (the surface slab, the explicit solvent layers, the implicit solvent region where the jellium counter charge is immersed and the vacuum layer) are employed, thus requiring a relatively large supercell along the dimension perpendicular to the electrode surface. This feature of the solvated jellium model may affect its overall performance in terms of computational efficiency. (2) Elimination of an artificially high electric field across the solvent region near the surface. Let us imagine that the surface slab is negatively charged (by adding more electrons in the self-consistent calculation) and a counter-charge plane is placed in the solvent region. An unphysically high electric field will exist between the two oppositely charged plates of this artificial capacitor without appropriate screening. Consequently, if adsorbate molecules or explicit solvent molecules are at the surface, the added electrons in the self-consistent calculations would prefer to localize on the adsorbates or solvent molecules instead of the surface slab. This unphysical behavior is the shortcoming of the counter electrode plane method173 due to the lack of screening within the solvent region. (3) A well-defined field-free region for referencing the electrode’s potential. The goal of a constant potential electrochemical simulation is to predict the energetics of a reaction at a given applied potential, in order to compare with experimental results. Therefore, it is essential to assign the electrodes’ absolute potential as a function of the amount of introduced electrons/ holes in the electrode slab. We could do this if there is a field-free region in the solvent space where the electrostatic potential is flat. The electrostatic potential in this field-free region corresponds to the potential of the bulk solvent. If we can establish the relationship between the bulk solvent potential with the potential of the vacuum level or SHE, the electrodes’ absolute potential at different charged states can be determined. Within the double-reference method,172 a homogeneous counter charge background fills the supercell. Hence, there is no region in the supercell where the electric field is screened

4.4. Limitations of the Current Electrochemical Models and Outlook

At the end of this discussion section, we summarize the key limitations and open questions of the current QM-based models for the simulation of electrochemical systems. The first issue that should concern most researchers is the balance between accuracy and computational efficiency. For example, the DFT-based periodic slab model is the most commonly used method for modeling heterogeneous catalytic reactions. Within the framework of DFT calculations, there are different levels of XC functionals (LDA, GGA, hybrid functionals, etc.). If we need to treat accurately an electrode surface containing open-shell TM ions, then DFT+U or hybrid functionals should be employed to mitigate SIE. But under PBCs, hybrid functionals are very computationally expensive, while DFT+U calculations have only a small computational overhead over DFT calculations. However, the choice of U depends on the TM elemental identity, oxidation state, and spin state. One might even be tempted to specify different U values for different properties being calculated. For example, it is very difficult to find a universal U value for a TMO that simultaneously can reproduce its experimental formation energy, band gap, and magnetic ground state. Overall, future exploration of heterogeneous catalytic reaction mechanisms still relies on the development of ever better XC functionals that accurately describe even strongly correlated electron systems with less demanding computational cost and good transferability for various properties and materials of interest. The second issue is the finite-size error that always exists in the simulation of electrochemical systems by a periodic slab with a finite supercell size. So far, even using a low-cost XC functional (e.g., LDA or PBE), the upper limit of the total number of atoms in an electrochemical surface model is typically around 300 (considering the computational efficiency). Since the supercell size is limited and the calculation is conducted under PBCs, one source of finite-size error is the spurious interaction of the adsorbate with its periodic images, which might be nonnegligible when the adsorbate has a relatively large dipole. This unphysical dipole interaction can be corrected in the calculations, as mentioned in subsection 1.3.1. However, if the adsorbate species contains a monopole (i.e., the supercell is charged overall), then the simulation under PBCs has to add a uniform compensating background charge to avoid divergence of the Ewald summation. The smeared-out artificial background charge causes an unphysical interaction with the adsorbate and substrate. We summarized schemes devised to minimize this error in subsection 1.3.1. Of course, getting rid of the constraint of PBCs yields the cluster model that provides an alternative way to model heterogeneous reactions on electrode surfaces, where a charged surface species will not be a problem. However, edge artifacts are a drawback of the cluster model, which is ameliorated by embedding in a potential that mimics the extended crystal. A detailed discussion regarding the cluster model is in subsection 1.3.1. Besides the interaction of the adsorbate with its periodic images, another type of finite-size error emerges if the modeled AD


Chemical Reviews

Review

As briefly mentioned in the above paragraph, AIMD simulation could provide us with an ensemble of microstates during the PCET reaction steps. However, the expense of AIMD limits the total simulation time less than ∼100 ps. If we need to model a PCET step with a relatively high barrier, then the simulation time of a standard AIMD method would not be sufficient. Accelerated sampling methods, e.g., metadynamics and constrained MD, provide the possibility to model the complete freeenergy diagram of a single or multiple PCET step(s). As discussed in subsection 1.3.6, the major drawback of these two accelerated MD methods is the choice/definition of the CV. Although the CV must represent a key reaction coordinate variation during the PCET steps, the exact mathematical expression of the CV is still arbitrary. The sensitivity of the modeling results with respect to the choice of CV could be uncertain for many electrochemical systems involving complicated PCET steps. Future theoretical development for judicious and validated CVs will provide researchers with well-founded insights into the multidimensional free-energy landscape of PCET reactions. The last limitation of current electrochemical models relates to kinetic barrier modeling approaches. Recent heterogeneous electrochemical modeling studies spend increasing effort to calculate PCET kinetic barriers, which is an improvement compared to the early CHE model. The general approach is to (1) use a TS search algorithm to identify the TS along a certain reaction pathway and (2) perform total energy calculations of the initial state (IS), TS, and final state (FS) to get the reaction activation energy. The computed results of the PCET steps’ freeenergy landscape assumes the PCET occurs adiabatically. However, PCET steps in certain systems may need to be modeled nonadiabatically if the diabatic coupling is small (

Theoretical Insights into Heterogeneous (Photo)electrochemical CO2

Recommend Documents