Two-dimensional honeycomb-kagome V2O3: a robust room-temperature magnetic Chern insulator interfaced with graphene

The possibility of dissipationless chiral edge states without the need of an external magnetic field in the quantum anomalous Hall effect (QAHE) offers a great potential in electronic/spintronic applications. The biggest hurdle for the realization of a room-temperature magnetic Chern insulator is to find a structurally stable material with a sufficiently large energy gap and Curie temperature that can be easily implemented in electronic devices. This work based on first-principle methods shows that a single atomic layer of V2O3 with honeycomb-kagome (HK) lattice is structurally stable with a spin-polarized Dirac cone which gives rise to a room-temperature QAHE by the existence of an atomic on-site spin-orbit coupling (SOC). Moreover, by a strain and substrate study, it was found that the quantum anomalous Hall system is robust against small deformations and can be supported by a graphene substrate.


I. INTRODUCTION
To advance electronics at the nanoscale, the exploration and exploitation of quantum degrees of freedom in materials becomes indispensable. Two-dimensional (2D) compounds form excellent systems wherein these quantum degrees of freedom can be exploited, towards electronic and spintronic applications. A particular group of interest are 2D Dirac materials, which have a Dirac cone with linear dispersion in their band structure [1]. These materials have been predicted to exhibit a large variety of exotic properties: massless fermions [2], ultrahigh carrier mobility [3], (fractional) quantum Hall effects (QHE) [4,5], etc. Furthermore, there have been great advances in the preparation and growth of 2D materials in recent years [6][7][8], which initiated the study of the largely unexplored group of 2D strongly-correlated Dirac systems (SCDS). As the Dirac cone has shown to be instrumental for the realization of many non-trivial topological phases [9,10], these SCDS uncover a new rich playground where relativistic dispersion, electron correlations, and topological ordering meet.
In this work, we perform an ab initio study on the 2D single atomic layer of V 2 O 3 with HK lattice structure. Firstly, we confirm the predicted [11] ground state properties of HK V 2 O 3 and establish that it is an excellent SCDS candidate with a room-temperature QAHE. The second part of this work mainly focuses on the experimental feasibility by studying the structure under compressive and tensile biaxial strain as well as the graphene supported V 2 O 3 monolayer system. * simon.mellaerts@kuleuven.be

II. COMPUTATIONAL METHODS
All spin-polarized calculations were carried out within density-functional theory (DFT), as implemented in the Vienna ab initio simulation package (VASP) [12]. The generalized gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof (PBE) [13], projectedaugmented-wave (PAW) potential [14], and the planewave basis with energy cutoff of 550 eV were used. For the structural relaxation, a force convergence criterion of 0.005 eV/Å was used with the Brillouin zone (BZ) sampled by a 12 × 12 × 1 Γ-centered k-point mesh, with a vacuum space of 17 Å adopted along the normal of the atomic plane. To account for the localized nature of the 3d electrons of the V cation a Hubbard correction U is employed within the rotationally invariant approach proposed by Dudarev et al. [15], where U eff = U − J is the only meaningful parameter. The self-consistent linear response calculation introduced by Cococcioni et al. [16] was adopted to determine U . In this way, a Hubbard correction of U = 3.28 eV is found (see Supplementary Material), which is close to the U Bulk = 3 eV value found in bulk V 2 O 3 system [17], which is then applied throughout the paper.
The magnetic and electronic self-consistent calculations were performed with a total energy convergence criterion of 10 −6 eV with the BZ sampled by a denser 24 × 24 × 1 Γ-centered k-point mesh. To test the sensitivity of the results to the choice of the functional, the local density approximation (LDA) [18] and the screened exchange hybrid density functional by Heyd-Scuseria-Ernzerhof (HSE06) [19] are employed. Additionally, for the substrate calculations, the van der Waals (vdW) interactions were taken into account by the use of Grimme's DFT-D3 method [20]. The phonon dispersion was calculated self-consistently on the basis of the density functional perturbation the- ory (DFPT) and with the use of the PHONOPY package [21]. The ab initio molecular dynamic (AIMD) simulations were carried out on a 4 × 4 × 1 supercell at 300 K in the canonical ensemble using the Nose-Hoover thermostat approach [22,23] with 3000 times steps of step size 2 fs. The Curie temperature was estimated by Monte Carlo (MC) simulations, as implemented in the VAM-PIRE package [24], and by the use of a fully quantummechanical method where the hamiltonian is solved by Green's functions [25]. For the former calculations a rectangular 50 × 50 √ 3 supercell was used, where the spins were thermalized for 10 4 equilibrium steps, followed by 2 · 10 4 averaging steps to calculate the thermal averaged magnetization for each temperature. The atomic structures were visualized by the VESTA program [26].

