CO₂ activation at atomically dispersed Si sites of N-doped graphenes: Insight into distinct electron mechanisms from first-principles calculations

Catalytic conversion of CO₂ is of great significance in reducing the emission of greenhouse gases and the optimal utilization of carbon resources. Over the past few decades, many transitional metals and their compounds have been demonstrated to be highly effective catalysts for CO₂ reduction in value-added chemicals.^1,2 However, metal catalysts often suffer from high cost, poor durability, and resource scarcity as well as causing environmental pollution, which may limit their large-scale applications.^3–5 Nowadays, it is highly required to seek efficient alternatives to metal-based catalysts for CO₂ transformation.

In recent years, nonmetallic catalysts have attracted considerable interest for their high earth abundance and environmental friendliness.^6–8 Among metal-free catalysts, silicon (Si) exhibits potential activity for CO₂ catalytic conversion, both experimentally^9–11 and theoretically.^12–15 At present, scientists have obtained an in-depth understanding of the activity of transition metals based on the well-known d-band center model proposed by Hammer and Norskov.¹⁶ Most recently, the p-band center model was also proposed for the p-block element.^17–20 Intriguingly, since the p-band is much more delocalized than the d-band, they may have different bonding interactions with the adsorbate.²¹ Specifically, a deeper p-band center may correspond to either large^19,22,23 or small^13,15,24 adsorption energies. Accordingly, it is important to figure out the bonding mechanism for the p-block element toward molecular adsorption.

Previous calculations by Zhou et al. indicate that higher Si-p or Si-p_z orbitals in energy lead to stronger binding to CO₂ for Ag (111) supported silicene and silicon nanocages.^13,15 Our previous calculations also revealed that CO₂ can be easily captured and further reduced to HCOOH, CH₃OH, and CH₄ at the single Si site of the SiN₄C₄ monolayer.¹² However, the working mechanism for CO₂ activation by Si and the essence of its bonding environment effect on the catalytic performance are still not thoroughly understood. Herein, we further designed a series of Si-coordinated N-doped graphene monolayers and explored the mechanism of CO₂ activation by different SiN_x moieties. Since SiN₃C₀ and SiN₄C₀ show good CO₂ capture ability, the CO₂ electroreduction performances on both monolayer sheets were further investigated.

Here, all density functional theory (DFT) calculations were performed by the Vienna ab initio simulation package (VASP).²⁵ The Perdew–Burke–Ernzerhof generalized gradient approximation (GGA-PBE) functional with the plane-wave basis set was adopted to treat electron–electron interactions.²⁶ The van der Waals (vdW) corrections were considered by using the DFT-D3 method proposed by Grimme et al.²⁷ An energy cutoff of 500 eV was adopted, and the convergence criterion for energy and force was set to 10⁻⁵ eV and 0.01 eV/Å, respectively. The Brillouin zone was sampled by a 5 × 5 × 1 k-point grid. A vacuum space of 20 Å was set to avoid the effect between periodic images. The ab initio molecular dynamics (AIMD) simulations with the NVT ensemble were performed, and the VASPKIT code²⁸ served as a postprocessing tool for the computational data of the VASP. The transition states of CO₂ adsorption on the Si-embedded N-doped graphene sheets were identified by the climbing image nudged elastic band (CI-NEB) method.²⁹ 

The binding energy (E_b) of single Si is calculated as

where $ESiNxCy$⁠, $ENxCy$⁠, and E_Si represent the electronic energies of Si-embedded N-doped graphene, only N-doped graphene (without Si), and the Si atom referenced to the bulk Si, respectively.

The CO₂ adsorption energy (E_a) on the SiN_xC_y sheet is defined as

where E_SiNxCy@CO2,E_SiNxCy, and E_CO2 represent the electronic energies of the CO₂ adsorbed SiN_xC_y sheet, the SiN_xC_y sheet, and the isolated CO₂, respectively.

The Si-p_z center may be used as a descriptor of CO₂ adsorption strength, and it is estimated by

where D(E) is the density of state (DOS) of the Si-p_z band at a given energy E.

