Highly stable two-dimensional metal-carbon monolayer with interpenetrating honeycomb structures

With the ongoing effort in proposing and realizing functional two-dimensional (2D) materials, we predict by first-principles calculations a family of 2D metal-carbon (M–C) crystals consisting of M–C trigonal lattice interpenetrated with the metal buckled honeycomb structure. We suggest by simulations that the 2D M–C crystals can be readily fabricated by a self-organizing lattice reconstruction process after placing metal atoms on hollow sites of γ-graphyne. In total, we found 12 members of the family and they exhibit a variety of electronic and magnetic properties. In this work, we highlight and focus on the Fe member of the family, 2D-Fe2C12. Each Fe in 2D-Fe2C12 has a magnetic moment of 1 μB due to the spin splitting of Fe E1 bands at Fermi surface, resulting in half metallicity and high catalytic activity with unusually high-density single-atom Fe active sites. Ab initio molecular dynamics simulations revealed that the 2D-Fe2C12 retains its structural integrity up to 700 K of simulated short duration annealing. We expect these results to stimulate experimental research for the 2D M–C crystals we proposed.


INTRODUCTION
Since the rise of graphene 1,2 , two-dimensional (2D) materials have attracted enormous interest in recent years owing to their unique properties and great potential in various applications. Discovering 2D functional materials beyond graphene 3,4 that have exotic physical and chemical properties lies at the heart of the current research effort along this direction. Examples of these 2D materials include h-BN 5,6 , g-C 3 N 4 7,8 , and 2D transition-metal dichalcogenides 9,10 that are semiconductors with sizeable bandgaps, and 2D-Cu 2 Si monolayer 11,12 that is a metal with high conductivity. With the ongoing effort of looking for new functional 2D materials, 2D carbon-based networks, such as porphyrin and phthalocyanine systems 13,14 , have been widely considered as an ideal template for realizing emerging 2D properties because of the high stability and flexibility of the network. The reason behind the stability of such systems is owing to the presence of N atoms that not only form strong bonds with C atom but also interact favorably with metal atoms in porphyrin and phthalocyanine. On the other hand, interactions between metal and C atoms tend to be much weaker. The evidence for this is that doping of carbon-based materials like graphene and carbon nanotubes with metal atoms result in poor stability of the dopants, leading to sintering and formation of metal nanoclusters, as indicated from various theoretical studies and experimental observations [15][16][17][18] . As a result, vacancies in carbon substrate like graphene and/or other impurity atoms are necessary for experiments to stabilize the 2D carbon-metal systems [19][20][21] . It is then of special significance that we, here in this paper, propose a carbon framework involving metal-carbon (M-C) bonds, which can give rise to a family of highly stable 2D crystals. The 2D M-C crystal has a structure that has two metals and 12 C in one unit cell (with a chemical formula M 2 C 12 ) and consists of both M-C 5 and C 6 hexagonal rings. Metal atoms themselves form a periodic honeycomb pattern embedded in a graphene honeycomb lattice, which we call in this paper the Honeycomb-in-Honeycomb (HIH) structure.
With different choices of metal atoms, 2D-M 2 C 12 exhibits a variety of electronic and magnetic properties. We show that the 2D-M 2 C 12 can be metal or semiconductor, and magnetic or nonmagnetic. In particular, in this paper, we highlight and focus on the Fe member of the family, 2D-Fe 2 C 12 , which is a half-metal with 1 μ B of magnetic moment on each Fe. Great research efforts have been made in recent years in inducing magnetism in 2D materials [22][23][24][25] . Here, we show that the 2D M-C crystals, for example 2D-Fe 2 C 12 , could be a avenue for making 2D magnetic materials that have potential for nanoscale spintronic applications such as spin filter and transistors. We will further show that the extremely high density of exposed Fe atoms in 2D-Fe 2 C 12 can function as catalytic active centers. 2D carbon-based materials have attracted considerable attention as a host material for impregnating catalytically active metal atoms as single-atom catalysts (SACs) [19][20][21]26 . For 2D-Fe 2 C 12 , unlike most of previously proposed SACs, the catalytic activity is intrinsic and not by modifying local properties of the system. Also, sintering of the catalysts resulting from the agglomeration of the dopants subjected to the harsh operating conditions can be avoided, therefore, the stability of 2D-Fe 2 C 12 manifests as an important advantage over many other SACs.