A. Structural stability and mechanical properties
The structurally optimized unit cell of HK V 2 O 3 is shown in Fig. 1a with a planar honeycomb-kagome lattice with lattice constant 6.193 Å. A reduction of 0.15 Å of the V-O bond -from 1.94 Å to 1.79 Å -is found com-pared to the shortest V-O bond in bulk V 2 O 3 [27]. A similar reduction is observed in other metal oxides [28][29][30]. This can be attributed to the reduction in coordination number of the V 3+ cation in this 2D form, where V 3+ has electron configuration [Ar]3d 2 with the two outer 4s electrons of vanadium ionically bonded to the oxygen.
To study the dynamical stability of the monolayer system, the phonon spectrum was evaluated by DFPT. As shown in Fig. 1b, the spectrum shows no imaginary frequency phonon modes, confirming the dynamical stability. Furthermore, the AIMD simulation at a fixed temperature of 300 K also confirmed the thermal stability, since total energy fluctuations of only 0.3 % are observed while planarity is preserved, as shown in Fig. 1c-d. In addition to the stability verification, the elastic properties were determined. By symmetry only two independent stiffness constants C 11 and C 12 exist, where we neglect any out-of-plane bending contribution to the elastic energy density, i.e. the vertical weight of the stifness constants vanishes. The stiffness constants were determined by evaluating the total energy of the system under uniform uni-and bi-axial deformations, over a range of −10 % to +10 %. For all deformations, the atomic positions were optimized and the corresponding energies were determined. In this way, the stiffness constants, in-plane Young's modulus Y s = (C 2 11 − C 2 12 )/C 11 and Poisson's ratio ν = C 12 /C 11 were calculated. The resulting values are summarized in Table I, and compared to other 2D materials.
The mechanical stability could also be confirmed by verifying the positivity of the strain energy density function U ( ) = 1 2 C ijkl ij kl for all | | ≤ 10 %. For the 2D hexagonal systems this translates into satisfying the conditions C 11 > 0 and C 11 > |C 12 |. Hence, it can be concluded that the single atomic layer of V 2 O 3 with HK lattice structure is stable on all three levels: dynamical, thermal and mechanical.

