Insights into the activity of single-atom Fe-N-C catalysts for oxygen reduction reaction

Single-atom Fe-N-C catalysts has attracted widespread attentions in the oxygen reduction reaction (ORR). However, the origin of ORR activity on Fe-N-C catalysts is still unclear, which hinder the further improvement of Fe-N-C catalysts. Herein, we provide a model to understand the ORR activity of Fe-N4 site from the spatial structure and energy level of the frontier orbitals by density functional theory calculations. Taking the regulation of divacancy defects on Fe-N4 site ORR activity as examples, we demonstrate that the hybridization between Fe 3dz2, 3dyz (3dxz) and O2 π* orbitals is the origin of Fe-N4 ORR activity. We found that the Fe–O bond length, the d-band center gap of spin states, the magnetic moment of Fe site and *O2 as descriptors can accurately predict the ORR activity of Fe-N4 site. Furthermore, these descriptors and ORR activity of Fe-N4 site are mainly distributed in two regions with obvious difference, which greatly relate to the height of Fe 3d projected orbital in the Z direction. This work provides a new insight into the ORR activity of single-atom M-N-C catalysts.

S ingle-atom Fe-N-C catalysts are considered a candidate to replace the platinum group metals in oxygen reduction reaction (ORR), due to its high activity, anti-toxicity, and metal atom utilization [1][2][3][4] . Advanced Fe-N-C materials show comparable ORR performance to that of benchmark Pt/C 3,5-7 . Both theoretical and experimental studies indicate that the Fe-N 4 site is the active species of Fe-N-C catalysts 3,4,6,[8][9][10] . However, the ORR activity of Fe-N 4 site on different carbon supports varies greatly, and Fe-N 4 ORR activity on curved carbon supports is usually orientationdependence [11][12][13][14][15] . The essential of the difference and orientationdependence of Fe-N 4 ORR activity remains unclear. Therefore, understanding the ORR activity of the Fe-N 4 site is crucial to further improve the ORR performance of Fe-N-C catalysts.
The intrinsic ORR activity of the Fe-N 4 site relates to the electronic structure of carbon supports 7,16 . Controlling the electronic structure of carbon supports alters the interaction between the supports and the Fe-N 4 site, resulting in the activity change of the Fe-N 4 site 6,7,13,[16][17][18] . However, the essential of the change in Fe-N 4 ORR activity is still puzzling due to the complex factors affecting it. Some works reported that the electronic structure of carbon supports (electron-withdrawing/donating property) can affect the d-band center or d-orbital level of Fe site 17,[19][20][21] . While, others claimed that it regulates the charge state of Fe site 6,22,23 . To understand the Fe-N 4 ORR activity, many efforts have been made to explore the descriptors of Fe-N 4 ORR activity 13,16,20,[24][25][26] . For example, the full width at half-maximum of C 1s photoemission spectra 13 , d-band center 20 , and the electronegativity of the nearest neighbor atoms 24 , are considered to describe the ORR activity of the Fe-N 4 site. However, these descriptors are difficult to understand the orientation-dependence and differences of ORR activity at the same Fe-N 4 site on different carbon supports 11,12,14 . Therefore, elucidating the origin of ORR activity and exploring the universal descriptors of ORR activity on the Fe-N 4 site are urgent.
In this work, we demonstrate that the hybridization among Fe 3dz 2 , 3dyz (3dxz), and O 2 π* orbital is the origin of ORR activity on the Fe-N 4 site by density functional theory (DFT) calculations. We established the relationship between the geometric, electronic structure, and ORR activity of the Fe-N 4 site. From the geometric structure, the Fe-O bond length (L Fe-O ) can accurately describe the ORR activity of the Fe-N 4 site. From the electronic structure, the d-band center gap of spin states (Δd), the magnetic moment of the Fe site (M Fe ), and *O 2 (M *O2 ) are nearly linear in relationship to the ORR activity of the Fe-N 4 site. We found that these descriptors and Fe-N 4 ORR activity are located in two regions. We provided a model to elucidate the origin and correlation of these descriptors from the spatial structure and level of the frontier orbital. This model reveals that descriptors essentially reflect hybridization between Fe 3dz 2 , 3dyz (3dxz), and O 2 π* orbital, which is related to the height of 3d orbital spatial distribution. This study offers new insights into the origin and descriptor of ORR activity on M-N-C catalysts and provides guidance for the design of high-performance ORR catalysts.
To evaluate the stability of single-atom Fe sites on carbon supports, we calculated the binding energy (E bind ) of single-atom Fe sites (Fig. 1c). For the perfect Fe-N 4 configuration, the binding energy of the single-atom Fe site is -1.01 eV, indicating that the configuration of the Fe-N 4 site on carbon supports is stable. Under the control of divacancy defects, the binding energies of Fe sites range from -1.88~0.63 eV. Obviously, divacancy defects on the carbon supports can affect the stability of a single-atom Fe site. Among them, the relatively stable systems (E bind <-1.01 eV) were selected as examples to study the ORR activity of Fe-N 4 sites, which include 5-8-5 (6,8,9), 555-777 (5a, 6a, 7b, 8a, 9a), and 5-7-7-5 defect (6I, 7I, 8I, 8II, 8III, 9I, 10I). The ORR activity on the upper (blue) and lower (red) sides of Fe-N 4 configurations were investigated due to surface asymmetry (Fig. 1d). We found the hydrogenation of *O 2 is a potential-determining step (PDS) at all Fe-N 4 configurations. Divacancy defect can regulate the free energy difference of this step on the Fe-N 4 site (Supplementary Table 1). The maximal free energy difference (ΔG max ) of ORR at the perfect FeN 4 site is -0.45 eV, which is consistent with previous studies 3,6,30 . Among them, divacancy defects, such as 6, 5a, 6a, 6I, 7I, and 8I, can enhance the activity of the Fe-N 4 site, with ΔG max less than -0.60 eV (Fig. 1d). For the asymmetric surfaces, the ORR activity of the Fe-N 4 site showed orientationdependence and difference, such as 6I, 7I, 8I, etc. Besides, divacancy defects can control the electron-donating capability of the carbon support, ultimately affecting the charge state and the d-band center of the Fe site ( Supplementary Fig. 4). However, these factors are difficult to elucidate the orientation-dependence of the ORR activity of the Fe-N 4 site on an asymmetric surface.
Factors of ORR activity at the Fe-N 4 site. We investigated the ORR process of Fe-N 4 configurations from energy, geometry, and electronic structure to explore the origin of Fe-N 4 ORR activity. Here, we highlight the results of 6 (5-8-5), 6a (555-777), and 6I (5-7-7-5) models because they can improve the ORR activity of the Fe-N 4 site with good thermodynamic stability (Supplementary Fig. 5). The free energy diagram of ORR shows that the ΔG value of *OOH formation on 6, 6a, and 6I are -0.62, -0.65, and -0.66 eV ( Fig. 2a and Supplementary Fig. 6), respectively. Compared with the perfect FeN 4 site (-0.45 eV), divacancy defects (6, 6a, and 6I) can facilitate *O 2 conversion to *OOH intermediate. We found that the adsorption energy of *O 2 with end-on1 and end-on2 configurations on the FeN 4 site are -0.53 and -0.41 eV, respectively, and the adsorption energy of *O 2 with side-on configuration on the FeN 4 site is -0.01 eV ( Supplementary Fig. 7). This result shows that the end-on adsorption is more favorable than the side-on adsorption at the FeN 4 site. Among them, the adsorption energies of *O 2 on 6, 6a, and 6I are -0.34, -0.33, and -0.38 eV (Fig. 2b), respectively, which are larger than that of the perfect FeN 4 site (-0.53 eV). These results suggest the weak adsorption of *O 2 on the Fe-N 4 site is the key to enhance the ORR activity of the Fe-N 4 site. Meanwhile, the O-O bond lengths on the O 2 @FeN 4 , O 2 @6, O 2 @6a, and O 2 @6I are 1.284, 1.297, 1.294, and 1.296 Å (Fig. 2b), respectively, which are longer than the bond length of O 2 molecule (1.234 Å). Besides, the Fe-O bond lengths on the O 2 @FeN 4 , O 2 @6, O 2 @6a, and O 2 @6I are 1.974, 1.745, 1.724, and 1.731 Å (Fig. 2b), respectively. These results reveal that divacancy defects (6, 6a, and 6I) can enhance the interaction between *O 2 and Fe site and promote the activation of *O 2 , thus increasing the Fe-N 4 ORR activity.
We studied the electronic structure of the Fe-N 4 site and O 2 to understand the activation mechanism of *O 2 . Before the adsorption of O 2 , the projected density of state (PDOS) of the Fe site indicates that the Fe site has spin polarization with a spin magnetic moment of 2 μB. The occupied states of the Fe site near the Fermi level consist of localized 3dxz, 3dx 2 -y 2 , and 3dz 2 ( Supplementary Fig. 8). Compared to the PDOS of the FeN 4 site, divacancy defects (6, 6a, and 6I) modified the intensity and distribution of Fe 3d states (especially 3dxz and 3dz 2 ). For the O 2 molecule, the PDOS of O 2 also shows spin polarization with a spin magnetic moment of 2 μB. The anti-bonding orbitals (π*) of O 2 are mainly composed of the O 2 2py and 2pz orbitals ( Supplementary Fig. 8). When O 2 adsorb on the Fe-N 4 site, the 2pz orbital in the O 2 π* orbital split into the discrete level, and the  intensity of localized Fe 3dz 2 states is significantly reduced (Supplementary Fig. 9). A new hybrid state formed near the Fermi level, due to the hybridization between O 2 2pz and Fe 3dz 2 orbital. In the O 2 @6, O 2 @6a, and O 2 @6I, the strong hybridization between Fe 3dz 2 , and O 2pz orbital lead to the significant discrete of the O 2 2pz orbital state (Fig. 2c), suggesting that the O 2 molecules are better activated. Besides, there is also hybridization between the 2py orbital of O 2 and the 3dyz orbital of Fe ( Supplementary Fig. 9). The d-band center of Fe spin-up and spin-down are close to each other due to hybridization ( Fig. 2c and Supplementary Fig. 4). The d-band center gap of Fe spin state (Δd) of O 2 @FeN 4 , O 2 @6, O 2 @6a, and O 2 @6I are 2.45, 1.76, 1.67, and 1.70 eV (Fig. 2c), respectively, meaning the Δd is related to the activation of O 2 .
Generally, the filling of electrons in the O 2 π* orbital is the beginning of the activation of O 2 . Charge density difference and Bader charge analysis reveal *O 2 gain 0.34, 0.38, 0.40, and 0.39 e from perfect FeN 4 , 6, 6a, and 6I (Fig. 2e), respectively. The electron-gaining region (yellow region) of the *O 2 mainly distributes in the π* orbitals of O 2 , which is consistent with the distribution region of the Wannier function ( Fig. 2d and Supplementary Figs. 13, 14), reflecting that the electron transfer between Fe 3dz 2 (3dyz) and O 2 π* orbitals. Furthermore, due to electrons transfer between Fe sites and O 2 , the Fe 3dz 2 (3dyz) and O 2 π* orbitals are broken, thus changing the magnetic moment of the Fe site and O 2 . Spin density analysis show divacancy defects (6, 6a, and 6I) hardly change the spin density of Fe sites before O 2 adsorption ( Supplementary Fig. 15). After O 2 adsorption, the spin density is mainly localized at the Fe site and *O 2 has the opposite direction (Fig. 2f). The spin density of the Fe site and *O 2 in the O 2 @6, O 2 @6a, and O 2 @6I is weaker than that in O 2 @FeN 4 (Fig. 2f). Therefore, the origin of O 2 activation is the hybridization of the Fe 3dz 2 (3dyz) and *O 2 2pz (2py) orbitals, which results in electron transfer from the Fe site to *O 2 , breaking the Fe 3dz 2 (3dyz) and O 2 2pz (2py) orbitals and reducing Δd and magnetic moments.
Descriptors of ORR activity at Fe-N 4 sites. As mentioned above, the adsorption energy of *O 2 , the bond length, the charge state, Δd, and the magnetic moment may be related to the ORR activity of the Fe-N 4 site (Supplementary Table 2). To gain accurate descriptors for the ORR activity of the Fe-N 4 site, we investigated the correlation between these factors and ORR activity ( Fig. 3 and Supplementary Figs. 16-18). According to the free energy diagram, the adsorption energy of *O 2 may be related to ORR activity. The adsorption energy of *O 2 at the Fe-N 4 site ranges from -0.6 to -0.1 eV under the regulation of divacancy defects ( Supplementary Fig. 16). The linear fitting results reveal that the coefficient of determination R-square between the adsorption energy of *O 2 and the ORR activity of the Fe-N 4 site is 0. 15 Fig. 17). When the O-O bond length is larger than 1.290 Å, and the free energy difference for hydrogenation of *O 2 is less than -0.6 eV. The coefficient of determination R-square between the O-O bond length and Fe-N 4 ORR activity is 0.81, with the relevant formula being ΔG = -9.93L O-O + 12.3 ( Supplementary Fig. 17). Besides, the length of the Fe-O bond is distributed in stages (Fig. 3a).  Supplementary Fig. 18), respectively, due to the regulation of divacancy defect. The linear fitting results show that the coefficient of determination R-square between the charge state of the Fe site (*O 2 ) and the ORR activity of the Fe-N 4 site is 0.10 (0.58), indicating that they are difficult to describe the ORR activity of Fe-N 4 site ( Supplementary Fig. 18). In addition, Δd (Fig. 3b), M Fe (Fig. 3c), and M *O2 (Fig. 3d) are also characteristic symbols for the change of Fe and *O 2 electronic structure. We found that Δd, M Fe , and M *O2 are also distributed in stages (Fig. 3b-d) To generalize the correlations between electronic structure and Fe-N 4 ORR activity, we calculated the electronic structure of O 2 @Fe-N 4 by using other theoretical approximations, including local density approximation (LDA (CA)) 32,33 , Perdew-Wang 91 (GGA (PW91)) 34 , strongly constrained and appropriately normed semilocal density functional (meta-GGA (SCAN)) 35 , and Heyd-Scuseria-Ernzerhof screened hybrid density functional (HSE06) 36 Supplementary Fig. 20), respectively. In addition, M *O2 and Δd were also highly correlated with ΔG ( Supplementary Figs. 21  and 22). Therefore, the correlation between the electronic structure and the Fe-N 4 ORR activity remains valid in other theoretical approximations.
Correlations between the descriptors of ORR activity at Fe-N 4 sites. We explored the correlation between these descriptors. Obviously, there is a nearly linear relationship between Δd, M Fe , M *O2 , and the Fe-O bond length, with the R-squares of 0.98, 0.99, and 0.94 ( Fig. 4a and Supplementary Fig. 23), respectively, indicating that these descriptors can reflect the origin of Fe-N 4 ORR activity. To reveal the essence of descriptors, we analyzed the projected magnetic moment of Fe, the spatial structure of Fe 3d and O 2 π* orbital, and the schematic frontier orbital level diagram.
The projected magnetic moment of Fe fluctuates obviously, due to the regulation of the divacancy defects (Fig. 4b). Only when the magnetic moments of 3dz 2 and 3dyz + 3dxz decrease at the same time, the Fe-N 4 site shows higher ORR activity. This result suggests that the ORR activity at the active Fe-N 4 site is attributed to the strong hybridization of Fe 3dz 2 -*O 2 2pz and Fe 3dyz-*O 2 2py (Fe 3dxz-*O 2 2px).
The spatial distribution of Fe 3d and O 2 π* orbital can elucidate the correlation between the Fe-O bond length and M Fe . As shown in Fig. 4c, the 3dxy (3dx 2 -y 2 ), 3dyz (3dxz) and 3dz 2 are different heights in the Z direction. When the Fe-O bond is long, the activation of O 2 attribute to the hybridization between Fe 3dz 2 and O 2 2pz orbital, which leads to the decrease of the projected magnetic moments of Fe 3dz 2 . When the Fe-O bond is short, the activation of O 2 is additionally improved by the hybridization of Fe 3dyz (3dxz) and 2py (2px) orbitals. Meanwhile, the projected magnetic moments of Fe 3dz 2 and 3dyz (3dxz) are significantly reduced. The Fe 3dxy and 3dx 2 -y 2 have weak hybridization with O 2 π* orbital due to their low spatial distribution in the Z direction. Thereby the projected magnetic moments of 3dxy and 3dx 2 -y 2 vary only slightly.
The schematic frontier orbital diagram of O 2 @Fe-N 4 can reveal the differences in O 2 activation (Fig. 4d) 38 . Before the adsorption of O 2 , the Fe 3d orbital splits into an empty 3dxy, singly occupied 3dyz, 3dz 2 type orbitals, and occupied 3dxz, 3dx 2 -y 2 orbitals, due to Fe 3d orbital hybridization with the carbon support. For the O 2 @Fe-N 4 , the differences in O 2 activation is attributed to the degree of hybridization of Fe 3dz 2 and 3dyz (or 3dxz) with O 2 π* orbital. For the active Fe-N 4 site, the strong hybridization of 3dz 2 and 3dyz (or 3dxz) with O 2 π* orbital promote the O 2 activation. For the inactive Fe-N 4 site, the strong hybridization of 3dz 2 with O 2 π* orbital and negligible hybridization of 3dyz (or 3dxz) with O 2 π* orbital lead to the poor O 2 activation. The change of Δd, M Fe , and M *O2 are also due to the above factors. Therefore, the Fe-N 4 ORR activity is origin from the activation of O 2 , that is, the hybridization between Fe 3dz 2 (3dyz or 3dxz) and *O 2 π* orbital. The Fe-O bond length, Δd, M Fe , and M *O2 are the specific symbols describing hybridization.

