Activation of Molecular O2 on CoFe2 O4 (001) Surfaces: An Embedded Cluster Study.

Dioxygen activation pathways on the (001) surfaces of cobalt ferrite, CoFe2 O4 , were investigated computationally using density functional theory and the hybrid Perdew-Burke-Ernzerhof exchange-correlation functional (PBE0) within the periodic electrostatic embedded cluster model. We considered two terminations: the A-layer exposing Fe2+ and Co2+ metal sites in tetrahedral and octahedral positions, respectively, and the B-layer exposing octahedrally coordinated Co3+ . On the A-layer, molecular oxygen is chemisorbed as a superoxide on the Fe monocenter or bridging a Fe-Co cation pair, whereas on the B-layer it is adsorbed at the most stable anionic vacancy. Activation is promoted by transfer of electrons provided by the d metal centers onto the adsorbed oxygen. The subsequent dissociation of dioxygen into monoatomic species and surface reoxidation have been identified as the most critical steps that may limit the rate of the oxidation processes. Of the reactive metal-O species, [FeIII -O]2+ is thermodynamically most stable, while the oxygen of the Co-O species may easily migrate across the A-layer with barriers smaller than the associative desorption.

Introduction

Dioxygen, O2, is the most abundant and inexpensive potential green oxidant for combustion, electrochemical reactions and total and selective catalytic oxidation processes. However, reactions with O2 are often spin‐forbidden and kinetically sluggish due to its triplet ground electronic state and strong O−O bond. To overcome this limitation, molecular O2 is activated into more reactive derivatives, such as: the dioxygen cation (O2 +), the superoxide ion (O2 −), the peroxide ion (O2 2−) and the monoatomic adoxygen species. The activation of O2 is, thus, an essential part of catalytic cycles and often the rate determining step in the oxidation processes. To effectively lower the kinetic barriers associated with O2 activation, reducible metal catalysts, rich in oxidation states and with variable coordination geometries, are required. More recent attention has focused on the mixed‐valence iron‐ and cobalt‐based oxide catalysts of the spinel and perovskite types, due to their availability, stability, and promising redox catalytic activity in many fields including the gas‐solid and liquid‐solid interfaces. Within this diverse class of compounds, cobalt ferrite, CoFe2O4 (CFO), is an attractive candidate.

CFO is an insulating material and an efficient catalyst for CO oxidation, solar water splitting, oxygen evolution (OER) and oxygen reduction (ORR) reactions. It crystallizes in a close‐packed inverse spinel cubic (Fd m) structure, where Co2+ (d7, high spin) ions occupy of the octahedral B‐sites, whereas the Fe3+(d5, high spin) ions occupy 1/8 of the tetrahedral A‐ and of the B‐sites. Depending on the sample preparation conditions, the degree of inversion x for CFO is not complete and can vary between 0.76 and 0.93, resulting in a mixed spinel, (Co Fe1‐)A[Co1‐Fe1+]BO4.

The catalytic activity of CFO and spinels in general, have been investigated using a variety of experimental techniques. Among others, these include X‐ray diffraction and transmission electron microscopy, X‐ray photoemission and absorption spectroscopy, Mössbauer spectroscopy, and X‐ray magnetic circular dichroism. Computationally, density functional theory (DFT) has become a powerful research tool for studying solid surfaces. Most commonly, the CFO crystal lattice has been theoretically modelled with supercells imposing periodic boundary conditions (PBCs) within the generalized gradient approximation with an onsite Hubbard term (GGA+U) to account for the strongly correlated 3d electrons. In this context, Hou et al. have addressed the structural, electronic and magnetic properties of bulk CFO at different degrees of inversion (x= 0.0, 0.25, 0.50, 0.75 and 1.0). Using the same approach, Hajiyani and Pentcheva, have investigated the effect of surface termination and cation substitution on the electrocatalytic performance of the CoxNi1‐xFe2O4 (001) surface (x=0.0, 0.5, 1.0) for water oxidation.

Recently, we have also employed periodic DFT+U to explore the reactivity of the pristine and reduced CFO (001)/(100) surfaces toward water.

Periodic DFT is by far the most widely used theoretical approach for treating metal oxide interfaces, however, it falls short due to high computational demands, in two scenarios. First, dealing with isolated reaction sites requires large supercells to avoid spurious interaction from the periodic pattern. Second, strongly correlated systems containing defects (e.g, oxygen vacancies) often require a fraction of exact exchange for accurate description of the electronic structure.[ , , ] To circumvent these issues, alternative methods based on embedded finite‐sized clusters are used.

The embedded cluster methodology is well‐established and several schemes have been developed and employed over the years. One such embedding scheme is the periodic electrostatic embedded cluster model (PEECM) which we employ in this work. The PEECM was developed in 2009 by Burrow et al. as a computational alternative for modelling point defects and local interactions in highly ionic crystals, such as calcium fluoride (CaF2) and ceric oxide (CeO2), that avoids problems which can arise when periodic environments are approximated with finite point charge embeddings. To date, it has been successfully applied to various metal oxides, including but not limited to actinide oxides (AnO2, An=U, Np, Pu), magnesium oxide (MgO), and ionic transition‐metal oxides with significant covalent character (e.g., ZnO and TiO2 ). PEECM divides the extended solid system into three parts: the inner part, the outer part and the isolating shell. The inner local part contains the catalytic reaction site in question as a quantum mechanically (QM) treated cluster. The outer part accounts for the lattice environment by a periodic array of point charges. At the boundary between the former two parts an isolating shell of effective core potentials (ECPs) at the cationic sites is included to prevent artificial leakage of electron density from the QM into the outer region. The PEECM is particularly efficient for the description of ionic systems, as it takes advantage of the periodic fast multipole method to include the Madelung potential in the embedded cluster scheme.