B. Magnetic ground state
To determine the magnetic ground state of the system various magnetic configurations have been studied. The total energy of the ion-electron system was determined for ferromagnetic (FM), paramagnetic (PM), and four types of antiferromagnetic (AFM) states, depicted in Fig. 2a. By comparison of the total energies in Table II, it can be concluded that the system is in a FM ground state, with the AFM stripy (AFM-ST) state as the second lowest energy state. As a first approximation, the following spin hamiltonian was considered, where i enumerates the magnetic ions. The first term is the nearest-neighbor exchange interaction with exchange parameter J 1 , while the second term corresponds to the single-ion anisotropy with the anisotropy parameter A. The interaction is ferromagnetic when J 1 > 0 and positive A will favor an out-of-plane spin component S z . The exchange parameter is evaluated by J 1 = [E AF M − E F M ]/2zS 2 with the total spin moment S = 1, and with z = 3 neighboring spins. This gives J 1 = 58 meV. On the other hand, the magnetocrystalline anisotropy en- , with square brackets indicating the orientation of the spins, equals 2 meV/f.u. (i.e. A = 1 meV). This strong MAE can be understood from degenerate perturbation theory. The SOC-induced interaction involves the coupling between states with identical spin, with the most important interaction between the unoccupied and occupied d states. As can be seen in Fig. 3b, these states at the conduction band minimum (CBM) at Γ and valence band maximum (VBM) at K, respectively consist of the (d xy , d x 2 −y 2 )-and (d xz , d yz )-orbitals. This SOCinduced interaction involves a change in the orbital angular momentum quantum number This interaction between these out-of-plane spins leads to a maximization of the SOC energy stabilization, as described in [31]. To obtain an estimate of the Curie temperature (T C ), the Curie-Bloch equation in the classical limit, with T the temperature and β a critical exponent, is fitted to the normalized magnetization obtained by the MC simulation (see Fig. 2b). The resulting T C is equal to 447 K, and can also be compared within the mean field approximation by using the formula with Boltzmann constant k B , which gives 449 K. This is in relative close agreement with the MC result, which can be attributed to the relatively high MAE.
These MC simulations treat the spins as classical vector quantities. Alternatively, the semi-classical Holstein-Primakoff approximation is often employed. However, both approximations remain only justified for high spin values and very low temperatures. Therefore, a fully quantum mechanical method is employed, where the quantum Heisenberg hamiltonian is solved by Green's functions [25]. Since the Néel AFM phase is not the second lowest energy state, it can be expected that second-and third-neighbor exchange interactions (J 2 , J 3 ) are playing a dominant role in the magnetic ordering. Therefore, the DFT total energies were mapped onto the Ising hamiltonian to obtain J 1 , J 2 and J 3 by the following four equations: (4) and The J 1 , J 2 and J 3 values were found to be 41.9 meV, −3.8 meV and 17.4 meV, respectively. Taking into account these second-and third-neighbor exchange interactions in the spin hamiltonian (1), the hamiltonian is solved by the Green's functions, giving a Curie temperature equal to 488.46 K.

C. Electronic structure
The electronic band structure and density of states (DOS) within GGA+U are depicted in Fig. 3. The electronic band structure shows two key features: a Dirac cone at the K point and the existence of two flat-bands at ±1 eV forming a degeneracy with the dispersive bands at the Γ point. From the DOS it becomes clear that there exists a strong spin-polarization with the two characteristic band features being derived from the spin-up V (d xz , d yz )-orbitals. This means that the Dirac cone is completely spin-polarized, making it a Dirac spin-gapless semiconductor (DSGS) [32]. By taking into account the intrinsic SOC, it is found that an energy gap is opened at the K point, resulting in an indirect energy gap of 0.45 eV. It turns out that the energy gap is strongly dependent on the Hubbard correction U (see Supplementary Material), which has led to the idea of the cooperative effect between electron correlations and SOC. The relatively large energy gap opened at the Dirac point can be understood from the fact that the orbital degeneracy of the occupied states allows an atomic on-site SOC, involving no hopping process [33]. This is in contrast to the second-order hopping SOC in the Kane-Mele model, or first-order hopping SOC in the Xenes [34]. Within this understanding, it is clear that an increase of U will force the electrons into a more localized state implying an increased effect of the atomic SOC [35,36].
In addition to the GGA+U , the LDA and hybrid exchange-correlation functional HSE06 were used to study the electronic properties of the system. The resulting energy gap values are summarized in Table III, and the calculated band structures can be found in the Supplementary Material. It is noted that there is a strong variation of the energy gap depending on the chosen functional. This can be expected as these functionals are approximating the localized nature of the 3d orbitals very differently, which will be crucial for the atomic SOC effect. From the DOS it can be inferred that the crystal field splits the d-orbitals into three orbital levels whose energy ordering can be derived from the orbital orientation. The d z 2 -orbital level has negligible overlap with the O 2p orbitals and thus has the lowest energy state, the (d xz , d yz )-orbitals oriented out-of-plane will form a π-bond with the oxygen O p z as bridging ligand, and the (d xy , d x 2 −y 2 )-orbitals will form an in-plane σ-bond with bridging orbitals O (p x , p y ). This latter bond involves a strong orbital overlap and therefore the involved orbitals will be energetically less favorable, forming the highest energy state. The resulting orbital ordering is shown in Fig. 3e. It should be emphasized that in the vicinity of the Fermi level, there is approximately no O 2p orbital weight, which confirms the localized nature of the two V 3d electrons. Based on Griffith's crystal field theory [37], the spin state of the V 3+ cation can be confirmed. By the relative strengths of the exchange (∆E ex ) and crystal field splitting (∆E CF ) with derived values of 3.34 eV and 0.8 eV -it can be concluded that the V cation has a high (S = 1) spin state with magnetic moment 2µ B .
To gain further insight in the nature of the chemical bond, the atomic charges were determined by the Bader charge analysis code [38]. It is found that the  atomic charges on V and O are resp. Q V = 1.59e and Q O = −1.07e, confirming the ionic bonding character where the electrons of the V cation are attracted towards the O anion. Nevertheless, the atomic charges are slightly lower than in the bulk parent, which suggests that there is an increased electron delocalization and thus stronger bonding covalency compared to the bulk. This behavior can be linked to the reduced V-O bond length. On the other hand, the relative strong ionicity of the bond might explain the energetic stability of the planar configuration of these TMO monolayer systems, since the buckling of the V-O bond would result in an energetically unfavorable large dipole moment normal to the atomic plane. In this way, the reduced bond length and enhanced covalency can be explained as a means to suppress the possible large dipole moment. These trends and their explanations were already pointed out for other TMOs [28,29].
To study the topology of the band structure the Berry curvature Ω(k) of the system was calculated directly from the DFT calculated wave functions by the VASP-BERRY code, which is based on Fukui's method [39]. The resulting Berry curvature, shown in Fig. 3c, becomes non-zero at the Dirac points K and K . By integrating the Berry curvature Ω n (k) of each n-th energy band over the whole BZ and summing over all occupied bands, the Chern number is found to be C = 1. Therefore, it can be concluded that the system is a room-temperature Chern insulator.

