Activating the FeS (001) Surface for CO₂ Adsorption and Reduction through the Formation of Sulfur Vacancies: A DFT-D3 Study

Abstract

As a promising material for heterogeneous catalytic applications, layered iron (II) monosulfide (FeS) contains active edges and an inert basal (001) plane. Activating the basal (001) plane could improve the catalytic performance of the FeS material towards CO₂ activation and reduction reactions. Herein, we report dispersion-corrected density functional theory (DFT-D3) calculations of the adsorption of CO₂ and the elementary steps involved in its reduction through the reverse water-gas shift reaction on a defective FeS (001) surface containing sulfur vacancies. The exposed Fe sites resulting from the creation of sulfur vacancies are shown to act as highly active sites for CO₂ activation and reduction. Based on the calculated adsorption energies, we show that the CO₂ molecules will outcompete H₂O and H₂ molecules for the exposed active Fe sites if all three molecules are present on or near the surface. The CO₂ molecule is found to weakly physisorb (−0.20 eV) compared to the sulfur-deficient (001) surface where it adsorbs much strongly, releasing adsorption energy of −1.78 and −1.83 eV at the defective FeS (001) surface containing a single and double sulfur vacancy, respectively. The CO₂ molecule gained significant charge from the interacting surface Fe ions at the defective surface upon adsorption, which resulted in activation of the C–O bonds confirmed via vibrational frequency analyses. The reaction and activation energy barriers of the elementary steps involved in the CO₂ hydrogenation reactions to form CO and H₂O species are also unraveled.

1. Introduction

The conversion of CO₂ into fuels and chemicals is a promising process towards achieving green energy [1,2,3,4]. However, efficient, and inexpensive catalysts are needed to facilitate the conversion of CO₂, which will also reduce the environmental impacts of carbon emissions [5,6]. In submarine hydrothermal environments, the surfaces of ferrous sulfide minerals are posited to have acted as heterogeneous catalysts for CO₂ reduction, catalyzing the reaction of CO₂ and H₂ to form small organic molecules in the primordial ocean [7,8,9,10,11]. It has been suggested by Russell et al. that iron sulfide phases such as mackinawite (FeS), greigite (Fe₃S₄), and violarite (FeNi₂S₄) may have played important catalytic roles in prebiotic chemistry [9]. The surfaces of these minerals are structurally similar to the primary active sites, (Fe,Ni)S clusters, of natural enzymes such as carbon monoxide dehydrogenase (CODH), which efficiently and reversibly catalyze the reduction of CO₂ to CO [12,13]. Huber and Wächtershäuser have shown that it was feasible to synthesize acetic acid on sulfide surfaces under the conditions of early Earth [14], whereas recent studies have shown that greigite (Fe₃S₄) and pyrrhotite (Fe_1−xS) can electrocatalytically reduce CO₂ to methanol, formic, and acetic acid under room temperature and pressure [15,16].

The activation of the CO₂ molecule is the key step for its catalyzed reduction reactions [17] and detailed understanding of the mechanisms involved in CO₂ conversion is indispensable in the design of efficient catalysts. Depending on the surface structure and electronic properties, the activity of catalyst surfaces may differ significantly with surfaces exhibiting low work functions showing more promise for CO₂ reduction as they favour electron transfer [17,18]. The nature of the surface structures, compositions, and electronic properties of iron sulfides, that facilitate CO₂ activation and reduction under hydrothermal conditions, remain poorly understood. Whereas synthetic mackinawite (FeS) is often reported to be stoichiometric, naturally occurring FeS is known to be non-stoichiometric, either due to the presence of S vacancies or interstitial metal atoms intercalated within its layers [19,20].

Compositional variations of mackinawite have been determined by several authors. Clarke et al. summarized some of the earlier results and suggested a composition in the range of MS_0.93–MS_0.96, where M = Fe²⁺ + Ni²⁺ + Co²⁺, with nickel making up to 10% by weight [21]. Results based on electron microprobe analysis of Ni and Co-free mackinawite, gave ratios of FeS_0.9–FeS_0.946 [22]. Ward et al. concluded that phase-pure mackinawite can be represented by the formula FeS_0.94 [23]. Taylor and Finger et. al. found crystallographic evidence confirming analytical suggestions that mackinawite is sulfur-deficient (FeS_1−x, typically 0 ≤ x ≤ 0.07) [24].