Discussion
In this work, we investigated the origin and descriptors of ORR activity on the Fe-N 4 site using DFT calculation. We explored divacancy defects on carbon supports that could modulate the ORR activity of the Fe-N 4 site. The first electron step (*O 2 + H + + e → *OOH) is the potential-determining step (PDS) on all Fe-N 4 configurations. Fe-N 4 site on 6, 5a, 6a, 6I, 7I, and 8I showed high ORR activity with ΔG max less than -0.60 eV. The ORR activity of Fe-N 4 sites on the asymmetric surface showed orientation-dependence. We demonstrated that the hybridization between Fe 3dz 2 , 3dyz (3dxz), and O 2 π* orbitals is the origin of Fe-N 4 ORR activity by studying the electronic structure. The Fe-O bond length, Δd, M Fe , and M *O2 show staged differences, which can be used as descriptors to accurately describe the ORR activity of the  177 μB). This work offers insights into the Fe-N 4 ORR activity and provides guidelines for designing the high-performance Fe-N-C catalysts.

Methods
DFT calculations. Our spin-polarized density functional theory simulation was calculated by using the Vienna ab initio simulation package (VASP) 39 . The PAW potentials and the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) were employed to describe the interaction of electron-ion and the electron-electron exchange and correlation functional, respectively 40,41 . After the test, the cutoff energy was set at 400 eV, and Monkhorst-Pack mesh with a 2 × 2 × 1 grid was used in our calculations after testing ( Supplementary Fig. 24). van der Waals (VDW) forces were corrected with the D2 method of Grimme 42 . The convergence criterion was set at 0.02 eV Å −1 for the force and 1 × 10 −5 eV per atom for energy. We used the correlation energy (U) of 4 eV and the exchange energy (J) of 1 eV for Fe 3d orbitals 43,44 . We built a super-cell (7 × 8) carbon substrate as a model, which includes 106 C atoms, 4 N atoms, and 1 Fe atom. The lattice parameters of this slab are a = 17.31 Å and b = 16.98 Å after optimization. The lattice constants are optimized again due to the introduction of divacancy defects in the system. The vacuum layer was set at 15 Å. Maximally localized Wannier functions (MLWFs) methods are applied by Wan-nier90 to calculate the spatial distribution of Fe 3d and O 2 2p orbital 31,[45][46][47] . We also considered the effect of solvent on the ORR activity of the Fe-N 4 site by employing an implicit solvent model 48,49 . To confirm the generalization validity of the correlation between electronic structure and Fe-N 4 ORR activity, local density approximation with the Ceperly-Alder functional (LDA (CA)) 32,33 , Perdew-Wang 91 (GGA (PW91)) 34 , strongly constrained and appropriately normed semilocal density functional (meta-GGA (SCAN)) 35 , and Heyd-Scuseria-Ernzerhof screened hybrid density functional (HSE06) were used to evaluate the electronic structure 36,37 . In the framework of LDA (CA) and GGA (PW91), we also consider the Coulomb interaction between Fe 3d orbital electrons, and the Hubbard U value is 3 eV. First-principles Born-Oppenheimer molecular dynamics (BOMD) simulation was performed in the canonical ensemble with Nosé-Hoover heat bath schemes 50,51 . During the simulation, spin polarization is not considered, and the Brillouin-zone integrations were performed using the Gamma-point grid. The total simulation time is 5 ps with a time step of 1 fs at 300 K. We also used the VESTA package to visualize the structures and charge density differences 52 .
The binding energy of the Fe site (E bind ) on the carbon supports can be calculated by where E Fe@sub is the total energy of the Fe atom embedded in the carbon supports; E sub and E Fe are the energies of the substrate and single Fe atom in the Fe bulk, respectively. The Gibbs free energy can be expressed as where ΔE, ΔZPE, and ΔS are the reaction energy calculated by the DFT methods, the changes in zero-point energies, and the entropy during the reaction, respectively [53][54][55] . T is the temperature (298.15 K, in our work). The d-band center gap of spin state (Δd) can be defined as where ε up is the d-band center of the Fe 3d spin-up projected density of states. And ε dn is the d-band center of the Fe 3d spin-down projected density of states 56 . The magnetic moment (M) is defined as where N(spin-up) and N(spin-down) represent the number of electrons in the spin-up and spin-down occupied state, respectively 57,58 . M(3dxy), M(3dx 2y 2 ), M(3dyz), M(3dxz), and M(3dz 2 ) represent the projected of magnetic moments on Fe 3dxy, 3dx 2 -y 2 , 3dyz, 3dxz, and 3dz 2 , respectively.

Data availability
The data that support the findings of this study are available from the corresponding author on reasonable request.

Code availability
The computational codes used in this work are available from the corresponding author on reasonable request.