RESULTS AND DISCUSSIONS Design of 2D-M 2 C 12
The conception of such material design is built upon the basis that graphyne (GY), a family of allotropes of 2D carbon consisting of a mixture of sp-and sp 2 -hybridized C atoms, offers various degrees of porosity 27,28 . γ-GY has been successfully synthesized via various bottom-up approaches 27,29 . Previous theoretical studies have shown that it is a good host material for metal atoms, with the triangular hollow sites shown in Fig. 1 as good-trapping sites 30 . The motivation of this work is that if the M-C interaction is sufficiently strong, a significant lattice reconstruction could occur after metal atoms being trapped, and then each highly reactive C-C triple bond can break one of its π bonds to facilitate M-C bond formation, resulting in a tri-coordinated metal atom and the formation of 2D-M 2 C 12 as shown in Fig. 1. To verify this, we performed density functional theory (DFT) based first-principles calculations. We first obtained the atomic and band structures of GY, as shown in Supplementary Figure 1 in supporting information, which agree very well with previous reports 28 . Then, each metal atom was placed at one of the triangular hollow sites and the ground-state structure of the monolayer was obtained using full relaxation without any symmetry constraints. We tested various metal species, for 12 of them (e.g., Ti, Cr, Mn, Fe, Cu, Mo, Ru, Rh, W, Re, Os, Ir), the carbon skeleton prefer undergoing a lattice reconstruction and the graphynic lattice is no longer visible, leading to a self-organized periodic honeycomb pattern of metal atoms in 2D-M 2 C 12 (Fig. 1). Note that in the lattice reconstruction process, there is no substitution of C atoms by metal atoms as evident in Supplementary Video 1.
Electronic structures of 2D-M 2 C 12 , M = Ti, Rh, Ru, Fe The 2D-M 2 C 12 monolayer shows a variety of electronic structures with different choices of metal atoms. Here, we show Ti, Rh, Ru, and Fe members as four typical cases. The resulting optimized atomic structures of all four systems resembles the graphenic network with a trigonal lattice of P-3 space group. Nevertheless, these materials exhibit very different electronic properties evident from the band structure of each system shown in Fig. 2, in which 2D-Ti 2 C 12 is a nonmagnetic metal, 2D-Rh 2 C 12 is a nonmagnetic semiconductor with a direct bandgap of 0.45 eV at M point, 2D-Ru 2 C 12 is a magnetic metal with a 0.31 μ B of magnetic moment on each Ru atom, and 2D-Fe 2 C 12 is half metal with a semiconducting spin-up channel of quasi bandgap 0.30 eV at M point, and~1 μ B of magnetic moment resides on each Fe atom. The highly tunable electronic structure of 2D-M 2 C 12 makes the material useful for various kinds of applications. In this work, we are particularly interested in half-metallic 2D-Fe 2 C 12 . Fig. 1 The possible mechanism for the synthesis of 2D-M 2 C 12 from GY. The triangular hollow sites in GY are good-trapping sites for the metal atoms. Upon adsorption, the strong metal-carbon interactions induce lattice reconstruction in a self-organizing process forming 2D-M 2 C 12 . Unit cells are enclosed in cyan dashed lines. Atomic structure of 2D-Fe 2 C 12 The atomic structure of the relaxed 2D-Fe 2 C 12 is shown in Fig. 3a, which is also very similar for the Ti-, Rh, and Ru-systems, forming a trigonal lattice with P-3 space group in the 2D plane. The optimized lattice constants are obtained by fitting using the third order Birch-Murnaghan equation of state 31 shown in Supplementary Figure 2a, giving a = b = 6.75 Å. There are two different types of C atoms, which we will denote as C 1 (bonded with Fe) and C 2 (the rest), respectively, although both are sp 2 hybridized. The calculated bond lengths for C 1 -C 1 , C 2 -C 2 , and C 1 -C 2 are 1.39 Å, 1.41 Å, and 1.43 Å, respectively. Each Fe is threefold-coordinated with C 1 atoms and the calculated C 1 −Fe bond length is 1.80 Å. The material shows an interesting HIH structure as aforementioned.
It is noteworthy that the configuration with two Fe atoms that protrude in the same direction (atomic structure shown in Supplementary Figure 2b) were also tested but is 10 meV higher in energy than the ground-state structure. Nevertheless, the electron localization function plot in Supplementary Figure 3 shows significant hybridization between C 1 and Fe atoms and these favorable interactions stabilizes the Fe atoms, which we will show shortly, and results in the magnetic moment largely localized on the Fe atoms, as supported from the spin differential density plotted in Fig. 3b.

