Realization of Kane-Mele Model in $\boldsymbol{X}\bf{N_{4}}$-Embedded Graphene ($\boldsymbol{X}$=Pt, Ir, Rh, Os)

Introduction—. Two-dimensional (2D) $\mathbb{Z}_{2}$ topological insulators (TIs) are characterized by a pair of spin-helical counterpropagating edge modes, which are robust against non-magnetic disorder due to the protection of time-reversal symmetry ref1 ; ref2 . This striking electronic transport phenomenon allows dissipationless electron transport without elastic backscattering, and thus offers a promising platform to realize novel low-power electronic and spintronic devices. 2D $\mathbb{Z}_{2}$ TI was first theoretically proposed by Kane and Mele in monolayer graphene ref4 , and then realized in HgTe/CdTe as well as InAs/GaSb quantum wells ref22 ; ref23 ; ref5 . Recent years have also witnessed abundant theoretically proposed candidates of 2D $\mathbb{Z}_{2}$ TIs ref55 ; ref6 ; ref24 ; ref25 ; ref7 ; ref26 ; ref27 ; ref28 ; ref29 ; ref30 ; ref31 ; ref32 ; ref33 ; ref34 ; ref35 ; ref36 ; ref37 ; ref38 ; ref39 , among which only Bi-based thin films and Jacutingaite have been experimentally confirmed ref44 ; ref45 ; ref46 ; ref47 . Nevertheless, observation of the 2D $\mathbb{Z}_{2}$ TI in graphene-based systems has not been experimentally realized yet. As the first and most appealing candidate for 2D $\mathbb{Z}_{2}$ TI, graphene has been attracting numerous focus for over fifteen years. However, although graphene exhibits various extraordinary properties ref8 , its negligible spin-orbit coupling (SOC) precludes a sizeable band gap to realize the Kane-Mele type $\mathbb{Z}_{2}$ TI.

Several rewarding decorating approaches have been proposed to enhance the SOC in graphene, e.g., atomic adsorption ref9 ; ref10 ; ref11 , atomic intercalation ref40 ; ref41 , proximity effect ref49 , or atomic embedding. All of them are challenging so far for experimental realization, except for embedding metal atoms, which are however mobile under irradiation to form a buckled lattice that breaks the flatness of graphene and thus loses its high mobility ref59 ; ref60 . Recently, an improved strategy, i.e., embedding metal atoms together with nitrogen atoms into graphene, becomes practical, as some experiments have demonstrated ref18 ; ref19 ; ref21 that the nitrogen atoms are able to anchor the metal atoms in the coplanar lattice of graphene, facilitating to form a stable 2D metal compound, named after $X$N₄-embedded graphene with $X$ being the metal atom. Compared with other types of modified graphene mentioned above, $X$N₄-embedded graphene has intriguing advantages that are demonstrated by both experiments and first-principles calculations, including that (i) it maintains a flat single-atom-thick structure as pristine graphene ref18 ; ref19 ; ref21 ; (ii) it covers twenty four kinds of metal elements with different atomic coverages ref20 ; (iii) its structural stability of CoN₄-embedded graphene has been ascertained via phonon spectra calculations and molecular dynamic simulations ref12 . As the heavy transition metal atoms are enabled to be embedded into graphene, it may exhibit considerable SOC to realize the expectant topological states, e.g. $\mathbb{Z}_{2}$ TIs.

In this Letter, we construct a representative structural model of $X$N₄-embedded graphene and investigate its topological properties by using the density functional theory and tight-binding model calculations. We demonstrate that four kinds of 2D materials for $X$N₄-embedded graphene, as three of which (i.e., $X=$Pt, Ir, Rh) have been experimentally fabricated, are 2D $\mathbb{Z}_{2}$ TIs. We show that in the absence of SOC, they exhibit spin-degenerate gapless Dirac dispersions around two Dirac points that are away from K/K’ points due to the atomic structural distortion. When SOC is further switched on, a sizeable band gap of 47/101 meV in (Pt/Os)N₄C₁₀ opens at the Fermi level, while a 46/25 meV gap in (Ir/Rh)N₄C₁₀ opens above the Fermi level. The $\mathbb{Z}_{2}$ index, edge states, spin Berry curvature and spin-Hall conductivity are calculated to confirm the topologically nontrivial properties. At the low-energy continuum limit, a modified Kane-Mele model is constructed to naturally captures the physical mechanism of the gapped Dirac cone in $X$N₄-embedded graphene. Finally, we illustrate that either the random distribution or the variable coverage of the embedding $X$N₄ units plays negligible role in the formation of Dirac cone structures, indicating that the observation of $\mathbb{Z}_{2}$ phase in $X$N₄-embedded graphene is experimentally feasible.