To the best of our knowledge, catalytic investigations on spinels in the limit of low coverage are still scarce. Moreover, while the reactivity of the cobalt oxide, Co3O4, surfaces towards O2 has been theoretically considered in many cases, the reactivity of the analogous Fe‐substituted spinel CFO, has not been studied in detailed yet. Motivated by this lack in research, we conduct here hybrid DFT embedded cluster calculations to elucidate the pathways for O2 activation on the (001) surfaces of spinel CFO, in contact with a gas phase. First, we model the (001) CFO surfaces with clusters of different sizes and validate them in convergence studies of adsorption potentials. We consider two distinct terminations of the (001) spinel surface, namely, the A‐ and B‐layers. The B‐layer exposes 5‐fold coordinated Co sites in a +3 valence state, whereas the A‐layer exposes 2‐fold coordinated Fe and 5‐fold coordinated Co sites in +2 valence states. To differentiate between the active Co in the A and B layer, respectively, we use the CoA and CoB naming convention. For dioxygen activation pathways, we present a comprehensive description of the interaction between O2 and the aforementioned metal active sites at the pristine surfaces, including energetic profiles, electronic structure characterization of the reactive oxygen species (ROS) and possible surface diffusion pathways. We, then, investigate the reactivity of an oxygen vacancy toward O2 activation, due to its importance in the Mars van Krevelen (MvK) mechanism.

Computational Methods

DFT calculations were carried out with the PEECM as implemented in the TURBOMOLE program package, To account for the correlated nature of the 3d electrons in CFO we used the hybrid Perdew‐Burke‐Ernzerhof exchange‐correlation functional (PBE0) in combination with the split valence basis with a single set of polarization functions (def2‐SVP), from the TURBOMOLE library. For the electronic Coulomb interaction term, we applied the resolution of identity approximation with auxiliary basis sets optimized for the corresponding orbital basis sets to improve upon the computational performance. Results obtained with the larger def2‐TZVP basis differ by only a few kJ/mol, whereas the inclusion of the dispersion correction is unreliable, with unusually large energy differences between the D3 and D4 corrected values (see Table S1 of the Supporting Information).

All single point calculations were performed with the augmented Roothaan‐Hall method for optimizing Hartree‐Fock and Kohn‐Sham density matrices. The self‐consistent field (SCF) convergence threshold was set to 1×10−6 a.u. with the m3 integration grid, whereas the cluster geometry optimization convergence criteria were set to 1×10−6 and 1×10−3 a.u. for the total energy and the maximum norm of the Cartesian energy gradient, respectively. For the O2 adsorption energies we additionally included zero‐point vibrational energy contributions calculated with the harmonic approximation, by considering the derivatives of the quadrature weights in combination with tighter convergence criteria, SCF energy: 1×10−8 a.u and DFT integration grid: m5. The Hessian matrix and the dipole gradients were computed by numerical differentiation, with a stepwidth of 0.02 a.u. for the Cartesian coordinate displacements. To obtain an initial guess for the transition states, reaction pathways were predicted using the chain‐of‐state method implemented in the woelfling module of TURBOMOLE. The transition states were further optimized using the trust‐region image minimization (TRIM) algorithm.

The 2D infinite array of point charges, aperiodic in the z‐direction, was set up by combining relaxed and bulk unit cells of 19 and 17‐layers for the (001)‐A and ‐B surfaces, respectively, into symmetric slabs. We chose to expose the two equivalent polar surfaces on both sides of the slab, while keeping the bulk‐like nature in the middle, in order to minimize the unphysical dipole moment, which will otherwise arise perpendicular to the supercell. To account for the relaxation of the (001) CFO surfaces we started from optimized geometries obtained in periodic calculations at the DFT(PBE)+U level.

Finite sized clusters, modelling the (001) CFO surfaces were treated quantum mechanically and embedded in a point charge field created by the formal charges of the respective ions, i.e., 2+ for the surface 2‐fold coordinated Fe and Co cations at the A layer, 3+ for the surface 5‐fold coordinated Co cations at the B layer, 2‐ for the oxygen anions, and 2+/3+ for bulk Co/Fe, respectively. Although CFO is strongly ionic (estimated iconicity 0.8631), previous published work has indicated that the bonding character in CFO has also a small covalent component which is essential to the cation ordering in the spinel structure. To make sure that we don't overemphasize the ionicity, we tested the clusters in embedding potentials created by fractional charges (3/4 and 1/2 of the formal charges) and found that the fastest and smoothest convergence of the results is already reached with the full formal charges.

The nearest shells of cations around the quantum clusters were augmented with 28 e− Zn/Ga ECPs. They have been chosen because the respective ions have similar radii as Co2+/3+ and Fe2+/3+. During optimization, the positions of the ECPs and the cluster atoms directly coordinated to the ECPs were kept fixed, allowing only the geometry of the most inner part of the quantum cluster to relax and adjust to the method and the introduced defects and adsorbates.

Results and Discussion

Cluster‐size validation

Three series of five hemispherical clusters, increasing in size, were built symmetrically around the Fe, CoA and CoB active sites of the A and B terminated (001) CFO surfaces, with the metal ion for the active center located as close to the cluster center as possible.

To validate the cluster models, we calculated the one‐dimensional potential energy curves for the perpendicular adsorption of carbon monoxide, CO, as a function of the metal‐C binding distance. For each cluster, we employed quadratic interpolation of the three energetically lowest points to determine the equilibrium metal‐C distance and adsorption energy. The convergence behavior of both quantities with the cluster size, for each set of clusters, is presented in Figure 1.

Figure 1

Convergence studies with the electrostatically embedded QM clusters modelling the Fe/CoA sites on the A‐layer and the CoB site on the B‐layer of the (001) CFO surface: a) convergence of the CO adsorption energy and b) convergence of the metal‐C distance with the cluster size. The size of the defected cluster with one VO1 on the (001)‐B surface has been validated in a c) convergence study of the formation energy for the vacancy. FM and AFM refer to the ferromagnetic and antiferromagnetic spin orderings of CFO.

A‐layer

The clusters with centered Fe range from 15 to 169 atoms and have the following stoichiometry: Co4Fe3O8, Co10Fe23O42, Co16Fe31O58, Co24Fe35O72, Co28Fe47O94. A separate set of clusters was built around the CoA metal site with cluster sizes between 35 and 159 atoms: Co4Fe11O20, Co12Fe15O34, Co18Fe23O52, Co22Fe35O74, Co22Fe47O90. The corresponding structures for each set of clusters are shown in Figures S1 and S2 of the Supporting Information. The aforementioned clusters are all non‐stoichiometric, as it is quite challenging to build both stoichiometric and symmetric clusters from spinel CFO.