Stability of 2D-Fe 2 C 12
The stability of 2D-Fe 2 C 12 was preliminarily tested by calculating the cohesive energy (E coh ) and the formation energy (E form ). The cohesive energy was estimated to be 6.86 eV atom −1 . For comparison, the E coh of other 2D materials were also calculated and tabulated in Table 1. The high value of E coh shows that the stability of 2D-Fe 2 C 12 is at least comparable to materials like silicene, 2D-Cu 2 Si and even monolayer MoS 2 . Furthermore, the formation energy of 2D-Fe 2 C 12 from GY and Fe atoms is calculated to be −0.56 eV atom −1 . Such exothermicity should be indicative of the facile formation of 2D-Fe 2 C 12 . The structural stability of 2D-Fe 2 C 12 was ascertained by the calculated phonon dispersion with all the vibrational modes being real in the entire Brillouin zone (BZ) as plotted in Fig. 3c. A noteworthy feature is that the phonon branches occurring above ν = 63.0 THz in GY 32 are absent because the C≡C were broken so that the C 1 atoms can hybridize with Fe states.  The temperature-dependent stability of the material was assessed by the ab initio molecular dynamics (AIMD) simulations. The corresponding trajectories and the snapshots of the geometries at different junctures are shown in Supplementary  Figure 4, with the variation of the Fe-C 1 bonds and the snapshots at 6 and 12 ps at 300 K shown in Fig. 3d. It is evident that the material is very stable, with the Fe atoms showing no signs of undergoing any form of significant migration. From a series of AIMD simulations (Supplementary Figure 4), the atomic network is well-maintained at the temperatures of up till 700 K, at which 2D-Fe 2 C 12 could survive a 12 ps annealing up to 700 K. At 800 K, the out-of-plane vibrations of Fe atoms are violent enough to break the Fe-C 1 bonds, resulting in the destruction of the network consisting of FeC 5 rings and the partial restoration of the GY carbon skeleton, indicating the decomposition of the material. It is worth mentioning here that the error of the AIMD calculations increases with temperature, but based on our results, it is safe to conclude that the material is highly stable under room temperature.
Magnetic coupling in 2D-Fe 2 C 12 We next study the magnetic coupling between the Fe atoms in calculating the total energies of the ferromagnetic (FM) and antiferromagnetic (AFM) configuration in one unit cell, and also two other AFM configurations in a 1 × 2 supercell (AFM2 and AFM3 in Supplementary Figure 5). It was found that the FM configuration is 88, 35, 13, and 103 meV lower in energy than the AFM, AFM2, AFM3 and the nonmagnetic configuration per unit cell, respectively, indicating that the Fe atoms prefer to couple ferromagnetically. The calculated magnetic moment for the FM configuration is~2 μ B per unit cell, with 1 μ B residing on each Fe atom, as shown from the isosurface plot of the spin density.
Origin of magnetism in 2D-Fe 2 C 12 The magnetism of 2D-Fe 2 C 12 can be elucidated from its band structure. The energy bands responsible for the magnetism is shown in Fig. 4a. Those bands show an intriguing half-metallic behavior as aforementioned. Half metals have great potential in spintronics applications such as spin filtering and transistors that have been extensively reported in literature [33][34][35] . The half-metallic behavior of 2D-Fe 2 C 12 provides more opportunities for future design of 2D spintronic devices. In the absence of relatively flat impurity bands in the BZ, the bands in the vicinity of Fermi energy (E F ) are largely dominated by those originated from the hybridization between the Fe d orbitals and the π bands of the carbon skeletal network. Despite this strong hybridization, the calculated magnetic moment would originate from the occupation of the largely non-bonding energy bands 30 , which we have identified and labeled by their symmetry as A, E1, and E2 in Fig. 4a. The symmetry of bands is identified by plotting the banddecomposed charge densities (Fig. 4b). There are five such nonbonding bands (two E2, two E1, and one A) in each spin channel near E F in the range of -1.75 to −0.20 eV for the spin-up channel and -1.00 to 0.50 eV for spin-down. The bands A is the highest in energy, followed by E1 and E2. Such order in energy can be understood by d orbital splitting and hybridization of threefold coordinated Fe atoms in a C 3v ligand field (see Supplementary Notes in supporting information for more detailed analysis).
To work out the occupation of those bands, we noted that for the two Fe with 16 valence electrons in consideration (with each Fe contributing six 3d and two 4 s electrons), six of these electrons are required to form σ bonds with the C 1 atoms and two more to maintain the π band. The remaining eight electrons would fill up the largely non-bonding bands mentioned A, E1, and E2. All five of them in the spin-up channel would be filled, with the quasivalence band (VB) mainly composed of the Fe d z 2 , as corroborated from the projected density of states (PDOS) and the banddecomposed density. The remaining three electrons would fill up the bands E1 and E2 in the spin-down channel, with bands E1 crossing E F and approximately half filled. As such, this will leave two unpaired electrons giving rise to the magnetic moment of 2 μ B per unit cell. The results can be summarized in a simple schematic of non-bonding orbitals (shown in Fig. 4c) that represents the largely non-bonding bands, with the order of the orbitals obtained from the energies of the non-bonding bands at Γ point. Note that the spin-down E1 bands cross E F so we put the two spin-down E1 orbitals right at E F and make them half filled, which finally lead to 2 μ B of magnetic moment in one unit cell. This non-bonding orbital model is also capable of deducing the Fig. 4 Electronic structure analyses for 2D-Fe 2 C 12 . a The spin-polarized band structure of 2D-Fe 2 C 12 . Fe bands near Fermi are labeled by their symmetry A, E1, and E2. b Band-decomposed charge density of the largely non-bonding bands at Γ point, which are responsible for the origin of the magnetic moment. c A simple schematic of non-bonding orbitals that represents the largely non-bonding bands A, E1, and E2. Note that the spin-down E1 orbitals are half filled. magnetic moments calculated in other members of the family. The band structure and the spin density suggest that the magnetic Fe centers can be excellent electron donors. The halffilled conduction bands E1 of the spin-down channel crosses the E F , which are ready to transfer electron to an acceptor, suggesting that besides spintronics applications, the material may also be catalytically active.