Methods—. Our first-principles calculations based on density functional theory are performed using the Vienna Ab initio Simulation Package (VASP) ref13 ; ref14 with the project augmented-wave method ref15 and exchange-correlation functional in the Perdew-Berke-Ernzerhofs form within the generalized-gradient approximation (GGA) ref16 . A vacuum buffer space over $16$ Åis applied to minimize image interactions between adjacent slabs. All atoms are allowed to be fully relaxed until the force on each atom is less than $10^{-6}$ eV/Å. The plane-wave energy cutoff is set to be 600 eV, with the first Brillouin-zone integration being carried out using a $15\times 15\times 1$ $\Gamma$-centered $k$-point grid. Phonon band dispersions are evaluated by using density functional perturbation theory based on the PHONOPY program ref42 . Maximally localized Wannier functions are constructed by using the WANNIER90 package ref17 to study their topological properties.

Structural and Electronic Properties—. At each $X$N₄ unit, the central $X$ atom is connected with four N atoms. There are usually different coverage configurations of $X$N₄ units in graphene. Here, we choose a representative structure named after $X$N₄C₁₀, where $X$ is a transition metal atom. Hereinbelow, we take PtN₄C₁₀ (see Fig. 1(a)) as a concrete example. In our calculations, all atoms in $X$N₄C₁₀ are coplanar to form a stable structure, indicating that it is a well-defined 2D system. Furthermore, the embedded PtN₄ units deform the pristine graphene and lead to the variation of the covalent bonds between carbon atoms, i.e., there are three inequivalent carbon atoms, ${\rm C_{1}}$, ${\rm C_{2}}$ and ${\rm C_{3}}$ in PtN₄C₁₀, with four different bonds as listed in Table 1. As a consequence, $X$N₄C₁₀ possesses the rhombic lattice with a plane symmetry group $Cmm$, affiliated to the point group $D_{2}$, which has the two-fold rotation symmetry along $x$/$y$/$z$ axis and is identical to the mirror symmetry of $z$-$x$/$y$-$z$ plane (denoted by $\mathcal{M}_{y/x}$) and the inversion symmetry, respectively.

Figure 1(c) displays the phonon spectrum. One can see that there is no trace of imaginary frequency in the Brillouin zone, except a small portion of negative value near $\Gamma$ point. The limited negative imaginary frequency may arise from the difficulty of achieving numerical convergence for the flexural phonon branch, which usually appears in first-principles calculations for 2D materials ref57 ; ref58 . Therefore, PtN₄C₁₀ can be considered to be dynamically stable.

The electronic structure of PtN₄C₁₀ shows a bulk band gap of tens of meV. In the absence of SOC, the valence and conduction bands touch at the Fermi level to form a distorted gapless Dirac cone, as displayed in Fig. 2(a). The density of states is vanishing at the Fermi level, as displayed in Fig. 2(b), which further confirms the existence of Dirac cone structures. In graphene, the two inequivalent Dirac cones at K/K’ points are protected by the three fold rotational symmetry. A distortion of the honeycomb lattice can cause a shift of the Dirac cone in the (distorted) first Brillouin zone with lower spatial symmetry, which is commonly reported in deformed graphene, e.g., phagraphene ref53 . For the electronic structure of PtN₄C₁₀, two Dirac points lie in the high symmetry lines between $\Gamma$ and Y, and $\Gamma$ and -Y, respectively. Therefore except along $k_{x}$ direction, the bands of PtN₄C₁₀ near the Dirac cone have anisotropic linear dispersions, meaning anisotropic Fermi velocities. Furthermore, the PtN₄ unit with $D_{2}$ symmetry generates a deformed square planar crystal field in the structural plane, which splits Pt $d$ orbitals into four distinct energy levels denoted by $E_{1g}(d_{zx},d_{yz}),A_{1g}(d_{z^{2}}),B_{1g}(d_{x^{2}-y^{2}}),B_{2g}(d_{xy})$. The C $p_{z}$ and N $p_{z}$ orbitals form $\pi$ bonds together with Pt $d_{zx}$ and $d_{yz}$ orbitals, because they have the same odd parity, as their projected density of states are plotted in Fig. 2(b). The other $d$ orbitals with the even parity, on the other hand, do not couple with the above orbitals and thus form three discrete flat bands, in which two ($d_{x^{2}-y^{2}}$ and $d_{z^{2}}$) lie below the Fermi level and one ($d_{xy}$) lies above the Fermi level. After introducing SOC, a bulk band gap around 47 meV opens at the Dirac point as displayed in Fig. 2(c), indicating an insulating state.