The (001) basal plane is by far the thermodynamically most stable surface of FeS and it encloses the largest area of its crystal morphology, which has a thin and tabular form [25,26]. Despite its large surface area, the stoichiometric FeS (001) surface is unreactive towards adsorbing species, such as H₂O, NO_x and CO₂ [26,27,28,29]. However, as the experimental material is sulfur-deficient, here we investigate the effects of sulfur vacancies in the FeS (001) surface on the activation of CO₂ and its hydrogenation reduction reactions. The existence of oxygen vacancies at metal oxide surfaces has been demonstrated to promote CO₂ adsorption and its activation [30,31,32]. Recent density functional theory (DFT) investigations, for example, have shown that bent CO₂⁻ species are formed at the Zn-terminated ZnO (0001) surface via inserting one of the O atoms of the CO₂ molecule into the vacancy [31,32,33]. The presence of surface oxygen vacancy sites on In₂O₃ (110) promote the formation of methanol via CO₂ hydrogenation reactions [34]. Also, the adsorption of CO₂ has been demonstrated to be thermodynamically more favourable at reduced CeO₂ (110) surface containing oxygen vacancies than at the stoichiometric CeO₂ (110) [35]. Compared with the numerous studies on the reduction of CO₂ on metal oxide catalysts [36,37,38,39,40,41,42], only a few studies have reported on metal sulfides [15,16,29].

The present study presents a comprehensive DFT mechanistic investigation of the activation and reduction of CO₂ on defective FeS (001) surface. We first examined the energetics of sulfur vacancy formation on the FeS (001) surface, subsequently investigating the adsorption and electronic properties, including charge transfer and vibrational frequency assignments for the CO₂ activation. Finally, we discuss the thermodynamics and calculated activation energy barriers of the elementary steps involved in CO₂ reduction through the reverse water–gas shift reaction (CO₂ + H₂ → CO + H₂O) on the defective FeS (001) surface.

2. Results and Discussion

2.1. Bulk and Surface Characterisation

FeS adopts the tetragonal structure (space group P4/nmm) with lattice parameters a = b = 3.674 Å, c = 5.033 Å, and c/a ratio = 1.370 Å [43,44,45] as shown in Figure 1. The structure is composed of vertically stacked two-dimensional FeS layers in the c-direction that are held together by weak Van der Waals forces [43,44,45]. A full unit relaxation yields the unit cell parameters a = b = 3.615 Å, c = 5.001 Å, and c/a ratio = 1.383 Å, all of which are in good agreement with commonly reported experiment values [43,44,45].

Displayed in Figure 2a is a (3 × 3) supercell of the perfect FeS (001) surface, showing that each iron atom is arranged in square-planar coordination with neighbouring irons, and the sulfur atoms are in an asymmetric one-sided four-fold coordination with iron. A single (V_S1) or double (V_S2) sulfur vacancy on the FeS (001) surface was created by removing one or two surface sulfur atoms from the perfect FeS (001) surface, as shown in Figure 2b,c. We observe small contractions in the exposed square-planar Fe–Fe distances in the defective (V_S1) surface site, calculated at an average value of 2.466 Å compared to the 2.554 Å in the prefect stoichiometric (001) surface. For the V_S2 surface, we observed both contraction and expansion in the exposed square-planar Fe–Fe distances, with the shortest and longest Fe–Fe distances calculated at 2.379 and 2.731 Å, respectively.

The formation energy per vacancy (E_form) of the surface containing n sulfur vacancies is defined as:

$$E_{form}=E_{F{e}_{36}{S}_{36-n}}+n{\mu }_S-E_{F{e}_{36}{S}_{36}}$$

(1)