Catalytic activity of 2D-Fe 2 C 12
To study the catalytic activity of 2D-Fe 2 C 12 , the CO oxidation reaction (COR) was chosen as a probe reaction. The adsorption of CO, O 2 , and CO 2 is crucial for understanding the catalytic activity. It was found that O 2 and CO molecules do not adsorb on the C atoms but can chemisorb strongly on Fe with adsorption energies being −1.32 and −0.95 eV, respectively, with the most stable adsorption configurations shown in Fig. 5a (O 2 ) and S8 (CO). O 2 takes a side-on adsorption on Fe, with a sizeable charge transfer of 0.78 |e| from Fe to the 2π* antibonding orbital of O 2 as indicated by the charge redistribution iso-surfaces (Fig. 5a). In Fig. 5b, we show the energy levels (near E F ) of the Fe and O 2 before adsorption. As suggested by the figure, the <1 |e| charge transfer shall happen from the partially filled Fe E1 spin-down orbital to the empty spin-down 2π* orbital of O 2 . Such charge transfer is verified by the PDOS of O 2 before and after adsorption plotted in Fig. 5c. Before adsorption, the spin-down 2π* orbital of isolated O 2 is unoccupied. After adsorption, the spin-down 2π* orbital split to two peaks, with one of them is pulled below Fermi level, indicating the partial charge transfer onto the 2π* orbital. The charge transfer decreases the magnetic moment of O 2 from 2.0 to 0.55 μ B and leads the O-O bond significantly elongated from 1.23 to 1.40 Å. Compared with O 2 , CO-adsorption is weaker with a smaller 0.29 |e| transferred from the Fe to CO. At last, CO 2 does not adsorb on any C sites and only recorded a very weak adsorption energy of −0.03 eV on Fe (see adsorption structure in Supplementary Figure 8b), implying that the CO 2 molecules can spontaneously desorb from Fe under room temperature.
We also considered the adsorption of CO with one O 2 molecule pre-adsorbed on Fe. We found that in this case, the adsorption energy of CO is~−0.32 eV. The co-adsorption configuration with the charge redistribution induced by both CO and O 2 is shown in Fig. 5d. The analysis of CO, O 2 , and CO 2 adsorptions suggests that the catalyzed COR can happen in both the Eley-Rideal (ER) and the Langmuir-Hinshelwood (LH) mechanisms. Indeed, our calculations show that both mechanisms are possible, whereas the LH mechanism is favored with a low reaction barrier <0.5 eV. For the LH mechanism, as shown in Fig. 6a, in the calculated first step of the reaction (CO + O 2 → CO 2 + O), the reaction proceeds initially with the rotation of the adsorbed O 2 molecule to form the peroxotype (OCOO) metastable transition state (TS), with a low reaction barrier (E b ) of 0.44 eV. Then the O-O bond increases, requiring only 0.38 eV to form the second TS. Subsequently, the Fe-C bond breaks, releasing a CO 2 molecule. This process is highly exothermic with a Δ r E = −3.29 eV. For the second step (O + CO → CO 2 , shown in Fig. 6b), another CO comes, reacts with the left O atom and forms CO 2 . This step has a reaction barrier of 0.48 eV. For the ER mechanism, the calculated energy profile is shown in Supplementary Figure 9, we noted that the gaseous CO molecule attacked the pre-adsorbed and activated O 2 , reached the TS with E b = 0.79 eV. Then, the second CO molecule attacks the left O to form the second CO 2 with E b = 0.48 eV. Therefore, 2D-Fe 2 C 12 is an excellent catalyst for COR which will predominantly happen via the LH path with a low reaction barrier of 0.48 eV, indicating that the reaction can happen under room temperature.
Via first-principles calculations, we proposed a family of 2D M-C monolayers with an intriguing HIH structure in which metal atoms form a periodic honeycomb pattern embedded in graphene honeycomb lattice. The 2D M-C monolayer has a chemical formula M 2 C 12 and can be obtained by a self-organizing lattice reconstruction process of GY with metal atoms in its hollow triangular sites. In total, we have found 12 members of the family. With different metal, 2D-M 2 C 12 show different electronic and magnetic properties, which makes the material useful for various applications. We focused our analysis on the Fe member of the family, 2D-Fe 2 C 12 , which is a half metal. The structural stability of 2D-Fe 2 C 12 were evaluated and we found that it can survive up till 700 K. Detailed analysis reveals that the half-metallic behavior originates from the spin-split of Fe E1 orbital at E F . The half metallicity clearly suggests that the 2D-Fe 2 C 12 can be useful for spintronic applications such as spin filtering and spin transistors. We also show that the exposed magnetic Fe sites in 2D-Fe 2 C 12 are highly catalytically active towards CO oxidation. The ultra-high  density of Fe atoms implies high performance of 2D-Fe 2 C 12 as a single-atom site catalyst. We foresee that such metal embedded carbon-based material can also be realized soon with the addition of the appropriate reagents in the reaction mixture and its application potential are awaiting to be verified by experiments.