Based on the computational hydrogen electrode (CHE) model proposed by Nørskov et al.,³⁰ the Gibbs free energy change (⁠$Δ$G) of each elementary step is calculated by

where $Δ$E denotes the electronic energy obtained by DFT calculations, $Δ$ZPE and $Δ$S represent the zero-point energy change and the entropy change, respectively. T is the room temperature (298.15 K), and pH is set to be zero. The limiting potential (U_L) is defined as U_L = −$Δ$G_max/e⁻. Thermodynamic correction data for the intermediates are shown in Table S1. In order to simulate an aqueous environment, the implicit solvation model implemented in VASPsol³¹ was adopted.

Based on the optimized pristine graphene sheet, whose lattice constants are a = 12.32 Å, b = 12.80 Å, and c = 20 Å and α = β = γ = 90°, two types of Si-embedded N-doped graphene nanosheets have been constructed. As Fig. 1(a) shows, the two newly designed structures are denoted as SiN_xC_3−x and SiN_xC_4−x, in which x is the number of N atoms bonded to Si, and correspondingly, 3 − x and 4 − x represent the number of C atoms bonded to Si in SiN_xC_3−x and SiN_xC_4−x, respectively. Their structural differences are shown in Fig. 1(a) and Table S2. In general, the Si atom in the SiN_xC_3−x monolayer protrudes from the modified graphene surface, while the SiN_xC_4−x sheet basically maintains a planar configuration for the Si-embedded N-doped graphene sheets. Similar spatial configurations have also been reported in metal-N₄ porous carbon monolayers.^32–35 

The active-site atoms of single-atom catalysts may migrate during the catalytic reaction, which results in the agglomeration of individual atoms as the catalytic center.^36,37 To evaluate the stability of single Si atoms anchored to the N-doped graphene, its binding energies in SiN_xC_3−x and SiN_xC_4−x sheets have been calculated and are shown in Fig. 1(b). Remarkably, the predicted binding energies for all structures with various Si-embedded moieties are negative, suggesting that they generally have high thermodynamic stability. Among them, SiN₃C₀ and SiN₄C₀ have relatively weak Si anchoring interactions with binding energies of −2.20 and −2.25 eV, respectively. Furthermore, the 9 ps AIMD simulations under 300 K have been further performed for SiN₃C₀ and SiN₄C₀, and selected snapshots of SiN₃C₀ and SiN₄C₀ at 9 ps and evolution of Si–N bond lengths are shown in Fig. S1. We note that both structural skeletons survive without any chemical bond breaking during the AIMD simulation, although the structures of SiN₃C₀ and SiN₄C₀ experience partial distortion and the Si–N bond lengths fluctuate near their equilibrium configurations. Accordingly, these Si-embedded N-doped graphene sheets are of high thermal stability.

The interactions between CO₂ and SiN_xC_3−x and SiN_xC_4−x have been studied by using first-principles calculations. For SiN_xC_3−x (x = 0–3), CO₂ can be chemically adsorbed on SiN_xC_3−x monolayer sheets except for SiN₀C₃, indicating that the doping of N atoms may significantly improve the ability of SiN_xC_3−x to activate CO₂. As Fig. 2 shows, for the CO₂-chemisorbed state on SiN₁C₂–CO₂, the Si–CO₂ bond length and the bond angle of CO₂ are predicted to be 2.01 Å and 141.3°, respectively. With the increase in N-doped atoms coordinated to Si, CO₂ activation is further strengthened, and amazingly, Si is bonded to both C and O atoms in SiN₂C₁–CO₂ and SiN₃C₀–CO₂. In order to understand the intrinsic mechanism of Si–O bonding, we calculated the highest occupied crystal orbital (HOCO) and the next HOCO (HOCO − 1) at $Γ$(0,0,0) points of SiN₁C₂ and SiN₂C₁. As shown in Fig. S2, HOCO is mainly distributed on the central Si atom for both SiN₁C₂ and SiN₂C₁, which accounts for the Si–C bonding through the electron donor–acceptor interaction. For SiN₂C₁, the Si atom has notable contribution to both HOCO and HOCO − 1, and such electronic features allow the Si center to behave not only as an electron donor for the Si–C bonding but also as an electron acceptor for the Si–O bonding, leading to the $η$²-CO₂-chemisorbed configuration, as shown in Fig. 2.