where $E_{F{e}_{36}{S}_{36}} \mathrm{and} E_{F{e}_{36}{S}_{36-n}}$ are the total energies of a (3 × 3) FeS supercell without and with n sulfur vacancies. The chemical potential of S atom (${\mu }_S$), is taken under S-rich condition to be equal to the ${\mu }_S^0$ calculated from orthorhombic S₈. The formation energy of a single sulfur vacancy (V_S1) in the (3 × 3) FeS supercell is found to be 1.53 eV, while that for two sulfur vacancies (V_S2) adjacent and far from each other is 1.71 eV and 1.85 eV, respectively. Although the sulfur vacancy formation energies are endothermic, the energy required is not excessively high, suggesting that sulfur vacancies can form in the FeS basal plane. Moreover, we found that the creation of S vacancies lowers the work function (Φ) of the perfect (001) surface from 4.72 eV to 4.33 eV for the single S vacancy surface and down to 4.06 eV for the double S vacancy surface. The reduced work function has implications for CO₂ activation as it favours electron transfer processes.

2.2. Adsorption of CO₂ on Perfect and Defective FeS (00) Surface

First, we investigated the adsorption of CO₂ on the perfect and defective FeS (001) surfaces by optimizing different initial adsorption configurations of the CO₂ molecule. At the perfect FeS (001) surface, the CO₂ molecule is physisorbed (Figure 3a) releasing an adsorption energy −0.20 eV. The layer of S⁻² ions terminating the surface shields the inner Fe⁺² ions, giving rise to repulsive interactions between the oxygen atoms in the CO₂ molecule and the surface S⁻² ions [29]. The adsorption energy (E_ads), selected structural parameters, charge transfers (Δq(CO₂)), and vibrational frequencies are summarized in

Table 1. Adsorption energies, structure, charges, and vibrational frequencies of CO₂ on perfect and defective (V_S1 and V_S2) FeS (001).

	Eads	C–Oa	C–Ob	C–Fe	Oa–Fe	α(OaCOb)	Δq(CO2)	υas	υs	υb
	(eV)	(Å)	(Å)	(Å)	(Å)	(o)	(|e|)	(cm−1)	(cm−1)	(cm−1)
Free CO2	-	1.176	1.176	-	-	180.0	0.00	2373	1323	631
(001)-P	−0.20	1.177	1.176	4.013	4.222	179.6	0.00	2349	1319	628
(001)-VS1	−1.78	1.417	1.230	2.044	1.913	118.2	−1.14	1547	809	672
(001)-VS2	−1.83	1.429	1.229	2.024	1.910	118.1	−1.23	1546	801	668

P = perfect stoichiometric surface, V_S1 = defective surface with single S-vacancy and V_S2 = defective surface with double S-vacancy.

. The CO₂ remained linear with no change in the C–O bond distances relative to the gas phase molecule, which is consistent with no charge transfer to the CO₂ molecule upon adsorption. Thus, we can conclude that the perfect FeS (001) does not activate the CO₂ molecule.

The CO₂ molecule on the other hand, chemisorbs strongly at the defective surfaces, binding to the exposed Fe sites (Figure 3b,c). The adsorption energy for CO₂ at the defective V_S1 and V_S2 surfaces is calculated at −1.78 and −1.83 eV, respectively. The CO₂ molecule gained a significant amount of charge from the interacting surface Fe ions, resulting in the formation of negatively charged bent species (CO₂^δ−) with elongated C–O bond distances (see Table 1). Bader population analysis shows that charges of 1.14 and 1.23 e⁻ were transferred from the interacting surface Fe ions to the CO₂ molecule at the defective V_S1 and V_S2 surfaces, respectively. The amount of charge gained is consistent with the lower work function predicted for the defective surfaces compared to the perfect stoichiometric surface. At the defective V_S1 surface, the C–O bonds are predicted to be 1.417 Å for C–O_a and 1.230 Å for C–O_b, compared to 1.174 Å in the gas phase. Similarly, at the defective V_S2 surface, the C–O bonds are elongated to 1.429 for the C–O_a and 1.229 Å for the C–O_b bonds. The α(O_aCO_b) bond angle for the adsorbed CO₂ molecule at the defective V_S1 and V_S2 surfaces is calculated at 118.2° and 118.1°, respectively, compared to the gas-phase linear (180.0°) molecule. The elongated C–O bonds are confirmed by vibrational frequency analyses, presented in Table 1. For instance, the asymmetric (υ_as) and symmetric (υ_s) C–O stretching and the bending (υ_b) vibrational modes for the adsorbed CO₂ at the defective V_S1 surface are calculated at 1547, 809, and 672 cm⁻¹, respectively. Compared to the gas phase C–O stretching modes (2373 and 1323 cm⁻¹), it is worth noting that both the υ_as and υ_s modes of the adsorbed CO₂ molecule are significantly red-shifted, a property which could be tested experimentally. These results are consistent with previous DFT studies which showed that the existence of oxygen vacancy defects on metal oxide ZnO (0001) [31,32,33], In₂O₃ (110) [34], CeO₂ (110) [35] and β-Ga₂O₃ (100) [42] surfaces promote CO₂ adsorption and activation.