Structure optimization and electronic structure analysis
All first-principles calculations were carried out under the spin-polarized DFT formalism. Structure optimization, electronic structure, and total energy calculations were done using the Vienna Ab initio Simulation Package (VASP) 36,37 employing the projector augmented wave pseudopotentials 38 and a plane wave basis set with a kinetic cutoff energy of 500 eV were used. In addition, the Perdew-Burke-Ernzerhof format 39 of the generalized gradient approximation was adopted for the exchange-correlation functional. All the structures were relaxed without any symmetry constraints and until the total energy and the residual Hellmann-Feynman force acting on each atom is <10 −5 eV and 0.01 eV Å −1 , respectively. To eliminate unphysical interaction between the replicas of periodic images, a vacuum space of at least 20 Å in the direction perpendicular to the 2D-Fe 2 C 12 atomic plane, i.e., the z-direction. The BZ was sampled by a 5 × 5 × 1 and 7 × 7 × 1 Γ-centered k-point grid in the Monkhorst-Pack scheme for the structural optimization and the electronic structure calculations, respectively. The effects of van der Waals' interactions were considered through the so-called DFT + D2 method 40 and we found that DFT + D2 gives essentially the same results as DFT.

Phonon and AIMD
Structural stability of the material is attested by calculating its phonon band dispersions using the density functional perturbation theory implemented in VASP and the Phonopy 41 , in which a 3 × 3 × 1 supercell was used sampled with a 2 × 2 × 1 k-mesh. The spin-polarized AIMD simulations in the canonical ensemble (NVT) with the Nosé-Hoover thermostat were performed to assess its dynamical stability. In the AIMD simulations, a 2 × 2 × 1 supercell was used, and each was performed for a total duration of 12.0 ps with each time-step of 1.5 fs at 300, 500, 700, and 800 K.

Minimum energy path calculations
In the simulation of the catalyzed chemical reaction, 2 × 2 × 1 supercells were used to model the 2D catalyst for the adsorption of the various molecules and the determination of the minimum energy pathway and the TSs using the climbing image nudged elastic band (CI-NEB) method 42 implemented in VASP. Here, the adsorption energy E ads of species i is defined by E ads i ð Þ ¼ E i@Catalyst À E Catalyst À E i where E i@Catalyst and E Catalyst are the total electronic energy of the 2D catalyst with and without the adsorbate species i, and E i is the electronic energy of the species i in the gaseous state.

DATA AVAILABILITY
The data presented in current study are available from the corresponding author on reasonable request.