For the A‐layer cluster models we considered two spin configurations: ferromagnetic (FM) and antiferromagnetic (AFM). The FM spin ordering of CFO with all unpaired electrons spin‐up is by 142.8 kJ/mol per formula unit less favorable than the AFM ordering with A sites spin‐up and B sites spin‐down.

Focusing on the results for the FM configuration in Figure 1a and the Fe active center, we notice that the adsorption energy changes more prominently (by 13.5 kJ/mol) between the first two clusters, indicating that the first cluster, Co4Fe3O8, is too small to correctly describe the chemical properties on this surface site. With increasing the cluster size the adsorption energy converges quickly, such that value on the largest cluster (Co28Fe47O94), −49.1 kJ/mol, is only by 0.3 kJ/mol more negative than the value on the second cluster (Co10Fe23O42), −48.8 kJ/mol. On the contrary, for the CoA active site, the adsorption energy on the smallest cluster (Co4Fe11O20) already agrees very well, within 1.5 kJ/mol, with the one on the largest cluster (Co22Fe47O90). Compared to the adsorption energy, the metal‐C binding distance in Figure 1b converges even faster and smoother for both series of clusters. The Fe−C and CoA−C bond lengths fall within the same regime, between 2.10 and 2.20 Å, and they change only little, between 0.00 and 0.07 Å, with increasing cluster size.

Adsorption of CO on the AFM clusters follows nearly identical pattern as for the FM case (see Figure 1a and Figure 1b). For the Fe active site, the adsorption energy on the second cluster, Co10Fe23O42, differs by 1.4 kJ/mol, whereas the Fe−C bond length differs by 0.04 Å, between both spin orderings. Turning next to the CoA site and the FM/AFM configurations, the differences of both quantities for the third cluster, Co18Fe23O52, are 0.69 kJ/mol and 0.01 Å. Since the FM and AFM spin orderings of CFO give similar results, we have limited the rest of our study to the lowest‐energy spin configuration (i.e., the AFM configuration: A↑/B↓).

Moreover, on the basis of this convergence behavior, we conclude that the second cluster, Co10Fe23O42, describes the chemical properties around Fe with sufficient accuracy and can be used for further investigation of local events. The surface structure of this cluster is presented in Figure 2a. For catalytic studies on the CoA site we have chosen the third cluster, Co18Fe23O52, with a sufficiently large surface area to account for the relaxation of the surrounding ions as well. The surface structure of this cluster is presented in Figure 2b and shows two CoA sites near the center. One of the CoA sites is, to some extent, ‘hindered’ by the proximity of the additional Fe cations from the A‐layer, while the other remains ‘free’ as the Fe ions relax away. To differentiate, we will use the ‘hindered’ and ‘free’ analogy and individually refer to these sites as CoAh and CoAf, respectively. It is, thus, important to mention at this stage, that the CO adsorption curves were calculated for the CoAf site.

Figure 2

Top view of the QM clusters: Co10Fe23O42, Co18Fe23O52, Co22Fe33O80, used to investigate the a) Fe2+ and b) Co2+ active sites on the (001)‐A, and the c) Co3+ active site on the (001)‐B surface terminations of CFO. For clarity, only the first and second layer atoms are represented by spheres (Fe=green, Co=nude, O=red). The O1 oxygen vacancy introduced at the (001)‐B surface is marked by a yellow dashed circle. The different types of oxygens are defined in the text.

B‐layer

The B‐layer exposing CoB sites was modeled with the following set of clusters: Co2Fe6O12, Co4Fe8O18, Co16Fe19O50, Co22Fe33O80, Co30Fe47O108 containing between 20 to 185 atoms. The cluster structures are given in Figure S3 of the Supporting Information.

In comparison to CoAf, the CO interaction with CoB is stronger, as indicated by the more negative adsorption energies in Figure 1a and the shorter Co−C bond lengths in Figure 1b. Convergence of the surface properties for the pristine cluster models of the B‐layer is achieved already with the second cluster, Co4Fe8O18.

To explore the influence of surface defects in the oxidation activity of CFO, we focused on an oxygen vacancy (VO) at the B‐layer and additionally considered a sixth cluster, Co40Fe57O136, with 233 atoms. To make sure that the aforementioned pristine quantum clusters, exposing the B‐layer, still provide a reasonable description of the surface properties when one VO is introduced, we calculated the formation energy for the defect as a function of the cluster size using the following relation:

where E(cluster+VO) and E(cluster) denote the energies of the clusters with and without a defect, respectively, and E(O2) is the energy of an isolated O2 molecule.

Of the two types of oxygens identified on the (001) surfaces of spinels, we studied only the most stable vacancy site, the O1 defect (hereafter denoted as O1B, to differentiate from the analogous O1A, in the A‐layer).[ , ] In CFO, this oxygen is bonded to an octahedral Fe in the third layer, differently from the second oxygen type, namely O2, which is bonded to a tetrahedral Fe cation in the second layer. The different types of oxygens are denoted in the CFO surface structures presented in Figure 2. Convergence of the formation energy for VO1 with the cluster size, in Figure 1c, reveals that a large surface area and cluster thickness is needed to describe this defect properly. Therefore, we chose the fourth cluster, Co22Fe33O80, to further investigate the reactivity of both, the pristine and defected (001)‐B surface, in the limit of low coverage. The surface structure of this cluster is presented in Figure 2c, in which the introduced vacancy is also marked.

The removal of an oxygen atom from the cluster leaves behind two excess electrons which are redistributed between the nearest CoB ions at the surface, reducing their valence state from +3 to +2. Finally, it is also worth to note that the density of states (DOS) for the pristine (001) – B surface as well as the VO1 formation energy calculated with PBE0/def2‐SVP and the embedded cluster approach agree rather well with those obtained recently in periodic DFT+U calculations using the PBE functional and a plane wave basis; the VO1 formation energy in the largest cluster, Co40Fe57O136, is −48.5 kJ/mol, whereas the DFT+U calculated value is 51.1 kJ/mol. The comparison of the DOS is shown in Figure S4 of the Supporting Information.

Dioxygen activation pathways