Insight into the charge density redistribution in the CO₂–FeS adsorption system at the defective FeS (001) surfaces, was ascertained from the differential charge density (Δρ) iso-surface plot, which is calculated as:

Δρ = ρ(FeS_surf + CO₂) − (ρ(FeS_surf) + ρ(CO₂))

(2)

where ρ(FeS_surf + CO₂), ρ(FeS_surf) and ρ(CO₂) denotes the electron density of the FeS–CO₂ systems, the naked defective FeS (001) surface and the isolated CO₂, respectively, with the atomic positions kept the same as those of the total FeS–CO₂ system. The electron density difference isosurface plot, shown in Figure 4 reveals electron density redistribution within the CO₂–FeS system at the defective V_S1 surface. A significant electron density accumulation within the newly formed C–Fe and O–Fe bonds is observed, which is consistent with chemisorption.

2.3. Adsorption and Dissociation of H₂ on Perfect and Defective FeS (00) Surface

H₂ adsorption and dissociation will potentially impact CO₂ binding via competitive adsorption at the active sites, and subsequently influence CO₂ reduction via hydrogenation reactions. Therefore, we have systematically investigated the structures and energetics of molecular and dissociative H₂ adsorption on the perfect and defective FeS (001) surfaces. At the perfect FeS (001) surface, molecular H₂ is physisorbed (E_ads = −0.12 eV) with the closest hydrogen to surface distance (H–S) predicted at 3.289 Å (see Figure 5a). The H–H bond distance remained 0.75 Å, as in the gas phase. The dissociative adsorption of H₂ is found to be endothermic (E_ads = +1.54 eV) at the perfect FeS (001) surface (Figure 5b), with the dissociated H atoms binding to the S sites at an average S–H distance of 1.375 Å. The large energy difference between the molecular and the dissociative H₂ adsorption suggests that H₂ will preferentially desorb molecularly from the perfect FeS (001) surface rather than dissociate.

The adsorption of an H₂ molecule at the defective V_S1 FeS (001) surface resulted in spontaneous exothermic dissociation, as shown in Figure 6. The dissociative adsorption structures are labelled according to the sites to which the H atoms bind. For instance, D(Fe–Fe) corresponds to H atoms binding to two Fe sites, and D(Fe–S) corresponds to the H atoms binding at Fe and S sites. The adsorption energies of the D(Fe–Fe), D(Fe–S), and D(S–S) structures at the defective surface are calculated at −1.76 eV, −0.48 eV, and +1.46 eV, respectively. The exothermic adsorption energies calculated for the D(Fe–Fe) and D(Fe–S) structures, compared to the endothermic adsorption energy of the D(S–S) structure on the perfect surface, show that the exposed Fe sites due to the creation of an S vacancy are the active sites for H₂ dissociation. The newly formed S–H bonds have the same bond length (1.374 Å), whereas the Fe–H bonds are slightly different in each adsorption configuration, ranging from 1.585–1.839 Å, as shown in

Table 2. Adsorption energies and structure of H₂ on perfect and defective FeS(001), where M and D denotes molecular and dissociative adsorption modes, respectively.

	Configuration	Eads (eV)	Ha–Hb (Å)	Ha–S (Å)	Hb–S (Å)	Ha–Fe (Å)	Hb–Fe (Å)
