All energy and electronic structure calculations were performed at the DFT level using the Vienna ab initio simulation package (VASP)^[25−26]. The electron-ion interaction was described with the projector-augmented wave (PAW) method^[27−28]. Electron exchange and correlation energies were treated within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA-PBE)^[29]. The cutoff energy was 400 eV and the Brillouin zone was sampled using the Monkhorst-Pack scheme^[30]. Energies and atomic forces were converged to within 10^–5 eV and 0.05 eV/Å, respectively. The Methfessel-Paxton method^[31] was chosen as the smearing algorithm. To give an exact description of the magnetic properties of iron and manganese, spin polarization was considered.

The MHM^[32−33] was employed to sample the configurational space of the surfaces and their reconstructions in the presence of adsorbate species. The MHM is a structure prediction method that relies on consecutive short molecular dynamics (MD) trajectories to escape from the catchment basin of the corresponding local minima, followed by geometry optimizations. The search aims at exploring low-energy structures by relying on a feedback mechanism based on the kinetic energy of the MD escape trials. The initial MD velocities are preferably aligned along soft mode directions of a local minimum, leading to crossing over low-energy barriers. Hence, the MHM utilizes the Bell-Evans-Polanyi principle^[34], which states that strongly exothermic reactions have overall low activation barriers. Recently, the MHM has been successfully employed for the exploration of low-dimensional systems, including 2D materials^[35] and surface reconstructions^[36−37].

All structural searches in this work were performed at the DFT level by coupling the MHM with VASP. The MHM was employed with an adaptive temperature ladder spanning 700–6700 K, allowing the system to overcome high-energy barriers and escape from local minima. The potential energy surface (PES) was initially screened with the MHM using coarse DFT settings with a cutoff energy of 270 eV. The searched low-energy structures were sorted according to DFT energy calculations with a cutoff energy of 300 eV.

To reflect the adsorption intensity and system stability of carbon, the average carbon binding energy (Fe + N_C → N_C/Fe) was calculated according to the following equation:where $ {{E}}_{\rm{N_C/slab}} $ is the total energy of the slab with N adsorbed C atoms in its equilibrium geometry, $ {{E}}_{\text{slab}} $ is the total energy of the bare slab, and $ {{E}}_{\text{C}} $ is the total energy of a free carbon in the gas phase. The adsorption energy of CO was calculated according to:where $ {{E}}_{\text{CO/slab}} $ and $ {{E}}_{\text{CO}} $ are the total energies of the slab with an adsorbed CO molecule and a CO molecule in the gas phase, respectively. The barrier for CO dissociation was calculated according to:where $ {{E}}_{\text{TS}} $ and $ {{E}}_{\text{IS}} $ are the total energies of the corresponding transition state and initial state, respectively. Transition states were located using the climbing-image nudged elastic band method^[38]. Zero-point energy corrections were included in the energetic data.

A 2 × 2 × 2 supercell containing eight body-centered cubic (BCC) Fe unit cells was built. Further, the Fe₁₅Mn₁ bulk was built by replacing one Fe atom with one Mn atom in the supercell (Figure S1), with equilibrium lattice constants of a = b = c = 5.733 Å. The bare surfaces were simulated with six-layer p(U × V) periodic slab models, in which p(2 × 2) for the (100) facet, p(2 × 1) for the (110) facet, and p(1 × 1) for the (111), (211) and (310) facets were adopted. The lattice parameters of the slab models are as follows: a = b = 11.47 Å for the (100) facet; a = 11.47 Å and b = 8.11 Å for the (110) facet; a = b = 8.11 Å and γ = 60° for the (111) facet; a = 9.93 Å and b = 8.11 Å for the (211) facet; and a = 5.73 Å and b = 18.13 Å for the (310) facet. The k-points of 2 × 2 × 1, 2 × 3 × 1, 4 × 4 × 1, 3 × 3 × 1 and 4 × 1 × 1 were set for the (100), (110), (111), (211) and (310) facets, respectively. The top four layers of these models were fully relaxed to explore the potential surface reconstruction caused by carbon deposition. When performing the structural search, considering the computational cost associated with large systems, slab models were constructed by retaining the top three layers of the obtained local minima geometries, and the top two layers were fully relaxed. The adsorption and dissociation of CO were calculated using the (100)-a and (310) surface structures generated through MHM searches, with the top two layers fully relaxed.