For the activation of a single O2 molecule on the (001)‐A and B CFO surfaces, we investigated two adsorption orientations: the end‐on one, with the O2 molecular axis perpendicular to the surface, and the side‐on one, with the O2 molecular axis parallel to the surface. On the pristine surfaces we considered the metal cations as activation sites, whereas on the defected B‐layer surface we focused on the activation of O2 at the anionic O1B vacancy. The adsorption energies reported in Table 1 were calculated for O2 as

Table 1

Equilibrium configuration parameters, the change in the IAO partial charges and adsorption energies for O2 adsorbed in a molecular form or as a superoxide on the pristine and defect (001) CFO surfaces.

	Pristine A				Pristine B[a]		Defect B[a]
Active site:	Monocenter Fe		Dicenter Fe−CoAh	CoAf	CoB−IS[b]	CoB−LS[c]	O1B vacancy
Conformation	End‐on	Side‐on	Side‐on	End‐on	End‐on	End‐on	Physisorption	Chemisorption
d O‐O [Å]	1.20	1.27	1.28	1.19	1.19 (1.24)	1.19 (1.25)	1.19 (1.23)	1.34 (1.38)
d metal‐O [Å]	2.24	2.09/2.04	1.97 (Fe)/2.17 (Co)	2.18	2.67 (2.21)	2.22 (1.95)	/	1.90/1.96 (Co)
Δq(O2)	−0.142	−0.512	−0.739	−0.067	−0.002	0.020	−0.011	−1.021
Δq (active site)	0.096	0.288	0.375 (Fe)/0.020 (Co)	0.017	0.008	−0.072	/	/
Δq (A site)	<10−3	<10−3	<10−3	0.002	<10−3	0.002	0.006	0.014
Δq (B site)	0.003	0.005	0.001	0.004	<10−3	<10−3	<10−3	0.148/0.238(Co) [d]
Δq O1	0.010	0.052	0.131	0.002	0.003	<10−3	<10−3	0.026
Δq O2	0.006	0.021	0.068	0.018	<10−3	0.023	<10−3	−0.001
ΔEads(O2) [kJ/mol]	−23.8	−49.1	−15.6	16.9	−16.4 (−21.0)	6.2 (−2.4)	−21.4 (−2.9)	26.7 (−41.7)

[a] values in parentheses correspond to periodic DFT+U results. [b] CoB in intermediate spin configuration. [c] CoB in low‐spin configuration. [d] values refer to the CoB ions directly interacting with O2.

/

/

/

where E(cluster+O2) and E(cluster) denote the energies of the QM cluster with and without O2 adsorbate, respectively, E(O2) is the energy of an isolated O2 molecule in the gas phase, and ΔZPVE is the zero‐point vibrational energy difference contribution. The energies for the individual monoatomic adoxygen adsorptions reported in Table 2 were calculated per oxygen atom with respect to half of the energy of molecular O2.

Table 2

Equilibrium configuration parameters, the change in the IAO partial charges and adsorption energies per oxygen atom for the monoatomic metal−O ROS on the A and B layers of the (001) CFO surface.

	Pristine A				Pristine B
Active site:	Monocenter Fe		Dicenter Fe‐CoAh		CoB
	Oadatom1	Oadatom2	Fe−O	CoAh−O	CoB−O
d metal‐O [Å]	1.59	2.00	1.86	1.77	1.62
Δq(O)	−0.646	−0.561	−0.634	−0.572	−0.263
Δq (active site)	0.279		0.361	0.280	−0.044
Δq (A site)	0.042[a]		0.001		−0.002
Δq (B site)	0.002		0.004		0.006
Δq O1	0.086		0.088		0.061
Δq O2	0.273		0.091		0.064
ΔEads(O) [kJ/mol]	83.4[b]		56.4	115.1	146.9

[a] value reported for the third layer Fe atom in tetrahedral site, interacting with one of the oxygen adatoms. [b] average adsorption energy per oxygen atom.

End‐on adsorption of O2 on the Fe site results in an energy gain of 23.8 kJ/mol and a monodentate Fe‐(η 1)O2 binding. The optimized geometry, depicted as a in Figure 3, is characterized with a tilted O2 molecular axis at approximately 50° from its initial perpendicular orientation. A slight elongation of the O−O bond from 1.19 Å (isolated O2 molecule) to 1.20 Å indicates a weak chemical interaction between iron and oxygen, with a bonding distance d Fe‐O=2.24 Å.

Figure 3

The two‐step activation pathway for O2 on the (001)‐A surface of CFO. The monocenter activation sequence on the Fe site is presented in bold and with small letters, whereas the dicenter activation on the Fe‐CoAh cation pair is denoted with capital letters. The transition states are marked with #. Adsorbed O2 in yellow is distinguished from lattice oxygen in red. For clarity, only the atoms from the topmost layers and those involved in direct interaction with the adsorbed O2 are represented by spheres.

From the initial end‐on configuration a, O2 can rearrange to structures where it is bonded in a side‐on mode and activated either on the Fe site (monocenter activation) or in tandem on the surface Fe and CoA cations (dicenter activation). Detailed analysis at the molecular level, presented in Figure 3, reveals that both activation pathways may be divided into two individual steps: elongation of the O−O bond and its subsequent dissociation. To differentiate the two pathways for rearrangement, we have presented the monocenter activation on the Fe metal center in bold and labeled the corresponding structural conformations with small letters. By contrast, the dicenter activation sequence on the Fe−CoA cation pair is labeled with capital letters. In both cases, the transition states are marked with #. Adsorbed oxygen – in yellow color is distinguished from the lattice oxygen ‐ red.

