Layered SiC Sheets: A Potential Catalyst for Oxygen Reduction Reaction

The large-scale practical application of fuel cells cannot come true if the high-priced Pt-based electrocatalysts for oxygen reduction reaction (ORR) cannot be replaced by other efficient, low-cost, and stable electrodes. Here, based on density functional theory (DFT), we exploited the potentials of layered SiC sheets as a novel catalyst for ORR. From our DFT results, it can be predicted that layered SiC sheets exhibit excellent ORR catalytic activity without CO poisoning, while the CO poisoning is the major drawback in conventional Pt-based catalysts. Furthermore, the layered SiC sheets in alkaline media has better catalytic activity than Pt(111) surface and have potential as a metal-free catalyst for ORR in fuel cells.

We then turn to consider the adsorption and reduction of O 2 on a single-layer SiC. The adsorption energies of adsorbates on layered SiC sheets are determined by E ad 5 E ads 1 E SiC 2 E ads/SiC , where E ads , E SiC and E ads/SiC are the total energies of the isolated adsorbate molecule, the layered SiC sheets and the adsorption systems, respectively. Positive adsorption energies correspond to exothermic adsorption processes, while negative adsorption energies do endothermic ones. Several possible high-symmetry adsorption sites are considered, including the top, bridge and hollow sites. The most stable adsorption structures and energies of O 2 and ORR intermediates are summarized in Fig. 1 and Table 1. The meta-stable adsorption structures and energies are displayed in Fig. S2 of supporting information. It is found that the most energetically favorable site for O 2 adsorption on the single-layer SiC is characterized by O 2 nearly parallel to the SiC bond with adsorption energy of 0.56 eV. The O-O bond is elongated by 0.29 Å due to the electronic charge transfer from SiC to the 2p* orbital of O 2 , indicating that O 2 can be dissociated easily on SiC. The charge transfer from the single-layer SiC to O 2 molecule is quantified as 0.34 e according to Hirshfeld population charge distribution. Furthermore, the maximum coverage of O 2 on SiC is also explicitly identified. As shown in Table S2 and Fig. S3 of supporting information, only two O 2 molecules can be adsorbed on the single-layer SiC with 2 3 2 supercell, corresponding to a coverage 0.5 ML (monolayer) (1 ML is defined as one O atom per surface atom). The adsorption energy of the first O 2 molecules on single-layer SiC is 0.48 eV. The second O 2 molecule has two distinct adsorption structures with adsorption energies of 0.61 and 0.13 eV, respectively. It is found that the energetically more favorable configuration of second O 2 forms two chemical bonds with Si atoms. This special adsorption structure results in the stronger adsorption of the second O 2 molecule than the first one. This is because the electronegativity of Si (1.90) is smaller than that of C (2.55) while the Si-O bond is stronger than the C-O bond. Similar results are also found in 4 3 4 supercell, where the coverage also reaches 0.5 ML. Furthermore, O 2 adsorbed as ordered structures in a 4 3 4 supercell is considered, as shown in Fig. S3 of supporting information, the corresponding coverage can reach 0.75 ML. However, the coverage of O 2 on Pt(111) surface can reach 1 ML. All the O 2 molecules adsorb on the bridge sites. This is consistent with previous theoretical work about O 2 adsorption on Pt(111) surface 27,28 . Although the coverage of O 2 on layered SiC is smaller than that on Pt(111) surface, the remaining adsorption sites on single-layer SiC is favorable for adsorption of other reactants, which may be beneficial for ORR.
The most stable adsorption site for O atom is the bridge site, which gives rise to a highly stable epoxide-like structure with adsorption energy of 4.18 eV. In addition, H and OH prefer to adsorbing at atop site of Si atom other than that of C atom with the adsorption energies of 1.32 and 2.80 eV, respectively, which differs from the case in Ndoped carbon materials 29 30,31 .
For CO adsorption, adsorption energy value changes from 20.08 eV to 0.07 eV after considering the effect of van der Waals bonding. Since the small adsorption energy value of CO adsorption, layered SiC sheets show excellent CO tolerance. In contrast, the adsorption energy of CO on Pt(111) surface (1.86 eV) is twice as large as that of O 2 (0.84 eV), as shown in Table S3 of supporting information. Therefore, CO blocks the active sites and hinders the ORR on Pt(111) surface 32,33 . From partial density of states (PDOS) (Fig. S5 of supporting information), it is found that the hybridization between CO molecule and Si atom is very small, suggesting that the interaction between CO and SiC is little. By contrast, the interaction between CO and Pt(111) surface is very strong. All the main orbitals of CO hybridize with Pt states. The renowned mechanism of COmetal interaction, namely, via donation from CO-5s to metals and back-donation from metals to CO-2p* orbital, is applicable for CO adsorbed on Pt(111) 34,35 . This is confirmed by the fact of electron depletion at CO-5s and partially occupied CO-2p*, which is empty  above the Fermi level in gas CO. Consequently, CO will poison the Pt(111) surface. It is known that ORR mechanisms at cathodes in acidic and alkaline solutions are different. In an acidic solution, the electrode reaction can be written as O 2 Both reaction pathways are considered here. In general, ORR can proceed in Langmuir-Hinshelwood (LH) or Eley-Rideal (ER) mechanisms. LH mechanism involves all the reacting intermediates on the surface, whereas ER mechanism involves species from the electrolyte reacting with a surface intermediate. Both mechanisms are also taken into account one by one.
The ORR following the LH mechanism is discussed firstly. The ORR mechanism on the single-layer SiC in LH mechanism under acidic media can be divided into five elemental steps: O 2 dissociation, two O atom hydrogenation steps and two OH hydrogenation steps, as shown in Fig. 2. Some possible ORR pathways with higher barrier energies are shown in Fig. S6 of supporting information. The activation energies for the five elemental steps are 0.29 eV for O 2 dissociation, 0.43 and 0.49 eV for first and second OH formation, 1.11 and 1.05 eV for first and second H 2 O formation, respectively. The corresponding reaction energies are 21.61, 21.13, 21.00, 0.44 and 0.51 eV. H 2 O formation is the rate-determining step (RDS) in the O 2 dissociation mechanism of ORR in the above process. The Si-H bond broken from the initial state to the transition state contributes to the most part of the activation energy of H 2 O formation. As the number of H atoms introduced to the O atom increases, the activation energy increases from 0.49 eV to 1.05 eV and the reaction energy increases from 21 eV to 0.51 eV, corresponding to the Brønsted-Evans-Polanyi relation between activation energy and reaction energy in the heterogeneous catalysis [36][37][38] . Activation energy is a almost linear function of reaction energy, and reactions belong to the same class even follow the same relation [36][37][38] . The decrease of the activity of O atom is due to the increase of the number of introduced H atoms.
ORR on the single-layer SiC in LH mechanism under an alkaline environment can be divided into the two elemental steps: Fig. 3. Differing from that in acidic environment, both steps are easy to cross with small activation energies (0.11 and 0.23 eV) and exothermic reaction  temperature. This is consistent with the small activation energies for ORR elemental steps in alkaline media. Thus, we predict that ORR in LH mechanism on the singlelayer SiC in the alkaline media is more favorable than that in the acidic media.
As N increases from 1 to 3, adsorption energies for ORR intermediates, activation and reaction energies for ORR elemental steps in LH mechanism on layered SiC sheets are almost unchanged as shown in Tables 1 and 2. Thus, layered SiC sheets with different N values should show similar ORR catalytic activities.
For a comparison purpose, ORR elemental steps in LH mechanism on Pt(111) surface are also calculated, as depicted in Fig. S7 and Tables S4 of supporting information. The five barriers in O 2 dissociation mechanism are calculated to be 0.85 eV for O 2 dissociation, 1.00 and 1.22 eV for two OH formation steps, and 0.38 eV for H 2 O formation. For the OOH association mechanism, it is found that the barriers are 0.61 eV for OOH formation, 0.03 eV for OOH dissociation, 1.22 eV for OH formation, and 0.38 eV for H 2 O formation. OH formation is RDS for both the O 2 dissociation and OOH association mechanisms. This is consistent with the recent literature 30 . The RDS of ORR on Pt(111) surface in acidic media is the O atom hydrogenation to OH with activation energy of 1.22 eV, which is similar to that of ORR on layered SiC sheets, suggesting that layered SiC sheets and Pt(111) surface may exhibit similar activities for ORR in LH mechanism under an acidic environment. In the alkaline media, however, the activation energies for the two elemental ORR steps in LH mechanism on Pt(111) surface are 0.43 and 0.55 eV, respectively, both of them are larger than that on layered SiC sheets since the bond length of Pt-Pt (2.77 Å ) is longer than that of Si-C (1.79 Å ). The energy cost for the longer distance diffusion of H atom on Pt(111) surface is larger than that on layered SiC sheets. In addition, the elemental step of O 1 H 2 O R 2OH, which is the RDS of ORR in LH mechanism on Pt(111) surface in alkaline media, is endothermic by 0.51 eV, while that on layered SiC sheets is exothermic with reaction energy of nearly 20.40 eV, suggesting that ORR in LH mechanism on layered SiC sheets is more advantaged than that on Pt(111) surface in alkaline environment.  Although energy diagram of ORR in ER mechanism is quite different from that in LH mechanism, both of which show similar results somehow. As shown in Fig. 4(a), ORR in ER mechanism exhibits exothermic reactivity at electrode potential U 5 0 V under an acidic medium. As the pH value increases, the energy of each net coupled proton and electron transfer (CPET) step is shifted due to the effect of concentration on the free energy of H 1 . At pH 5 1, the reaction energies of five elemental steps calculated are 21.61 eV for O 2 dissociation, 21.44 and 21.32 eV for first and second OH formation, 20.06 and 20.34 eV for first and second H 2 O formation, respectively. When pH $ 4, the reaction energy of first H 2 O formation changes from exothermic to endothermic, while other steps remain exothermic. As U increases from 0 V to the ideal electrode potential of 1.23 V, energy levels are shifted for each net CPET step 31,39 . The reaction energies of the reactions in the LH mechanism and O 2 dissociation in the ER mechanism remain unchanged. This is because no CPET step involves in these steps and the electrode potential only affects the energy levels of the electrons in the CPET steps 31,39 . At U 5 1.23 V, OH formation in ER mechanism is still exothermic, while two H 2 O formation steps change to 1.17 and 0.89 eV endothermically at pH 5 1, respectively. Consistent with ORR in LH mechanism, H 2 O formation is also the RDS of ORR on single SiC layer in ER mechanism at U 5 1.23 V. The high concentration of H 2 O in electrolyte can improve the rate constant for the backward reaction and may have a negative effect on the entire process. Just as any coin has two sides, the issue being discussed here is no exception. H 2 O solution also has some positive effects on the ORR. As shown in Table S5 of supporting information, the adsorption of O 2 is enhanced in H 2 O solution, resulting in the increase of reactant concentration on catalyst surface and forward reaction rate.
A substantially different picture is obtained for ORR in ER mechanism under alkaline environment. Similar with that under acidic media, pH imports similar effect and shifts up the energy level of every CPET step. As shown in Fig. 4 ORR in ER mechanism on Pt(111) surface is also calculated. As shown in Fig. S8(a) of supporting information, all elemental steps involved in ORR via ER mechanism on Pt(111) are exothermic at pH Based on Eyring's canonical transition state theory, our calculations can be incorporated into a reduced kinetic model that should report qualitative features of the reaction mechanisms at different applied potentials 39,42,43 . As displayed in Fig. 5(a), the RDS of ORR in ER mechanism is O 2 dissociation at low potential region, while it changes to H 2 O formation when U . 0.63 V at pH 5 1. As shown in Fig. 5(b), O 2 reaction with H 2 O is the RDS of ORR in ER mechanism under alkaline environment and the reaction rate decreases as electrode potential increases. Comparing Fig. 5 with Fig. S9 of supporting information, we can predict that the catalytic activity of layered SiC sheets is better than that of Pt(111) surface in alkaline media. Furthermore, due to the higher doped ratio, SiC can provide more ORR active sites than present doped carbon materials and may act as a feasible low-cost metal-free ORR electrocatalyst.