During the carbon deposition process on Fe-Mn surfaces, the position of the C atoms adsorbed/permeating on/into the slab model was determined by the respective atomic layers they entered/presented. To address the complexity introduced by the geometrical variety of surface structures, atomic layers are defined throughout the rest of the article. The atoms in the 1st, 2nd, and 3rd layers of each slab model were represented by different colors and denoted as 1L, 2L and 3L, respectively. Given the existence of certain “open” surfaces, there may be more than one sublayer within a defined layer. As carbon permeates into the interstices of the Fe-Mn alloy in atomic form, no particular carbon source is chosen in this work.

In accordance with our previous studies^[39−40], three low-index facets, namely (100), (110), and (111), and two facets with higher Miller indices, namely (211) and (310), were chosen to investigate carbon adsorption and permeation for Fe-Mn catalysts. For a comprehensive analysis, the models of the (100), (110) and (111) facets were categorized into two types, denoted as a and b for Fe-terminated and FeMn-terminated surfaces, respectively, to distinguish whether the Mn atom is a surface termination atom. The models of these surfaces are presented in Figure S2. The energetically stable surface structures and corresponding average carbon binding energies ($ {E}_{\text{B}}^{\text{ave}} $) at various carbon coverages during the carbon deposition process are shown in Figures S3–S10, which illustrate the sequential incorporation of carbon atoms into Fe-Mn surfaces.

Figure 1 illustrates the surface structures with the first C atom adsorbed on the 1L layer and permeating into the 2L and 3L layers on various surfaces. We first focus on the Fe-Mn(100)-a surface to illustrate the detailed mechanism of carbon adsorption and permeation. As shown in Figures 1 and S3, the first C atom adsorbs on the 4-fold hollow site in the 1L layer of the slab model, and additional C atoms adsorbed on the 1L layer continue to enter the 4-fold hollow site, mostly along diagonal directions relative to the first adsorbed C atom. The corresponding bonding geometry and carbon separation tendencies are consistent with both theoretical and experimental observations on the Fe(100) surface^[39,41]. The $ {E}_{\text{B}}^{\text{ave}} $ of the 14th C atom (10.6 C/nm²) becomes more stable in the 2L layer than in the 1L layer by approximately 0.03 eV, indicating that the C atom prefers to permeate into the 2L layer at this point. Proceeding in this way, it is found that the 2L layer can accommodate only one C atom before additional C atoms begin to permeate further into the 3L layer. For the Fe-Mn(100)-b surface with the Mn atoms exposed on the surface instead of located in the 2L layer, although it is also the 14th C atom that begins to enter the 2L layer, the 2L layer can accommodate up to 11 C atoms (18.3 C/nm²) before additional C atoms permeate further into the next layer. In contrast, the (111) facet exhibits significant differences. As shown in Figure 1, after introducing 8 C atoms into the first layer (14.1 C/nm²), the subsequent C atom prefers to directly enter the 3L layer instead of the 2L layer. For (110)-a, (110)-b, (211), and (310) surfaces, C atom permeation into the 2L layer occurs at 6.5, 8.6, 12.4 and 12.4 C/nm², respectively, and permeation into the 3L layer occurs at 10.8, 11.8, 14.9 and 16.4 C/nm², respectively. It is found that as the carbon coverage increases, the (211) facet undergoes reconstruction, during which an exposed Mn atom and an adjacent Fe atom migrate to nearby sites (Figure S9). These results demonstrate that crystallographic orientation and surface termination jointly govern the carburization behavior of Fe-Mn alloy surfaces.

It is worth noting that, for these results, each C-deposited surface structure discussed above represents a local minimum. This is achieved by “shooting in” C atoms through specified sequential processes and allowing them to reach the optimum site via a single geometry optimization. Surface reconstructions can cause phase transformations and form carbide phases. To determine if there are reconstructed structures with lower energy on the PES that can be achieved by overcoming activation barriers, we conducted minima-hopping structural searches for C-deposited surface structures.