Topological Properties—. The gapped PtN₄C₁₀ is a 2D $\mathbb{Z}_{2}$ TI. The $\mathbb{Z}_{2}$ topological number can be evaluated from the Wannier charge center (WCC) ref50 ; ref51 . The evolution curves of WCC are obtained by using the approach implemented in WannierTools packages ref43 [see Fig. 3(b)], where the curves cross the reference line odd times, indicating $\mathbb{Z}_{2}=1$. We further calculate the edge densities of states by means of Green’s function technique in a semi-finite ribbon with armchair edges, as displayed in Fig. 3(a). The helical edge states connect the bulk conduction and valence bands, which characterizes the 2D $\mathbb{Z}_{2}$ TIs.

Underlying Physics from Tight-binding Model—. To further understand the physical mechanism of the $\mathbb{Z}_{2}$ phase of PtN₄C₁₀, we construct a tight-binding model by using Slater-Koster approximation on the orthogonal basis of {$|p^{\rm C}_{z}\rangle$, $|p^{\rm N}_{z}\rangle$, $|d_{zx}\rangle,|d_{yz}\rangle$} $\otimes$ {$|\uparrow\rangle$  $|\downarrow\rangle$}. The tight-binding model Hamiltonian can be written as:

where $c^{\dagger}_{i\alpha}=(c^{\dagger}_{i\alpha\uparrow},c^{\dagger}_{i\alpha\downarrow})$ denotes the creation operator of an electron at the $i$-th atomic site with $\uparrow$/$\downarrow$ and $\alpha$/$\beta$ representing spin up/down and different orbitals, respectively. The first term stands for the on-site energy. The second term is the nearest-neighbor hopping with an amplitude of $t_{ij,\alpha\beta}$, which are functions of the Slater-Koster parameters listed in Table 2. The third term represents the atomic SOC of Pt $d$ orbitals with a strength of $t_{\rm SO}$. $\bm{s}=(s_{x},s_{y},s_{z})$ and $\bm{l}=(l_{x},l_{y},l_{z})$ are the Pauli matrices and orbital-angular-momentum operators, respectively. Due to the lattice deformation, the hopping amplitude between carbon atoms changes with the interatomic distance and can be determined by $t_{ij}=V^{CC}_{pp\pi}e^{-\beta(d_{ij}/d_{0}-1)}$, where $V^{CC}_{pp\pi}$ is the Slater-Koster parameter and $d_{ij}$ is the distance between two carbon atoms. Here, we set $\beta=2.20$ and $d_{0}=1.42$ Å. Using these fitting parameters that are extracted from first-principles calculations, one can obtain the band structure that agrees quantitatively with the first-principles results near the Dirac point, as shown in Figure 2. When $t_{\rm SO}=0$, a gapless Dirac cone band is formed; when $t_{\rm SO}\neq 0$, a band gap opens at the Dirac point.

To further verify the topologically nontrivial band gap, we first calculate the spin Berry curvature $\Omega^{s}(\textbf{k})$, which is well-defined for flattened single-atom-thick systems with the mirror symmetry $\mathcal{M}_{z}$ about the $x$-$y$ plane:

where $\psi_{n\textbf{k}}$ and $\varepsilon_{n}$ are respectively the Bloch eigenstate and eigenvalue of the $n$-th band at k, $f_{n}$ is the Fermi-Dirac distribution function, and $j_{x}$ is the spin current operator defined as $\frac{1}{2}(s_{z}v_{x}+v_{x}s_{z})$ with $v_{x(y)}$ being the velocity operator and $s_{z}$ the spin operator. Figure 3(c) plots the spin Berry curvature of all the occupied bands. The sharp peaks appear at the positions of Dirac gap in the Brillouin zone. The spin-Hall conductivity directly related to the spin Berry curvature can be thus evaluated via the Kubo formula ref48 :

Figure 3(d) plots the spin-Hall conductivity $\sigma_{xy}^{s}$ as a function of Fermi energy ($E_{\rm F}$) for PtN₄C₁₀. One can observe that $\sigma_{xy}^{s}$ keeps quantized whenever $E_{\rm F}$ lies inside the band gap.