Discussion
Since O 2 adsorption and dissociation are the key points for ORR on layered SiC sheets, the origin of the interaction between O 2 and the single-layer SiC is investigated by considering their electronic structure changes.  In summary, our DFT calculations suggest that novel metal-free layered SiC sheets with N 5 1 , 3 can exist stably and possess potential ORR catalytic activity due to two advantages of layered SiC sheets compared to Pt(111) surface: (i) free from CO poisoning, and (ii) lower activation energies for the RDS of ORR on layered SiC sheets in alkaline media than that on Pt(111) surface. In addition, the corresponding electronic structures are analyzed. The results show that the layered SiC sheets are candidates for practical applications in fuel cells.

Methods
All calculations are performed within density functional theory (DFT) framework as implemented in DMol 3 code 46,47 . The All Electron Relativistic core treat method is implemented for relativistic effects, which explicitly includes all electrons and introduces some relativistic effects into the core. The generalized gradient approximation (GGA) with Perdew2Burke2Ernzerhof (PBE) functional is employed to describe exchange and correlation effects 48 . PBE is the most common density functional in materials and surface science. However, it cannot accurate describe the van der Waals forces [49][50][51][52] . This is mainly because the GGA-PBE fails to describe non-local dispersion forces, which are expected to be relevant in layered inorganic compounds and weak adsorption systems. The lack of inclusion of long range correlation in the local density functional (LDF) calculations prevents the accurate calculation of the van der Waals bonds 46 . Therefore, a Grimme approach is adopted for dispersion corrections. Dmol 3 uses a double set of numerically tabulated basis functions 46 . A more precise term would be ''double numerical basis'' to be contrasted with double zeta basis, where the radial functions are defined as Slater zeta functions. The basis set can be significantly improved by adding higher angular momentum valence polarization functions and also by core polarization functions 46 . Total energy for the double basis is quite uniformly higher than the ones for the extended basis sets. It has shown that the basis set produces errors in self-consistent eigenvalues and total energy of <0.00003 Ha. In this work, the double numerical atomic orbital augmented by a polarization p-function (DNP) is chosen as the basis set 46 . The accuracy of DNP is comparable to a Gaussian 6-31(d) basis 53 . The double basis set may be considered as a large basis especially for the larger molecules. Recent theoretical works based this basis set have shown excellent consistency with experiments 54,55 . We have also compared our results with that obtained based on triple numerical plus polarization (TNP), as shown in Tables S6 and S7. It can be found that the adsorption energies based on TNP are slightly larger than that based on DNP, while the activation energies and reaction energies of ORR in LH mechanism based on both basis sets are almost the same. A smearing of 0.005 Ha (1 Ha 5 27.21 eV) to the orbital occupation is applied to achieve accurate electronic convergence. The spin-unrestricted method is used for all calculations. To ensure high-quality results, the real-space global orbital cutoff radius is chosen as high as 4.6 Å in the computations. The k-point density is set as 2/3 3 2/3 3 1 for per unit cell. The convergence tests for k-point density are shown in Table S8 of supporting information, which shows that the energy of single SiC layer is converging when k-point density is larger than 2/3 3 2/3 3 1. The convergence tolerance of energy is 1.0 3 10 25 Ha, maximum force is 0.002 Ha/Å , and maximum displacement is 0.005 Å in the geometry optimization. The transition states for ORR elemental steps are obtained by LST/QST tools in DMol 3 code. Frequency calculations are performed to confirm the location of the transition states. A conductor-like screening model (COSMO) is used to simulate a H 2 O solvent environment throughout the whole process. All data, except for clearly stated, are obtained under this method. The COSMO is a continuum model in which the solute molecule forms a cavity within the dielectric continuum of permittivity 56 29 . In addition, some other results also show that this implicit solvation model is an effective method to describe the solvation 30,31 . The solvation energies of ORR intermediates on SiC sheets are shown in Table S9 of supporting information and adsorption energies in a gas phase environment are shown in Table  S5 of supporting information, from which we can find that O 2 and OH are stabilized on SiC sheets by solvation. For layered SiC sheets, all simulations are performed in a 3 3 3 supercell except for coverage examination of O 2 adsorption, where 2 3 2 and 4 3 4 supercells are used. The minimal distance between the SiC layers and their mirror images is set as 15 Å , which is sufficiently large to avoid the interaction between them. The Si-C bond length is 1.79 Å . The Pt(111) surface is modeled by a periodic array of Pt slabs with a vacuum width in excess of 15 Å . The p(3 3 3) unit cells with three layer Pt slabs are used throughout in the ORR calculations. The bottom two Pt layers are fixed at their bulk positions and the top Pt layer is allowed to relax fully. Tests performed with four and five atomic layers, top two of which are allowed to relax, do not show significant differences in structural parameters. As shown in Table S10 of supporting information, the adsorption energy variation for adsorbed molecules is smaller than 0.09 eV.
In this work, both LH and ER mechanisms are expected to proceed. The LH reaction referrs to the reaction of adsorbed hydrogen atoms with another adsorbate, and the ER reaction does the reaction of a proton from solution interacting with an adsorbate. The complete electrochemical ORR in ER mechanism involves four CPETs to O 2 molecule at the cathode 31,39 . Electron transfers are coupled with a proton transfer as well. Barriers for electrochemical proton-transfer have been calculated for the reduction of O 2 to OOH and OH to H 2 O on Pt 57,58 . In both cases, the protontransfer reaction barrier calculated is small (0.15 eV to 0.25 eV) at the low potential that elementary steps are exothermic, and diminishes with higher applied voltages 57,58 . Similarly, as a first approximation, we expect also that barriers for electrochemical proton transfers to adsorbed species will be small and be easily surmountable at room temperature. As a result, we only calculated reaction energy for ORR in the ER mechanism. Free energies of the ORR intermediates are calculated based on a computational hydrogen electrode (CHE) model suggested by Nøskov et al 40,59 . Free energy change (DG) is determined by DG 5 DE 1 DZPE 2 TDS 1 DG pH 1 DG U , where DE is the reaction energy, DZPE is the zero point energy, T is temperature and DS is the change in the entropy. DG pH and DG U are the free energy contributions due to variations in H 1 concentration and electrode potential, U, respectively. In this work, we consider the contributions of DE, DG pH and DG U to free energy and neglect the effects of other terms 60 . Effects of other terms will be included in a forthcoming publication. We assume pH 5 0, 1, 3 and 5 for acidic medium and pH 5 9, 11 and 14 for alkaline medium. The pH effect is very hard to be considered directly in the electrolyte. Generally, it can be treated following the method directed by Nørskov et al 40 . At a pH differing from 0, the free energy of H 1 ions can be corrected by the concentration: G pH 5 2k B Tln[H 1 ] 5 pH 3 k B Tln10. This pHdependence effect does not enter the COSMO-approach.
Naturally, ORR is a highly complicated process. Undertaking a kinetics analysis using rate constants derived from first-principles calculations would be an ideal way to determine which pathways are relevant under different conditions. Such a kinetics analysis is unfortunately difficult, since it would require rigorous double-layer effects, such as intermediate concentrations, surface coverages, or the potential drop within the interface. At here, the potential dependent rate constants k(U) are obtained based on Eyring's canonical transition state theory: k(U)~k B T h exp {DG(U)