Herein, the Fe-Mn(100)-a and (310) surfaces at medium and high carbon coverages were selected to represent the low-index and high-index surfaces, respectively. After removing duplicate configurations, the number of candidate configurations generated by the MHM was 57 and 56 for the (100)-a surface at medium and high carbon coverages, respectively, and 49 and 49 for the (310) surface at medium and high carbon coverages, respectively. Geometry optimizations were then performed on the top ten most stable structures predicted by the MHM for each system to access the global minima. As shown in Figure 2(a), the surface structure of the (100)-a surface at medium carbon coverage (6.1 C/nm²) remains almost unchanged before and after applying the MHM. At high carbon coverage (11.4 C/nm²), carbon diffusion occurs after exploring the PES using the MHM. The atom A diffuses from the 4-fold hollow site (1L layer) above the Mn atom to another 4-fold hollow site (1L layer) above the Fe atom. The atom B in the interstitial site diffuses from the 2L layer with Mn atoms to the 3L layer without Mn atoms. These findings suggest that, for the (100) facet at increasing carbon coverages, C atoms may not tend to be located near Mn atoms. For the (310) surface at medium carbon coverage (8.7 C/nm²), the surface structure searched through the MHM does not undergo the rearrangement of atomic sites (Figure 2(b)). At high carbon coverage (16.4 C/nm²), for carbon atoms, the atoms A and B in the 2L layer diffuse to the 3L layer, and the atom C diffuses from the 2L layer to the junction of the 1L and 2L layers. This indicates that, for the (310) surface at high carbon coverages, the permeated C atoms are prone to diffuse, forming a structure with lower energy. No segregation of Mn atoms from the subsurface to the outermost layer was observed for either facet, indicating that Mn atoms doped into the Fe lattice remain in a stably separated state in our models.

To deepen our understanding of the structural differences between local minima and global minima, we investigated the radial pair distribution function (RDF) of these surface structures. The RDF provides the density probability of finding neighboring atoms at a given distance from a reference atom. For (100)-a and (310) surfaces at medium carbon coverages, as shown in Figure 3(a) and Figure 4(a), the shapes and positions of the RDF curves for local minima closely resemble those of the global minima obtained through the MHM. This is consistent with our previous discussion that the surface structure at medium carbon coverages remains nearly unchanged before and after applying the MHM (Figure 2). When a high carbon coverage is reached, the shapes and positions of their RDF curves remain similar, except for the highest peaks of the Fe-C and Fe-Mn curves (Figure 3(b) and Figure 4(b)). This discrepancy can be attributed to the diffusion of carbon atoms and the accompanying geometry relaxations of Fe and Mn atoms (Figure 2). Overall, our local minima, obtained by adding carbon atoms with a specified sequence, closely approximate the global minima geometries obtained from the MHM. For both (100)-a and (310) surfaces, although Fe and Mn atoms in the 1L and 2L layers undergo noticeable relaxations at high carbon coverages, no surface reconstructions occur, which indicates that the initial stage of carbon permeation is insufficient to induce a phase transition.

The $ {E}_{\text{B}}^{\text{ave}} $ represents the thermodynamic stability of C atoms deposited on the Fe-Mn surfaces. The more negative the value of $ {E}_{\text{B}}^{\text{ave}} $, the more energetically favorable the C-deposited surface structure is. Figure 5 illustrates the $ {E}_{\text{B}}^{\text{ave}} $ as a function of the number of introduced C atoms per unit area (N_C/nm²).

At the lowest carbon coverage, the $ {E}_{\text{B}}^{\text{ave}} $ of different Fe-Mn surfaces follows the sequence (100)-a < (100)-b < (110)-b < (310) < (110)-a < (111)-a < (211) < (111)-b, indicating that carbon adsorption on the (100) facet is initially most favorable compared with other surfaces.