D. Biaxial strain
TMOs are known to exhibit a rich phase diagram as a function of strain, in particular, for bulk V 2 O 3 the roomtemperature metal-insulator transition (MIT) between the paramagnetic insulating and metallic state can be realized by the application of epitaxial strain [40]. Therefore, a biaxial strain study on the 2D atomic layer of V 2 O 3 was performed. The atomic positions were optimized for all biaxial strained configurations from −10 % to +10 %. It was found that under compressive strain, the lattice structure undergoes an in-plane buckling of the V-O bond (HK in ), as can be seen Fig. 4a-b. This structural distortion is found by ionic relaxation with high force convergence criterion of 0.005 eV/Å. On the other hand, by lowering the force convergence criterion to 0.05 eV/Å, two metastable structures were found in the compressive region; one preserving the HK lattice structure (HK), while the other undergoes an out-of-plane buckling of the V-O bond (HK out ). Although these lattice configurations are un/meta-stable, they could become stabilized under certain conditions and depending on underlying substrates. Moreover, a combination of inand out-of-plane buckling can also be expected and was observed in an ab initio calculation of TMOs on metal substrates [41].
Identifying these different structural configurations under compressive strain in the freestanding V 2 O 3 monolayer allows a further study on the structural stability and provides a better understanding of the involved chemical bonds governing the stability of this planar single atomic layer. From the structural in-and out-of-plane distortions, it is noted that the deformation of the V-O bond develops as to prevent a further reduction of this bond length (see Fig. 4c). This behavior can be linked to the ionic character of the bond. On the other hand, from the energetic point of view, the system prefers to undergo an in-plane deformation rather than the out-of-plane defor- mation. This can be attributed to the energy cost of the formation of a dipole moment under a buckling of the V-O bond.
In addition to this structural study, the electronic and magnetic properties as a function of the biaxial strain were studied. Since the stabilization of the metastable states (HK and HK out ) will depend on external conditions, only the most stable HK in configuration is considered for the compressive region. The energy gap and the magnetic exchange interaction are shown in Fig. 4. It is clear that within the strain range of −5% to +5%, the electronic and magnetic properties are sufficiently preserved to maintain the room-temperature magnetic Chern insulating phase. It was also confirmed that the Chern invariant C remains equal to 1, even in the presence of the in-plane distortion of HK in . This can be expected as the V cations still form a honeycomb lattice, crucial for the existence of the Dirac cone [42]. From these results, it can be concluded that the room-temperature QAHE in 2D V 2 O 3 is robust against small structural deformations.