For SiN_xC_4−x (x = 0–4), only SiN₃C₁ and SiN₄C₀ sheets can chemically bind CO₂, and their chemisorption configurations are shown in Fig. 2. We note that predicted Si–C bond lengths and ∠OCO bond angles in SiN₃C₁–CO₂ and SiN₄C₀–CO₂ are 2.09 Å/1.97 Å and 141.2°/136.6°, respectively. Clearly, SiN₄C₀ shows stronger ability toward CO₂ activation than SiN₃C₁. In general, the more N atoms coordinate with Si in both SiN_xC_3−x and SiN_xC_4−x sheets, the more remarkable the CO₂ activation is.

The predicted transition-state (TS) configurations for formation of CO₂-chemisorbed states on SiN₁C_2, SiN₂C₁, SiN₃C_0, SiN₃C₁, and SiN₄C₀ monolayer sheets are depicted in Fig. S3, and the corresponding energy barriers relative to their physical adsorption configurations are calculated to be 0.09, 0.12, 0.23, 0.22, and 0.05 eV, respectively, revealing the facile CO₂ capture and chemisorption on these five single-atom Si-embedded N-doped graphene sheets. Herein, based on the unit cell of SiN₄C₀, a 2 × 2 supercell with four atomically dispersed Si sites has been constructed. Then, 16 CO₂ molecules were randomly put into the supercell, and 5 ps AIMD simulations were further performed. Figure S4 shows the final configuration after 5 ps simulations and the evolution of the Si–CO₂ distance over time, where one CO₂ molecule is chemically adsorbed at the Si site at about 4 ps at room temperature, showing that CO₂ is easily captured and activated by the SiN₄C₀ sheet.

It was known that the bending of CO₂ may lower the LUMO energy level,³⁸ which is beneficial for accommodation of excess electrons. Accordingly, the injection of electrons into the LUMO orbital of CO₂ could activate and bend CO₂. Figure S5 depicts the electron spin density of graphene, SiN_xC_3−x, and SiN_xC_4−x sheets. One can see that the introduction of Si and N atoms makes the electron spin density localize over different atoms. However, except that SiN₁C₂ has a net magnetic moment of 0.96 µ_B, the overall magnetic moment of other nine structures is zero. It is worth noting that the spin density is partially distributed on the Si atom for all sheets here, which is favorable for the bonding between Si and CO₂.

The valence electron (VE) populations of the Si atom and the number of electrons occupied by the Si-p_z orbital in SiN_xC_3−x and SiN_xC_4−x are compiled into Table I. One can see that as the number of N atoms increases, the valence electrons of the Si atom in both SiN_xC_3−x and SiN_xC_4−x increase first and then decrease. However, the electron population occupied by the Si-p_z orbital always increases for SiN_xC_3−x, while for SiN_xC_4−x, it increases at first and then decreases. Such discrepancies may result from the strong conjugation interaction between Si-p_z orbitals and their surrounding atoms for SiN_xC_4−x with a planar configuration. Here, charge transfer plays an important role in CO₂ activation. In the process of CO₂ chemisorption on SiN₁C₂, SiN₁C₂, SiN₁C₃, SiN₃C₁, and SiN₄C₀, the number of gain and loss electrons of Si and CO₂ is presented in Table I. We note that for SiN_xC_3−x, the electrons acquired by CO₂ are mainly from the Si atom, while for SiN_xC_4−x, the Si atom almost has no electron loss, and the electron transfer to CO₂ basically comes from the N-doped graphene substrate. In brief, the Si atom donates electrons to CO₂ for SiN_xC_3−x, while it acts as an electron shuttle to transfer electrons to CO₂ from the SiN_xC_4−x framework. Overall, the net effect results in the electron injection into CO₂ and in its activation and chemisorption.