With increasing carbon coverage, the investigated surfaces can be classified into three types according to the evolution in their $ {E}_{\text{B}}^{\text{ave}} $ before C atoms permeate into the 2L layer. The $ {E}_{\text{B}}^{\text{ave}} $ of (100)-a, (100)-b, and (310) surfaces show an initial gentle and subsequent rapid rise, while the $ {E}_{\text{B}}^{\text{ave}} $ of (110)-a and (110)-b surfaces experience a drastic increase. Conversely, the $ {E}_{\text{B}}^{\text{ave}} $ of (111)-a, (111)-b, and (211) surfaces gradually decrease with oscillation. We observe a strong correlation between the diverse trends in the thermodynamic stability of carbon adsorption on the 1L layer of these three types of surfaces and the degree of surface undercoordination ($ {{d}}_{\text{suc}} $) of these surface structures. This descriptor quantifies the surface undercoordination of arbitrary metal surfaces by measuring the disruption of intermolecular bonds^[40,42]. The (100)-a, (100)-b and (310) surfaces with high $ {{d}}_{\text{suc}} $ values of 0.034, 0.033 and 0.044, respectively, show a higher affinity for C atoms compared with other surfaces. Until the carbon coverage reaches approximately 16 C/nm², C atoms preferentially deposit on (100)-a, (100)-b, and (310) surfaces. The (111)-a, (111)-b and (211) surfaces with low $ {{d}}_{\text{suc}} $ values of 0.013, 0.011 and 0.009, respectively, display poor thermodynamic stability of carbon adsorption due to their high coordination. The decreasing trend in $ {E}_{\text{B}}^{\text{ave}} $ as the carbon coverage increases indicates that the binding strength of carbon on these surfaces strengthens with the introduction of more C atoms into the system. The (110)-a and (110)-b surfaces, which exhibit stability between the two types of surfaces mentioned above, are characterized by moderate $ {{d}}_{\text{suc}} $ values of 0.026 and 0.027, respectively, compared with other surfaces.

As C atoms begin to permeate into the 3L layer, the number of introduced C atoms per unit area follows the order (100)-a ≈ (110)-a ≈ (110)-b < (211) ≈ (111)-a ≈ (111)-b ≈ (310) < (100)-b, and the corresponding $ {E}_{\text{B}}^{\text{ave}} $ follows the order (100)-a < (310) ≈ (211) ≈ (111)-a ≈ (111)-b ≈ (100)-b < (110)-a ≈ (110)-b (Figure 5). Consequently, from the perspective of thermodynamics, the (100)-a surface exhibits the highest stability for C atoms permeating into the 3L layer at lower carbon coverage compared with other surfaces. Additionally, the (100)-b surface shows the most significant capability to accommodate C atoms in the 1L and 2L layers, while the (110)-a and (110)-b surfaces possess the strongest resistance to the permeation of additional C atoms into the 3L layer.

Several studies have reported that the formation of a solid solution compound over the coprecipitated Fe-Mn catalysts could influence the reduction and carburization behaviors of Fe catalysts^[12−13,43]. Furthermore, the manganese promoter was found to enhance the catalyst’s resistance to deactivation induced by carbon deposition^[44]. These studies uncover the significant impact of the manganese promoter on the carburization of iron-based catalysts. To understand the impact of manganese on the initial stage of carbon permeation during the carburization process from an atomic-scale perspective, we compared the calculated results for carbon permeation on Fe-Mn surfaces with those on pure Fe surfaces reported in our earlier study^[39]. As illustrated in Figure 6(a), when the first C atom permeates into 2L and 3L layers, the $ {E}_{\text{B}}^{\text{ave}} $ of various Fe-Mn surfaces (except for (100)-a surface with the first C atom permeating into the 3L layer) are higher than those of Fe surfaces with the corresponding Miller indices. This indicates that the thermodynamic stability of carbon permeation on Fe-Mn surfaces is lower than that on Fe surfaces. Additionally, as shown in Figure 6(b), the number of C atoms that can be introduced per unit area on Fe-Mn surfaces (except for (100)-a and (110)-a surfaces with the first C atom permeating into the 3L layer) is similar to or higher than that on Fe surfaces. Thus, in comparison to Fe surfaces, Fe-Mn surfaces may accommodate more C atoms before the C atoms permeate into the bulk region. These results suggest that the addition of manganese is thermodynamically unfavorable for carbon permeation on iron surfaces.