(001)-perf	M(phys)	−0.12	0.750	3.297	3.289	-	-
-	D(S-S)	+1.54	3.653	1.376	1.375	-	-
(001)-VS1	D(Fe-Fe)	−1.56	1.838	-	-	1.635	1.637
-	D(Fe-S)	−0.48	4.229	-	1.374	1.839	-
-	D(S-S)	+1.46	3.658	1.375	1.374	-	-

perf = perfect stoichiometric surface, V_S1 = defective surface with single S-vacancy.

.

2.4. Coadsorption and Reaction of CO₂ and H₂ on Defective FeS (00) Surface

The reduction of CO₂ is one route for its removal and sequestration or utilization, as the reduced species may prove to be a valuable feedstock for other processes. Generally, CO₂ reduction may occur through the reverse water-gas shift reaction (CO₂ + H₂ → CO + H₂O), a major side reaction in methanol synthesis. Therefore, we have analyzed the reactions between CO₂ and H adatoms and determined the reaction pathway on the defective V_S2 FeS (001) surface. Different possible CO₂ and H₂ coadsorption structures were explored to determine the lowest-energy coadsorption structure (IS, Figure 7). In the coadsorption structure, the adsorption structures of CO₂ and H adatoms are identical to that of the isolated individual species. The coadsorption energy with respect to free CO₂ and H₂ was calculated at −3.68 eV, which is more exothermic than the sum of the isolated adsorption energies of CO₂ and 2H (−3.34 eV). The increased stability of the coadsorbed system indicates an attractive interaction between the CO₂ and H adatoms on the defective FeS (001) surface. By allowing the H adatoms to react with the adsorbed CO₂, we have calculated the reaction energies and activation barriers of the hydrogenation reactions of the activated CO₂ molecule on the defective FeS (001) surface.

Attachment of the first hydrogen to the carbon centre to form the formate (HCOO) intermediate structure was not possible as the carbon centre is fully saturated, preventing the attachment of hydrogen atoms. Hydrogenation of the surface-bound O_a ion first (COOH) results in the dissociation of the C–O bond to form CO and OH species adsorbed at adjacent bridging Fe sites (Figure 7). This reaction step is exothermic by 0.20 eV with an activation energy barrier of 1.43 eV. The formation of a COOH intermediate by hydrogenating the unbound O_b ion first was found to be less exothermic (E_R1 = −0.08 eV) with a higher activation barrier of 1.65 eV, suggesting that hydrogenating the surface bound O_a ion first represents a more favourable route. In the transition state (TS1), the distance between the surface bound O_a ion and the approaching H atom (O–H) was calculated at1.498 Å, whereas the Fe–H bond was calculated at 1.462 Å. The dissociating C–O bond in the TS1 structure was calculated at 1.895 Å compared to 2.542 Å in the intermediate structure (IntS: CO* + OH* + H*). The hydrogen atom of the formed OH species in the intermediate state is oriented towards the O atom of the CO fragment at a distance of 1.841 Å.

The second hydrogenation step involves the diffusion of the second H atom towards the formed OH species from the first hydrogenation step to form CO and H₂O. The formation CO and H₂O is found to be an exothermic process with a reaction energy of −0.47 eV and a lower activation energy barrier of 1.07 eV for the second hydrogenation step, compared to 1.43 eV for the first step. The higher barrier of the first hydrogenation step of CO₂ to form CO and OH suggest that it is rate-limiting step in the reduction of CO₂ to CO and H₂O on the defective FeS (001) surface. In the second transition state (TS2), the distance between the approaching H atom and OH species to form water was calculated at 1.798 Å, whereas the Fe–H bond was calculated at 1.568 Å. The resulting water molecule interacts with an iron site through the oxygen atom at an Fe–O distance of 1.954 Å, with one of its hydrogen atoms oriented towards the O atom of the CO fragment at a distance of 1.481 Å.

Attaching the second hydrogen atom to the surface bound O_a ion of the less stable COOH intermediate from the first hydrogenation step, results in the formation of dihydroxymethylidene (C(OH)₂), but this reaction is found to be highly endothermic (E_R2 = +1.41 eV) with an activation barrier of 1.33 eV. The defective FeS (001), therefore, favours CO and H₂O formation over a dihydroxymethylidene intermediate or formic acid. The isolated H₂O molecule binds to the defective surface with an adsorption energy of −0.48 eV. The smaller adsorption energy of an isolated H₂O compared to isolated CO₂ (−1.83 eV for the double S-vacancy surface), suggests that when both species are present in the gas phase or solution, the CO₂ molecules will outcompete H₂O molecules for the active sulfur vacancy sites.