Monocenter activation on Fe: Starting from the molecularly adsorbed O2 in an end‐on mode, the suprafacial activation of O2 on Fe has to initially surpass an energy barrier of 11.7 kJ/mol before reaching a local minimum of −49.1 kJ/mol (sequence a→a→b). In this scenario, the O−O bond length (d O‐O) is stretched to 1.23 Å upon reaching the transition state geometry (a) and 1.27 Å for the resulting minimum (b), indicating the formation of a superoxide, O2 −. The superoxide anion in side‐on geometry is stabilized almost symmetrically on top of Fe featuring a bidentate chelating Fe‐(η 2)O2. Upon further stretching, the O−O bond is weakened and finally broken, producing two monoatomic oxygen species (sequence b→b→c). This step, however, requires a high activation energy of 250.6 kJ/mol as the transition state configuration (b) shows interaction between the third layer Fe cation in tetrahedral geometry and one of the oxygen adatoms, d Fe‐O=2.52 Å. The final product (c) is characterized by two clearly different Fe−O bond lengths and is energetically less stable than its parent superoxo adduct (b) by 224.2 kJ/mol. One of the O adatoms is pulled strongly by the hosting Fe site forming a metal‐oxo species with a bonding distance of 1.59 Å, whereas the other O atom is additionally stabilized by the interaction with the surrounding ions. The oxygen adatoms are separated by 2.74 Å. The reverse step of oxygen recombination features an energy barrier of 34.8 kJ/mol.

Dicenter activation on Fe−Co: The two‐step activation of O2 on the dual Fe−CoA center, in Figure 3, follows a similar pathway as for the monocenter activation on Fe, however, the corresponding adsorption geometries lie higher in energy. In the first step (sequence a→A#→B), the bond in molecular oxygen is stretched to 1.24 Å and O2 is shifted towards the closest CoAh in the transition state, before relaxing to a minimum with a superoxo nature of the adsorbed O2 (d O‐O=1.28 Å). In the minimum (B), the bidentate bridging Fe−(μ)OO−CoAh interaction is characterized by d Co‐O=1.97 Å. This structure is by 8.2 and 33.5 kJ/mol less stable than the starting molecular conformation (a) and the superoxide formation on Fe (b), respectively. The calculated barrier height is 12.9 kJ/mol. Further stretching of the O−O bond and subsequent dissociation into monoatomic species (sequence B→B#→C) requires a high activation energy of 243.9 kJ/mol. In the transition state (B#), stronger Fe−O (d Fe‐O=1.77 Å) and CoAh−O (d Co‐O=1.80 Å) interactions take place while the O adatoms are located 2.21 Å apart. This intermediate product relaxes to the final geometry (C) in which two metal‐oxyl species, [FeIII−O]2+ (d Fe‐O=1.86 Å) and [CoIII−O]2+ (d Co‐O=1.77 Å), are formed. In comparison, stabilization of an oxygen adatom on the Fe ion is an endothermic process with an energy cost of 56.4 kJ/mol, more favorable than the formation of a [CoIII−O]2+ species, which requires 115.1 kJ/mol. Recombination of these monoatomic oxygen species into dioxygen requires an activation energy of 40.9 kJ/mol.

From these findings, it is clear that the adsorption of dioxygen, from a gas phase, has a much lower activation energy than the subsequent dissociation into reactive monoatomic species. For the latter to occur, the temperature of the system should be high enough to overcome the large energy barriers presented in Figure 3. As a result, the activation of O2 on the (001)‐A surface of CFO may be the thermodynamic rate‐controlling step in the oxidation processes occurring via suprafacial, Langmuir‐Hinshelwood and Eley‐Rideal, mechanisms.

Activation on Co: As dioxygen approaches the surface in a perpendicular end‐on fashion on top of the CoAf metal center, some interaction takes place, which, in its final state, is characterized as monodentate CoAf‐(η 1)O2 binding with d Co‐O=2.18 Å (Figure S5 of the Supporting Information). Such association, however, is energetically not favorable as indicated by the corresponding positive adsorption energy of 16.9 kJ/mol.

Starting from side‐on geometries with elongated O−O bond length close to that of superoxide or peroxide species, on top of a single CoAf center or bridging a CoA−CoA cation pair, geometry optimizations always converged to the previously described CoAf−(η 1)O2 attachment. Tandem activation on the Fe−CoAf dicenter (i.e., the Fe−(μ)OO −CoAf adsorption conformation) could not be found, most probably due to the large distance between the two cations, d Fe‐Co=3.86 Å, as a consequence of the strong relaxation pattern in the A layer.

Electronic structure of surface oxygen species: For a thorough description of the nature of the identified reactive oxygen species formed on the A layer, we computed intrinsic bond orbitals (IBOs). We employed the approach in Ref. [39] in combination with spin density analysis to determine the oxidation states. For completeness, we also calculated the intrinsic atomic orbital (IAO) partial charges at the structures for the intermediates reported above. We note, however, that partial charges are poor predictors of the oxidation states and cannot be solely taken as indicative of the binding motive especially when some covalent character is present, as reported in many previous cases. Table 1 lists the equilibrium configuration parameters, the change in the IAO partial charges and adsorption energies for O2 adsorbed in a molecular form or as a superoxide on the pristine and defect (001) CFO surfaces. For the monoatomic metal−O interactions, these parameters are summarized in Table 2. The change of the IAO partial charges for the O2/O adsorbates, adsorption site, the surface Fe and Co sites in tetrahedral and octahedral geometries, are denoted by Δq(O2)/Δq(O), Δq (active site), Δq (A site) and Δq (B site), respectively, whereas Δq O1, Δq O2 denote the respective changes for the different surface oxygen types, namely, O1A/B and O2 (see Figure 2).

The bonding region around the O−O superoxo moiety on top of the Fe monocenter (structure b) is presented in Figure 4a. It shows six and seven singly occupied valence molecular orbitals, in α and β shells, respectively, located at the O−O moiety. Although d O‐O=1.27 Å of adsorbed O2 − is by 0.05 Å smaller than d O‐O=1.32 Å of the isolated O2 − in the gas phase, the overall configuration reveals a two‐center, two electron σ ‐ bond and a two‐center, one electron π‐bond, clearly confirming the presence of a one‐electron reduced superoxide anion. Bonding of this O2 − to the surface is coordinative in nature, such that electron density from the initially non‐bonding and orthogonal pα orbitals of O2 is partially donated to the Fe adsorption site. On the other hand, only five spin β valence IBOs fully localize at the Fe center suggesting that this cation acts as an actual electron donor and is oxidized to FeIII with a high‐spin configuration (see Figure S6 of the Supporting Information). The oxidation process is associated with an increase of the IAO charge of Fe by Δq Fe=0.288, whereas the superoxide develops a negative charge of Δq(O2)=−0.512.

Figure 4

Singly occupied IBOs localized around the molecularly activated dioxygen at the a) Fe monocenter and b) Fe−CoAh cation pair on the A‐layer. Isosurface value is ±0.05.