Low-Energy Effective Model—. In the low-energy limit, the electronic structure can be captured by the Kane-Mele model constraint by the symmetry of the system, which is even under mirror, parity and time reversal ref3 ; ref4 ; ref55 ; ref56 . Nevertheless, due to the decrease of the spatial symmetry, the fermi velocity is no longer isotropic. We thus modify the Kane-Mele model as $H_{\rm eff}=H_{0}+H_{\rm SO}$ around the Dirac points with $H_{0}=Ak_{x}\sigma_{x}\tau_{z}+Bk_{y}\sigma_{y}+Ck_{y}\tau_{z}$ and $H_{\rm SO}=\Delta_{\rm SO}\sigma_{z}\tau_{z}s_{z}$. Here, $\vec{\sigma}$ denotes the Pauli matrices, $\tau_{z}=\pm 1$ describes the states at two Dirac points in the first Brillouin zone [Fig. 3(c)], and $s_{z}$ represents electron’s spin. A, B and C are parameters in the unit of eV. The pair of gapless Dirac cones are protected by $\mathcal{M}_{x}$ and $\mathcal{M}_{y}$, which thus allow an additional term that is proportional to $k_{y}\tau_{z}$. $H_{\rm SO}$ is the intrinsic SOC in the Kane-Mele model, which generates a topologically nontrivial indirect gap $2\sqrt{1-C^{2}/B^{2}}\Delta_{\rm SO}$, as long as $C^{2}<B^{2}$. By fitting parameters, the bands exhibit good agreement with first-principles calculations (see Supplemental Materials ref52 ). Since our calculations of PtN₄-embedded graphene with different coverages suggest that the Dirac cone retains with the decrease of coverages ref52 , the effective model offers a very precise description of the Dirac cone in the low-energy range, and captures its odd $\mathbb{Z}_{2}$ phase. Thus, we show that PtN₄-embedded graphene belongs to the modified Kane-Mele type 2D $\mathbb{Z}_{2}$ TI.

Other Materials Candidates—. In Table 3, we list the parameters of dozen systems of $X$N₄C₁₀, with $X$ being the transition metal atoms. One can find that $X$N₄C₁₀ ($X$= Ir, Rh, and Os) with the same rhombic crystal structure also have sizeable nontrivial bulk gaps. The corresponding structural stability, electronic properties and topological properties are further provided in Supplemental Materials ref52 . It is worth mentioning that the structure of OsN₄C₁₀ is mechanically stable via our phonon spectrum calculation, and has a large topologically nontrivial gap around 101 meV at the Fermi level, though the experimental synthesis of OsN₄-embedded graphene has not been reported yet ref20 . Furthermore, other $X$N₄C₁₀ ($X$= Ni, Zn, Pd, Ag, Cd, and Au) also host the locally gapped Dirac cone structures, which are however overwhelmed by the first valance bands that cross the Fermi level. Such $\mathbb{Z}_{2}$ metallic states, however, could be tuned into the normal $\mathbb{Z}_{2}$ insulating states via doping non-magnetic disorder.

Random embedding distribution—. In realistic synthesized compounds in experiments ref18 ; ref19 ; ref21 ; ref20 , graphene is no longer subjected to periodically distributed embedding $X$N₄ units. It is thus reasonalbe to explore the effect of random distribution on their topological properties. From above analysis, the $X$N₄ units here play an essential role as enhancing intrinsic SOC, which is very similar to that of transition metal atoms adsorbed on graphene ref11 ; ref54 . In that case, random adatom distribution helps weaken inter-valley scattering while it has a negligible effect on the SOC. Besides, in our calculations ref52 , the Dirac cone structures are still retained at low coverages of the $X$N₄ units, hence the odd $\mathbb{Z}_{2}$ phase is expected to be stable in the presence of randomly distributed embedding $X$N₄ units.

Summary—. We systematically investigate the topological properties of $X{\rm N_{4}}$-embedded graphene. We demonstrate that ${X\rm N_{4}C_{10}}$ ($X$= Pt, Ir, Rh, Os) are 2D $\mathbb{Z}_{2}$ TIs. The embedded transition metal maintains the mechanic stability of graphene and modifies the Dirac cone structures of graphene to foster the stronger SOC. Without SOC, there are shifted Dirac bands with a vanishing band gap near the Fermi level. After introducing SOC, a topologically nontrivial band gap opens, characterized by the odd $\mathbb{Z}_{2}$ and topological edge states. A tight-binding model is constructed to reveal that the atomic SOC can induce such a nontrivial band gap. The low-energy bands of the gapped Dirac cone are described by a modified Kane-Mele model with the lower symmetry due to the deformation of graphene, indicating $X{\rm N_{4}}$-embedded graphene is 2D Kane-Mele type TIs. The bulk gap retains with a lower coverage of metal atoms, which is experimentally feasible. Our study of this system opens up routes towards exploring more graphene-based 2D $\mathbb{Z}_{2}$ TIs.