The binding strength between metal surfaces and adsorbates essentially depends on the electronic structure of the surface. To elucidate how carbon deposition and permeation regulate the electronic properties of Fe-Mn surfaces, the projected density of states (PDOS) of different C-deposited Fe-Mn surfaces were calculated in this work. The (100) facet is selected as a representative case to illustrate the electronic structure of C-deposited Fe-Mn surfaces. For the bare surface, as shown in Figure 7(a) and 7(d), the spin-up total density of states (TDOS) peaks at energy level around −4.0 to −2.5 eV and −1.5 to −0.9 eV appears clearly. Besides, there are obvious spin-down TDOS peaks at energy level around −2.3 to −0.6 eV and 1.3 to 2.5 eV. These four peaks are contributed by d electrons of Fe and Mn atoms, suggesting that there is considerable overlap between Fe 3d and Mn 3d orbitals. They collectively constitute the main part of the energy bands near the Fermi level. While there are wide and low-intensity regions at high energy level, which are composed of the overlap between Fe s-p and Mn s-p orbitals, mixing with Fe 3d and Mn 3d orbitals. Compared with Mn, the increased nuclear charge of Fe will be less effectively screened by the tight-bound 3d electrons. Therefore, the 3d bands of Fe, feeling more attraction from the nucleus, are slightly lower in energy than the 3d bands of Mn, as indicated by the low energy region of approximately −4.5 eV.

The permeation of C atoms into vacancies within the Fe-Mn alloy can be regarded as a perturbation to the alloy DOS. With carbon permeation, the four main peaks mentioned above are weakened to varying degrees but remain prominent, indicating that the electron density is picked up by C atoms from the d bands of metals due to the electronegativity difference between metals and carbon (Figure 7(b) and 7(e)).

The PDOS of Fe-d and Mn-d tend to shift to a lower energy level, indicating the potential downshift of the overall d bands. The spin-up TDOS at energy level around −7.5 to −5.0 eV and the spin-down TDOS at energy level around −7.5 to −3.0 eV intensify clearly. As shown in Figure 7(c) and 7(f), we find that these intensified states are mostly the 2p states of C, and the density of these states increases with the increase of the number of introduced C atoms. Moreover, there are isolated states of TDOS at energy level around −14.0 to −10.5 eV, which consists mainly of s and p electrons of C, with little contribution from Fe-d. The features of PDOS discussed above are basically the same as all Fe-Mn surfaces we investigated, so we do not elaborate on the other surfaces (Figures S11−S16). The difference in electronic structure between different surfaces is mainly reflected in the d-band center, thus we next focus on the d-band centers of these surfaces.

The d-band center is the average energy of the d-band and can serve as an effective descriptor of the ability for electron transfer between metal surfaces and adsorbates^[45−46]. Therefore, the d-band centers of the entire metal on different C-deposited surface structures were calculated to further elucidate the influence of the addition of carbon atoms on the catalytic-activity-related electronic properties. As shown in Figure 8, the d-band centers of different C-deposited surfaces move away from the Fermi level with the increase of the number of introduced C atoms (except for Fe-Mn (111)-b), which indicates that there is an interaction between the C s-p orbitals and metal 3d orbitals to modulate the d-band centers of Fe-Mn surfaces. Compared with Fe surfaces^[39], the d-band centers of Fe-Mn surfaces are closer to the Fermi level and exhibit a more gentle decline, suggesting that manganese incorporation retards the evolution of surface electronic structures toward those characteristic of bulk iron carbide phases.

To further elucidate how changes in surface geometry and electronic structure induced by varying carbon coverage affect catalytic behavior, we employed CO as a probe molecule and calculated its adsorption and dissociation on the (100)-a and (310) surface structures identified through MHM searches. These processes play a critical role in FTS^[5].

As shown in Figure 9(a) and 9(b), compared with the bare surfaces, CO adsorption energies on the C-deposited surfaces become more positive, indicating weakened adsorption. This behavior is consistent with the downward shift of the d-band center away from the Fermi level, which weakens interactions with the antibonding orbitals of the adsorbate and thus leads to reduced adsorption strength^[45−46]. In addition, the activation barriers for CO direct dissociation on C-deposited surfaces are higher than those on the corresponding bare surfaces (Figure 9(c) and 9(d)). These trends are particularly pronounced on the (310) facet, which provides a higher density of active sites; as a result, geometric effects become less dominant, allowing the influence of changes in electronic properties to be more clearly manifested.