For the O2 − moiety stabilized on the Fe−CoA pair (structure B in Figure 3) the IBO analysis is very similar to the previous superoxo case described above. The pictorial representation of the localized IBOs, in Figure 4b, features a two electron σ‐bond and one electron π‐bond, which is reflected in an elongation of the O−O bonding distance to 1.28 Å. Once again, the coordinative association is attributed to the donation of the O2 electron density from the pα orbital to the Fe center. Upon bonding, one electron from Fe is transferred to the O2 molecule and Fe is oxidized to FeIII, as revealed by five IBOs localized at this Fe atom and consisting mainly of its 3d orbitals (Figure S7a of the Supporting Information). At the same time, the occupation numbers of seven 3d spin‐orbitals belonging to CoAh in the coordinated environment remain at n≈1, indicating that this site holds on to its original CoII (3d 7) state (Figure S7b of the Supporting Information). The superoxo moiety is negatively charged (Δq(O2)=−0.739), and the majority of electrons is donated from the hosting Fe and the O1A surface sites (Δq Fe=0.375 and Δq O1=0.131).

A detailed representation of the electron spin density distribution for the metal‐O monoatomic species formed upon O2 − dissociation is shown in Figure 5, whereas the bonding nature of these species in terms of localized IBOs is collated in Figure 6. Focusing on the Fe monocenter, in Figure 5a, the spin repartition on oxygen in the Fe−O interaction with shorter bond length, d Fe‐O=1.59 Å, is almost fully quenched, indicating an O2− species with all electrons paired. Correspondingly, eight valence IBOs (four with α and four with β spin) localize at the adsorbed O (see Figure 6a). The overall coordination features both, a σ covalent character and π‐bonding. The major contribution (with over 70 %) in the orbital overlap of σα character comes from the O adatom, different from the σβ counterpart with the major contribution from the Fe site (over 60 %). Additionally, the O adatom acts as a π‐donor, as indicated by the electron density flow from two orthogonal pα orbitals to two Fe d orbitals. By contrast, in the Fe−O attachment with d Fe‐O=2.00 Å, the suprafacial O carries some spin density, corresponding to its reduced O− state. The IBO analysis, in Figure 6b, reveals an orbital overlap of σ symmetry between O and the hosting Fe site. The nearest oxygens from the lattice and the Fe(Td) cation form the third layer also engage in stabilizing this adoxygen. The O adatom interacts with O2 from the lattice covalently and the Fe(T ) orbital contribution in this mixing (characterized with 12 %) results from the ligand π ‐ donation. This bonding motive is reflected in the IAO partial charges as well, as inferred from the charge loss on the O2 lattice oxygens (Δq O2=0.273) in Table 2. The O2− moiety (Δq O=−0.646) is more negatively charged than O− (Δq O=−0.561), and most of this charge is transferred from the Fe adsorption site (Δq Fe=0.279).

Figure 5

Spin density contours for the antiferromagnetic cluster models of CFO with monoatomic adoxygen species formed upon reductive O2 dissociation on the: a) Fe monocenter and b) Fe−CoAh cation pair at the A‐layer. c) Spin density contour for the antiferromagnetic cluster of CFO, exposing the B‐layer with associatively adsorbed O adatom on the CoB site.

Figure 6

IBOs localized around the O adatom in the metal‐O interactions on the Fe monocenter: a) Fe−O (d Fe‐O=1.59 Å), b) Fe−O (d Fe‐O=2.00 Å) and on the Fe−CoAh dicenter: c) Fe−O (d Fe‐O=1.86 Å), d) Co−O (d Co‐O=1.77 Å) at the A‐layer, and on the CoB center e) CoB−O (d Co‐O=1.62 Å) at the B‐layer. For each molecular orbital the major contribution is presented by the atomic symbol and the corresponding value in parentheses. Isosurface value is ±0.05.

Turning next to the O2 dissociation on the Fe−CoAh cation pair, we note that in both interactions, Fe−O and CoAh−O, the suprafacial oxygens carry spin density (Figure 5b). In the Fe−O coordination, the oxygen atom does not satisfy the octet rule featuring seven IBOs, as shown in Figure 6c. The bonding is primarily of σ character, which originates from mixing of the oxygen 2pz and iron dz 2 orbitals. Additionally, a side π‐donation from oxygen to the hosting Fe is also present. In comparison to the previously described iron‐oxo species with d Fe‐O=1.59 Å, the π‐bonding interactions between Fe and oxygen are reduced, which is reflected in an elongation of the Fe−O bond distance by 0.27 Å, as expected for the induction of metal‐oxyl character, [FeIII−O]2+.

The CoAh−O interaction is characterized with eight IBOs presented in Figure 6d. Again, the coordination has both, σ and π character. The σ character is attributed to the overlap between the oxygen 2pz and cobalt dz 2 orbitals. In one of the π‐donations the electron density flows from the orthogonal pβ orbital corresponding to oxygen (total contribution around 60 %) to the CoAh center. By contrast, in the second case we observe a π‐back donation from the CoAh metal center to the O adatom. The overall CoAh−O interaction can be better described as a resonance hybrid of the Co‐oxo and Co‐oxyl species. Because cobalt is located behind the ‘oxo‐wall’ (group 9) its 3d electron configuration in tetragonal geometry does not allow the formation of a terminal Co=O double bond.[ , ] As a consequence, spin distribution is localized on the oxygen ligand, suggesting a higher degree of oxyl character, i.e., [CoIII−O]2+ species.

In both metal‐oxyl species, [FeIII−O]2+ and [CoIII−O]2+, the majority of electrons in the negatively charged O adatom is donated by the hosting metal center (Δq Fe=0.361 and Δq Co=0.280).

Surface diffusion of monoatomic reactive oxygen species (ROS): Once produced, the monoatomic adoxygen species can migrate and reside on different adsorption sites on the surface lattice. The number of successive jumps depends on the nature of the surface. For the A layer, we investigated three possible routes which are illustrated in Figure 7a. In route 1, we focused on the successive hopping of the O adatom to adjacent metal sites, whereas in routes 2 and 3, we studied the diffusion of adoxygen from the CoA centers to the neighboring O2 and O1A lattice oxygens, respectively.