E. Graphene substrate
Motivated by the successful deposition of Y 2 O 3 on graphene [43], the feasibility of using a graphene substrate for the growth of the V 2 O 3 monolayer was studied. A (2 × 2) V 2 O 3 supercell on a (5 × 5) graphene supercell was structurally optimized and that configuration is shown in Fig. 5a (see Supplemental Material). Similar to Al 2 O 3 and Y 2 O 3 monolayers on graphene [28,29], the V cations are located near the top sites of the C atoms to maximize the potential bonding. The planar HK lattice structure is conserved with a minimal distance between graphene and V 2 O 3 of 3.43Å, while the V-O bond length remains approximately preserved.
To obtain a first estimate of the interfacial interaction strength, the adsorbtion energy E ad is determined, where E G , E V2O3 and E G+V2O3 are the total energies of the isolated graphene, isolated V 2 O 3 , and the combined hybrid structure, respectively. The adsorbtion energy equals 77 meV per C atom indicating a weak interaction between both monolayers.
To preserve the electronic properties it is important to prevent charge transfer between both monolayers, therefore the atomic charges of the systems were determined for the hybrid system. The average atomic charges of both V and O remain unchanged, while the average excess charge of the C atom is only 0.004e. This indicates negligible electron transfer between both monolayers, which can also be confirmed by the calculation of the charge density difference ∆ρ = ρ G+V2O3 − ρ G − ρ V2O3 , shown in Fig. 5b. The charge density difference shows a negligibly small inhomogeneous charge distribution between both monolayers, which is mainly located between V and C atoms. This indicates a small orbital hybridization of the π-bonded orbitals of V 2 O 3 and graphene. This is similar to earlier studies of HK TMOs on graphene where it was shown that vdW interactions and orbital hybridization play an important role for the electronic properties at the interface [28,29].
The electronic properties of the V 2 O 3 monolayer on the graphene substrate are shown in Fig. 5c-d. From the band structure and DOS, it is noted that the conduction band minimum of graphene is about 0.2 eV below the Fermi level, indicating that graphene is slightly n-doped. The DOS contribution of the V 2 O 3 remains approximately preserved, however from the calculated band structure, it is clear that there is some hybridization of the π-bonded orbitals. Nonetheless, the Dirac cone of V 2 O 3 remains preserved with a reduction in the SOCinduced direct energy gap to 66 meV.

IV. CONCLUSIONS
In this work, we have confirmed that the single atomic layer of V 2 O 3 with HK lattice structure is a structurally stable room-temperature magnetic Chern insulator. This system features the coexistence of topological order and strong electron correlations, and might therefore be an excellent SCDS candidate. It was shown that the system can undergo small structural deformations which preserve the honeycomb lattice of the V cations, such that the Dirac cone and the corresponding room-temperature QAHE remain unaffected. Furthermore, this 2D V 2 O 3 can be further stabilized by the support of a graphene substrate, where there is only a small interfacial interaction mainly attributed to orbital hybridization, but preserving the Dirac cone of both materials. These observations together with earlier studies on TMO with HK lattice structure [28][29][30][43][44][45] should encourage the experimental investigation of this group of 2D TMOs.

DATA AVAILABILITY
The data that supports the findings of this study are provided here and the Supplementary Material.

Stiffness constants
From continuum mechanics, a general equation for the elastic energy for a deformation of a 2D material can be obtained: where E 0 is the equilibrium energy, A 0 the equilibrium unit cell area, σ ij the rank 2 stress tensor, ij the rank 2 strain tensor, and C ijkl the rank 4 stiffness tensor. To simplify this equation with high-rank tensors, the Voigt notation [? ] will be adopted. Furthermore, by the symmetries of the hexagonal lattice, the stiffness matrix reduces to a 3 × 3 matrix of the form [? ] with three independent stiffness constants C 11 ,C 12 and C 33 . Since only in-plane deformations are considered, C 33 is set to zero. To determine the stiffness constants, we consider a uni-and biaxial deformation respectively given by the strain tensors with δ a dimensional quantity, representing the stength of the deformation, defined by the relative change in the lattice constant along the axis of deformation, where a 0 and a are the lattice constants along the deformation direction of respectively the equilibrium and deformed structure. Using Eq. 8, the energy change associated with both deformations is given by The calculated total energies of the deformed configurations are shown in Fig. 6 and the corresponding functions in Eq. 12 are fitted to these energies.