Herein, the Si-p_z centers [$ε$(p_z)] as well as the adsorption energies of CO₂ were also calculated and are shown in Table I. We note that the more the doping-N atoms, the lower the Si-p_z center energy level. For SiN_xC_4−x, the low-energy Si-p_z center corresponds to relatively strong chemisorption of CO₂, while for SiN_xC_3−x, the chemisorption strength of CO₂ depends on the balance between the Si-p_z center energy level and the electron population in the Si-p_z orbital and Si atom, as shown in Table I. Here, the Si center in SiN_xC_3−x acts as an electron donor in CO₂ activation, and a high electron density at the Si site and a high-energy Si-p_z orbital benefit the electron transfer to CO₂, whereas the Si atom in SiN_xC_4−x serves as the electron shuttle to deliver the electron transfer from the SiN_xC_4−x substrate to CO₂. The lower the Si-p_z center energy level, the more favorable the electron capture is. In the transition state for CO₂ chemisorption on SiN₄C₀, the single-atom Si has a valence electron population of 1.83, higher than that of the pristine SiN₄C₀ by 0.47 e⁻. Presumably, with the increase in electrons in the Si-p_z orbital, its energy level increases, and the excess electrons may transfer to CO₂. As shown in Fig. 2, SiN_xC_4−x sheets maintain an approximately planar configuration, and there should be strong conjugation interactions between Si and the SiN_xC_4−x framework, which are beneficial to the electron transfer.

In order to understand the intrinsic relationship between the Si-p_z center and the CO₂ adsorption energy, a detailed density of states analysis was further carried out. Figure 3(a) shows the projected density of states (PDOS) of Si-p_z and N-p_z of SiN_xC_3−x. As for SiN₀C₃, one can see that the Si-embedded N-doped graphene opens the bandgap of pristine graphene and the Si-p_z orbital has a notable share of electronic state distributions in the valence and conduction bands near the Fermi level. As the number of N atoms increases, the conduction band around the Fermi level, contributed by the Si-p_z orbital, is gradually occupied, as shown in the PDOS of SiN₁C₂. For SiN₂C₁, the conduction band with a big share of the Si-p_z orbital is almost completely occupied near the Fermi level, and correspondingly, SiN₂C₁ possesses the strongest electron-donating ability and thus shows the strongest binding to CO₂. As for SiN₃C₀, since the energy band contributed by the Si-p_z orbital is already fully occupied and its energy level is further lowered, as shown in Fig. 3(a) and Table I, the extra electrons are populated into the conduction band contributed by the N-p_z orbital, which accounts for the sharp change in the Si-p_z center of SiN₃C₀ and relatively weaker binding ability toward CO₂ than toward SiN₂C₁.

Figure 3(b) shows the PDOS of Si-p_z and N-p_z orbitals of SiN_xC_4−x. Apparently, there is a strong conjugate interaction between Si-p_z and N-p_z orbitals. Similar to SiN_xC_3−x, as the number of N atoms increases, the Si-p_z unoccupied electronic states in SiN_xC_4−x are gradually occupied. It is worth noting that the Si-p_z orbital of SiN₄C₀ still has electron state distributions in the conduction band near the Fermi level, indicating that more electrons can be accommodated. On the whole, the shift of the Si-p_z center to the low-energy region results from its population of more electrons.

Since the SiN₃C₀ and SiN₄C₀ sheets show good CO₂ activation ability, the electroreduction of CO₂ (CO₂ER) on SiN₃C₀ and SiN₄C₀ was further investigated by the first-principles calculations. Figure S6 depicts the proposed CO₂ER pathways on SiN₃C₀ and SiN₄C₀ sheets. For the production of CO, the pathway is denoted as * + CO₂ → *CO₂ → *COOH → *CO → * + CO, which is the same for both SiN₃C₀ and SiN₄C₀. However, for the production of HCOOH, CH₃OH, and CH₄, there exist multiple pathways.

For both SiN₃C₀ and SiN₄C₀, the *OCHO intermediate may be involved in CO₂ER. However, the predicted free energy change (⁠$Δ$G) of the *OCHO → *HCOOH step on SiN₄C₀ is 2.05 eV, and thus, the CO₂ER pathway including *OCHO as a precursor is infeasible. Figure S7 depicts the optimal pathway for CO₂ER on SiN₄C₀. As shown in Fig. S7, the hydrogenation of *COOH determines the speed of the entire CO₂ER. Unfortunately, the limiting potentials (U_L) are calculated to be −1.23 V for the production of CO and −1.15 V for generating HCOOH, CH₃OH, and CH₄, and such relatively high limiting potentials arise from the high stability of *COOH on SiN₄C₀. Accordingly, the SiN₄C₀ sheet may not be a kind of high-efficient electrode material, although it has an excellent ability toward CO₂ activation.