Figure 7

a) Possible diffusion routes of the monoatomic reactive oxygen species across the A‐layer of the (001) CFO surface. b) Detailed diffusion pathway for the O adatom migration at adjacent metal cationic sites. c) Detailed pathway for the cationic – anionic surface diffusion of the O adatom. The transition states are marked with #.

The stepwise diffusion pathway for route 1 is displayed in Figure 7b. The hopping process starts with adoxygen on the Fe center, which, as we previously noted, is also the most stable metal‐oxyl species on the A layer. Hopping of the O adatom onto the closest CoA (sequence a→a→b) is hindered by a large activation barrier (ΔE=100.1 kJ/mol), which can be attributed to both, the pronounced stability difference between the initial and final states and the distance between the metal sites (d Fe‐Co=3.12 Å). In the transition state (a) the Fe−O bond is stretched to 2.62 Å, followed by strong CoAh−O interaction (d Co‐O=1.84 Å) until the O adatom is completely shifted to CoAh (b). Subsequent diffusion of the adoxygen onto the CoAf site in the b→b→c→c→d scenario, is not direct but via an additional local minimum (c), in which the O adatom bridges the two CoA sites, CoAh−(μ)O−CoAf. This structure (c), is by 6.1 kJ/mol less favorable than the parent state (b). From the thermodynamic point of view, the CoAh−(μ)O−CoAf conformation and the stabilization of the O adatom on the CoAf site are nearly identical with similar adsorption energies of 121.2 and 120.8 kJ/mol, respectively. Consequently, the CoAf−O specie (d) is less stable than CoAh−O in structure (b) by 5.7 kJ/mol. The energy barriers associated with this two‐step process are significantly smaller than for the Fe→CoAf jump, indicating that mobility of the adoxygen from one CoA site to the other is relatively facile. In the final step (sequence d→d→e), the O adatom is shifted from CoAf to the Fe site. The energy barrier (ΔE=74.0 kJ/mol) is smaller than for the Fe→CoAf transition, as the structural tension, to some extent, is compensated by the stability of the Fe−O species. The corresponding transition state (d), however, lies higher in the potential energy profile due to the larger d Fe‐Co=3.86 Å distance that the O adatom needs to overcome in order to diffuse.

The detailed order of mixed cationic‐anionic adoxygen diffusion in route 2, is presented in Figure 7c. Initially residing on the CoAh site, the adoxygen is shifted toward O2 and captured in a local minimum by the proximity of the Fe cations (sequence A→A#→B). The resulting structure (B) is by 14 kJ/mol more favorable than O residing on CoAh solely, which can once again be attributed to the electrostatic interaction with the Fe site (d O‐Fe=2.07 Å). The activation energy needed for this diffusion step is only 2.8 kJ/mol, additionally suggesting that the O adatom prefers bridging the Fe and CoAh metal centers, rather than shifting towards CoAf. Further relocation of the O adatom to the O2 anionic site is associated with a large energy barrier of 96.6 kJ/mol and the formation of an Oadatom−O2lattice species with a bond length of 1.39 Å, as presented in sequence B→B#→C. Among all the possibilities investigated so far, stabilization of adoxygen on the anionic O2 site is least favorable as it requires an energy cost of 154.1 kJ/mol. Furthermore, the O adatom does not reside on the O1A anionic site (route 3), instead it readily migrates to the closest Fe from the A layer.

As we can see, the barriers associated with the diffusion across the (001)‐A surface of CFO of the monoatomic ROS produced upon reductive O2 dissociation are smaller than those corresponding to the dissociative adsorption and the associative desorption. The monoatomic ROS may, therefore, migrate and reside at different sites on the A layer, of which the Fe metal center remains more favorable.

Activation on the pristine surface: On the (001)‐B surface of CFO, dioxygen adsorbs only in an end‐on fashion on top of the CoB site. The corresponding optimized structure, presented in Figure 8a, shows a CoB‐(η 1)O2 association with a binding energy of −16.4 kJ/mol and d Co‐O=2.67 Å. The CoB underneath O2 remains in an intermediate spin configuration, whereas adsorption on a CoB in a low spin state is less favorable by 22.6 kJ/mol (see Figure S8 of the Supporting Information).

Figure 8

Geometric structures of: a) dioxygen adsorbed in a molecular form on top of CoB with intermediate‐spin configuration, b) O adatom adsorbed on the CoB site at the pristine (001)‐B surface. c) Physisorbed molecular O2 and d) chemisorbed dioxygen as superoxide in the O1 vacancy at the defected B‐layer. The O1 vacancy is marked by a yellow dashed circle.

Similar as for the activation on CoAf at the A‐layer, we do not observe any molecularly activated O2 − or O2 2− species on top of the CoB monocenter or bridging a CoB−CoB cation pair. Binding of an O adatom to CoB at the bare (001)‐B surface is only favorable at elevated temperatures as indicated by the positive adsorption energy of 146.9 kJ/mol. The optimized structure is presented in Figure 8b. In comparison to the metal‐O association at the A layer, CoB−O is less stable than the reactive CoAh−O and Fe−O species, by 31.8 and 89.8 kJ/mol, respectively.

The IBO analysis around the O adatom at the B‐layer, is presented in Figure 6e and features eight molecular orbitals. The σ character in both, α and β counterparts, is dominantly oxygen‐based, whereas the bonding π(CoB−O) orbitals are dominantly CoB‐based d orbitals. The spin density (in Figure 5c) and charge analysis (Table 2), reveal that the O adatom is polarized and negatively charged (Δq O=−0.263), whereas most of the spin repartition on the hosting CoB site is quenched and this site obtains electrons (Δq Co=−0.044). The majority of charge is transferred from the surface oxygens (Δq O1=0.061 and Δq O2=0.064). Differently from the A‐layer, the stabilization of adoxygen in this case, can be attributed to the electrons flowing within the surface ions, rather than the reducing ability of the hosting metal site (i.e., Co3+). The overall CoB−O coordination is better described as CoIII‐oxyl with d Co‐O=1.62 Å and cobalt in a low‐spin configuration. The excess unpaired electron responsible for the spin density distribution around the O adatom, is located in an out‐of‐plane p orbital with respect to the Co−O axis obeying the ‘oxo‐wall’ paradigm.