3. Computational Details

The DFT calculations were carried out within the VASP-Vienna ab initio simulation package (VASP) [46,47,48,49], using the projector augmented wave (PAW) potentials [50]. The electronic exchange-correlation potential was calculated using the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) functional [51,52]. Non-local dispersion forces were accounted for using the Grimme DFT-D3 scheme [53]. An energy cut-off of 500 eV was used the expand the electronic wave functions and is found to be sufficient enough to converge the total energy and the Hellman–Feynman forces to within 10⁻⁶ eV and 0.001 eVÅ⁻¹, respectively. A Monkhorst-Pack [54] k-point mesh of 11 × 11 × 11 and 5 × 5 × 1 were used to sample Brillouin zone of the bulk and (001) surface, respectively.

A (3 × 3) simulation supercell of the (001) surface was employed for the adsorption calculations, which was large enough to minimize the lateral interactions between the adsorbate molecules in periodic cells. It consists of 36 iron atoms and 36 sulfur atoms, distributed in two FeS layers. To avoid interactions between periodic slabs, a vacuum size of 15 Å was added in the z-direction. The adsorption energies of adsorbates were defined according to equation:

E_ads (M) = E_M+surf − (E_surf + E_M)

(3)

where E_M+surf, E_surf, and E_M, respectively, represent the total energy of the adsorbate-surface systems, of the naked perfect or defective FeS (001) surface with sulfur vacancies, and of the adsorbates (CO₂ and H₂). The coadsorption energies were computed with respect to the sum of the total energies of the free molecules. Negative and positive adsorption energies indicate exothermic and endothermic processes, respectively. All energies (adsorption, reaction, and activation barriers) are corrected by the zero-point energy (ZPE), calculated according to the relation:

$$\Delta \mathrm{ZPE}={\left({\displaystyle \sum _{i=1}^{3n}\frac{h{v}_i}{2}}\right)}_{surf}-{\left({\displaystyle \sum _{i=1}^{3n}\frac{h{v}_i}{2}}\right)}_{gas}$$

(4)

where h and ν_i denote the Planck constant and vibrational frequencies, respectively. Charge transfer between the surfaces and adsorbates is quantified using the Bader charge analysis code developed by Henkelman and co-workers [55]. Transition states (TS) between reactants (initial, IS) and products (final, FS) were determined using the climbing-image nudged elastic band (CI-NEB) method [56,57]. Vibrational frequency calculations were performed to confirm an identified transition state, which is characterised by only one imaginary frequency. The activation energy barriers (E_A) are calculated as energy difference between TS and IS, whereas the reaction energies (E_R) were calculated as the energy difference between FS and IS.

4. Conclusions

In summary, a comprehensive analyses of the structural geometries, electronic properties, and the reaction mechanisms associated with CO₂ adsorption and its hydrogenation reactions on the sulfur-deficient FeS (001) surface was performed, using dispersion-correction density functional theory calculations. It is demonstrated that the presence of sulfur vacancies promotes CO_2, and H₂ adsorption and activation on the FeS (001) basal plane. The exposed Fe sites resulting from the creation of sulfur vacancies provide highly active catalytic sites for CO₂ activation and reduction. Accordingly, a high sulfur vacancy density is expected to improve the catalytic activity of FeS-containing catalysts for CO₂ activation and conversion. The results presented also show that CO₂ will outcompete H₂O for the sulfur vacancy sites when both species are present on the surface. The overall reaction energy for the hydrogenation of CO₂ to CO and H₂O is found to be exothermic by 0.47 eV, with the activation barrier for the rate-determining step calculated at 1.43 eV. The chemical picture and molecular-level insights derived from this work suggest that defective FeS (001) surfaces containing sulfur vacancies can catalyze the necessary reactions of CO₂ and H₂ to form small organic molecules.