Figure 4 shows the predicted optimal pathways for the production of CO, HCOOH, CH₃OH, and CH₄ on SiN₃C₀. For CO₂ER on SiN₃C₀, the rate of the whole reaction is determined by the first three protonation steps, and the limiting potentials for the production of CO, HCOOH, and CH₃OH/CH₄ are calculated to be −0.95, −0.54, and −0.59 V, respectively; the corresponding rate-determining steps are *COOH → CO, *OCHO → *HCOOH, and *HCOOH → *CHO. Remarkably, the limiting potentials for products CO, HCOOH, CH₃OH, and CH₄ on SiN₃C₀ are comparable to those for CO₂ER on Cu–C₃N₄, and such low limiting potentials were proved to have higher efficiency than Cu–NC and Cu(111) surfaces in selective electroreduction of CO₂ to these C1 products.³⁹ Intriguingly, the Si atom is bonded with both C and O of CO₂ in the intermediates *HCOOH and *CH₂O, as also observed in the configuration of CO₂ adsorbed SiN₃C₀. Such bonding features can enhance the interaction between the Si atom and intermediate species and further produce more value-added chemicals, compared to the catalytic conversion by the metal-N_x porous carbon.^35,40 In general, SiN₃C₀ has low limiting potentials for CO₂ electroreduction and promises the selectivity for the formation of CO, HCOOH, and CH₃OH/CH₄ to some extent.

In the present work, newly designed Si-embedded N-doped graphenes are predicted to have high thermodynamic stability by the first-principles calculations. The increase in N atoms in coordination with the Si atom can enhance the electron-donating ability of SiN_xC_3−x and SiN_xC_4−x, and CO₂ chemisorption can be achieved through the Si–CO₂ bonding interaction. In chemisorption configurations of CO₂ on SiN₁C_2, SiN₂C₁, SiN₃C_0, SiN₃C₁, and SiN₄C₀, the three-coordinated Si atoms in SiN_xC_3−x act as the electron donor, while the four-coordinated Si atoms in SiN_xC_4−x only behave as the electron shuttle to account for the electron transfer from the support framework to CO₂. Accordingly, CO₂ activation and chemisorption on SiN_xC_3−x and SiN_xC_4−x follow different electron mechanisms. For SiN_xC_3−x, the activity of a single Si site toward CO₂ activation depends on its electron density and the Si-p_z band center. The more the N-doping atoms, the larger the electron population at the Si site, while the deeper the Si-p_z band center, and the CO₂ chemisorption requires a balance among these factors. For SiN_xC_4−x, the low-energy Si-p_z center is a prerequisite for the Si site to capture the electron from the support framework, which is beneficial to the electron transfer to CO₂. The shift of Si-p_z band centers originate from the change in the electron occupation of the Si-p_z-based conduction band near the Fermi level. With the increase in N-doping atoms in coordination with Si, the electron population on the Si-p_z-based conduction band increases, and the Si-p_z band center gradually moves to the low-energy region. CO₂ electroreduction to HCOOH and CH₃OH/CH₄ on SiN₃C₀ has low limiting potentials of −0.54 and −0.59 V, respectively, indicating that SiN₃C₀ as a novel electrocatalyst can be applied for CO₂ capture and electroreduction.

See the supplementary material for the evolution of the Si–N distance and CO₂ capture during AIMD simulations, the highest occupied crystal orbital, the transition-state configurations of CO₂ chemisorption, the electron spin densities, and free energy profiles for CO₂ electroreduction on SiN₃C₀.

This work was supported by the National Science Foundation of China (Grant Nos. 21873078, 21673185, and 21933009).

The authors have no conflicts to disclose.

The data that support the findings of this study are available within the article and its supplementary material.