Activation at a vacancy: Physisorption of dioxygen on top of a VO1 vacancy at the B layer results in ΔEads=−21.4 kJ/mol for the end‐on conformation (Figure 8c) and ΔEads=−18.3 kJ/mol for the side‐on conformation. Adsorption in the vacancy is not favorable (ΔEads=26.7 kJ/mol) and it is accompanied by an elongation of the oxygen bond to d Co‐O=1.32 Å. During this process, one of the oxygen atoms binds to a single surface CoB site, while the other binds to CoB from the top layer and Fe(O ) from the third layer (Figure 8d). The IBO analysis around the O2 moiety is presented in Figure S9 of the Supporting Information and corresponds to a superoxo species which is partially reduced to a peroxide. Subsequent dissociation of dioxygen into monoatomic species, such that the oxygen pointing outwards is transferred on top of the surface CoB, while the other refills the vacancy, is energetically expensive; it requires an energy cost of ΔEads=106.5 kJ/mol (Figure S10 of the Supporting Information).

As a consequence, competing molecules in the gas phase may readily adsorb and dissociate at vacancies, preventing the interaction of O2 with the anionic defects. These results indicate that surface re‐oxidation at low temperatures in the gas phase is difficult and may be the thermodynamic limiting step for oxidation reactions based on the Mars van Krevelen (MvK) mechanism.

For catalytic applications, it is always interesting to compare the reactivity of different analogous surfaces. Therefore, in the following we will discuss how the reactivity of the (001) surface toward dioxygen adsorption in the gas phase differs among the spinel family of compounds. Focusing on the pristine B layer, adsorption of O2 on top of the 5‐fold coordinated Co of normal spinel Co3O4 is characterized with ΔEads=−10.6 kJ/mol, which is comparable to the adsorption on CoB of CFO (−16.4 kJ/mol). On the Ni and Co metal sites of inverse spinel nickel cobaltite (NiCo2O4), the adsorption of dioxygen is not favorable with ΔEads of 16.4 and −1.0 kJ/mol, respectively.[ , ] Turning next to the defected (001) – B layer with one VO1, adsorption of O2 at the vacancy is an exothermic process characterized with the formation of a superoxide and an energy gain of 40.5, 37.6 and 21.2 kJ/mol for the reduced Co3O4, Ni terminated NiCo2O4,[ , ] and NiFe2O4 surfaces, respectively. By contrast, our calculated adsorption energy for dioxygen at the O1 vacancy in CFO is positive, ΔEads=26.7 kJ/mol.

Such striking difference in the reactivity of the reduced CFO in comparison to other spinels, is most probably related to the computational approaches used. The cited data on the abovementioned spinels is derived from periodic DFT calculations with on‐site Hubbard U correction, whereas we employ hybrid DFT with 25 % Hartree‐Fock exact exchange. To verify that this discrepancy originates from the two computational methods, we used our previous setup and calculated the adsorption energies at the DFT(PBE)+U level for dioxygen on the pristine and reduced (001) ‐ B surfaces of CFO. These results are incorporated in Table 1 as well, whereas the relaxed geometries are given in Figure S11 of the Supporting Information.

We observe similar energetic trends for the end‐on CoB−η1−O2 adsorption on the pristine B layer. The periodic DFT+U predicts the adsorption of O2 on a CoB with an intermediate spin configuration to be favorable and by 18.6 kJ/mol more stable than on a CoB with low‐spin, completely in line with the hybrid DFT findings. However, the situation is different for the reduced B surface. While the optimized structures agree well between both approaches, the DFT+U setup gives a negative adsorption energy of −41.7 kJ/mol for the formation of a superoxide at the vacancy, while the O2 physisorption on top of the reduced surface is characterized with −2.9 kJ/mol. It is clear that reduced CFO toward O2 behaves in a similar manner as the other spinel‐related systems at the DFT+U level and do not reflect the hybrid DFT results. Therefore, the controversy in the reactivity can be attributed to the electronic structure method rather than the cationic substitution in the spinel crystal lattice. It is worth noting that various studies[ , , , , ] so far have shown that DFT+U does not always mimic the hybrid density functionals. It seems that this is mostly related to the choice of the optimal effective U parameter, which is very dependent on the particular system and/or the property of interest and there is only a little experience on the transferability of the U parameter from a nondefective to a defective system.[ , ]

Conclusion

The direct activation of molecular oxygen, as the most abundant and green oxidant, is often an elementary step of oxidation processes in heterogeneous catalysis. To effectively lower the kinetic barriers associated with O2 activation, open‐shell metal oxides with easily accessible electrons are required. In this work, we investigated the catalytic activity of the A and B terminations of the (001) surface of spinel CFO using DFT(PBE0) within the PEECM, in contact with a gas phase.

On the A‐layer, O2 is activated either on the Fe monocenter or bridging the Fe−CoAh cation pair. A molecular level description reveals a two‐step pathway: elongation of the O−O bond to form a superoxide and subsequent dissociation into monoatomic ROS. The latter is always associated with high activation energies and may be the rate‐determining step in the oxidation processes. By contrast, our computational findings indicate that the activation of O2 on the B‐layer occurs at the most stable anionic vacancy ‐ the O1B defect. Dissociation and refilling of the defect is quite expensive energetically, thus, restoring the original stoichiometry and composition of the CFO catalyst may be the rate‐controlling step for oxidation processes based on the intrafacial MvK mechanism.

By analyzing the detailed electronic structure of the reactive adoxygen species we identify the activation as one electron processes, in which the electrons are provided by and transferred from the hosting d metal centers. The interaction of the O adatoms with the Fe, CoA and CoB sites is coordinative in nature based on electron density donation from O to the metal and vice versa. Of the monoatomic ROS produced upon reductive O2 dissociation at the A‐layer, the [FeIII−O]2+ species is thermodynamically more stable (56.4 kJ/mol), whereas the [CoIII−O]2+ species are relatively mobile with barriers smaller than those of the associative recombination.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.