Estimate of the Hubbard correction
This linear response calculation is based on [16] and is divided into three steps. First, a self-consistent electronic calculation is performed without the presence of any perturbation α. Secondly, a non-charge-self-consistent calculation is performed with the presence of a non-zero pertubartion α. In this calculation, the charge density is not allowed to relax and thus a bare (non-interacting) response is obtained. Lastly, a charge-self-consistent claculation with nonzero α is employed where the charge density is allowed to relax and to screen the perturbation. This provides the interacting response of the system. Both bare and interacting response can be respectively quantified by the response density functions [16]: where I, J are the V-sites in the supercell, and the KS superscript stands for the non-interacting (Kohn-Sham) response. In order to obtain an accurate estimate of U , the calculation of the response density functions should be repeated for increasing supercell dimensions, until there is convergence. However, a 2 × 2 × 1 supercell containing 24 atoms should suffice for this calculation, with larger supercells giving minor improvements. Subsequently, the effective interaction parameter U can be determined by where IJ represents the evaluation of the response density for a perturbation a site J and its response measured at site I. These calculations were performed for a perturbation range of 0.4 eV with stepsize 0.05 eV. By linear fitting, illustrated in Fig. 7a, the response density functions were determined from the slope and using Eq. 15 a Hubbard correction of U = 3.28 eV was found.
Electronic structure in function of exchange-correlation functional

Effect of Hubbard correction
The electronic band structure was calculated for a Hubbard correction ranging from 2.5 eV to 4.5 eV. Without the inclusion of the spin-orbit coupling (SOC), it is observed that there is less overlap of the different orbital levels for increasing U , nonetheless, the overall effect is relatively small. By the inclusion of SOC, it is clear that the energy gap opened at the Dirac point strongly depends on the Hubbard correction U , as shown in Fig. 7b.

Exchange-correlation functional dependence
To test the sensitvity of the results of the electronic band structure to the choice of the functional, the three most common exchange-correlation functionals were employed. The generalized gradient approximation (GGA) in the form of Perdew-Burke Ernzerhof (PBE) [13], the local density approximation (LDA) [18] and the screened exchange hybrid density functional by Heyd-Scuseria-Ernzerhof (HSE06) [19]. The calculated band structures for LDA+U and HSE06 are shown in Fig. 8 and compared with the band structure obtained by GGA+U .

Topology of the band structure
In addition to the topological energy gap at the Fermi level with a Chern number equal to one, the Chern number was also calculated for the gap opening at the Γ point around −0.9 eV below the Fermi level. The calculated Berry cruvature is shown in Fig. 9 and was integrated by the VASPBERRY code, which is based on Fukui's method [39].

Graphene substrate
For the ionic relaxation of the graphene supported V 2 O 3 monolayer system, several configurations of the V 2 O 3 monolayer with respect to the graphene layer have been considered, as shown in Fig. 10. For all configurations, the supercell was constructed with completely relaxed structures of graphene and V 2 O 3 , with freestanding lattice constants of 2.477 Å and 6.193 Å, respectively. By comparing the total energies of the different configurations, it was confirmed that the other considered configurations are very close in energy (< 10 −3 eV), and as the superstructure of graphene on V 2 O 3 is incommensurate it is expected that for any configuration the V cations will be located approximately on top of the C atoms. Furthermore, it was verified that the density of states remains approximately preserved in these other configurations. At last, the initial distance between both monolayers was varied to confirm absolute stability. The most stable configuration was found to have unchanged structural parameters in both monolayers, i.e. the vdW epitaxy does not induce any significant strain within the system.  Figure 10. Two different superstructures with the (2 × 2) V 2 O 3 monolayer with respect to the (5 × 5) graphene layer. The atom colors gray, cyan, and red corresponding to carbon, vanadium and oxygen respectively.