Realizing Haldane Model in Fe-based Honeycomb Ferromagnetic Insulators

The topological Haldane model (THM) on a honeycomb lattice is a prototype of systems hosting topological phases of matter without external fields. It is the simplest model exhibiting the quantum Hall effect without Landau levels, which motivated theoretical and experimental explorations of topological insulators and superconductors. Despite its simplicity, its realization in condensed matter systems has been elusive due to a seemingly difficult condition of spinless fermions with sublattice-dependent magnetic flux terms. While there have been theoretical proposals including elaborate atomic-scale engineering, identifying candidate THM materials has not been successful, and the first experimental realization was recently made in ultracold atoms. Here we suggest that a series of Fe-based honeycomb ferromagnetic insulators, AFe2(PO4)2 (A=Ba,Cs,K,La) possess Chern bands described by the THM. How to detect the quantum anomalous Hall effect is also discussed.

Introduction.-In 1988, F.D.M. Haldane introduced an idea of the quantum Hall effect without Landau levels, and a simple tight-binding model of spinless fermions on a honeycomb lattice was suggested as an example 1 , which was dubbed the topological Haldane model (THM). It features a chiral edge spectrum with a Chern number without external magnetic field which is a prototype of the quantum anomalous Hall (QAH) insulator [2][3][4] . Although THM was "unlikely to be directly physically realizable", as Haldane stated in his paper 1 , yet his vision of intrinsic topological state of matter in condensed matter systems inspired later discoveries of timereversal symmetric topological insulators (TI) and promoted other topological phases [5][6][7] .
The THM is a spinless fermion model in a honeycomb lattice with nearest neigbor (n.n.) and complex next nearest neighbor (n.n.n.) hopping integrals 1 : where t 1 and t 2 are real and represent n.n. and n.n.n. hopping integral terms, respectively. Φ ij breaks time-reversal symmetry (TRS), and its sign differs for two sublattices (A, B), i.e., Φ for A and −Φ for B. Realization of the THM requires spinless fermions hopping on honeycomb lattice with spatially alternating flux yielding Aharonov-Bohm phase iΦ ij . Due to these difficult requirements of the THM, realization of QAH effect in materials was achieved only after the discovery of quantum spin Hall (QSH) insulator. Since each spin component of electrons in QSH insulators is regarded as a QAH state, one can obtain a QAH effect if one spin component QAH state is removed. The discovery of TI 8,9 , together with recent advancement of atomic-scale engineering techniques, then revived the interest for the QAH phase. There have been a surge of theoretical proposals in various system including magnetic-ion-doped HgTe quantum well 10 and TI surfaces 11 , engineered graphene 12,13 , transition metal oxides 14,15 and their heterostructures 16,17 . On the other hand, experimental observations of QAH effect was reported only in Cr-and V-doped (Bi,Sb) 2 Te 3 film 18,19 , confirming the idea that a TI with magnetic impurities removes one spin QAH state and reveals the QAH effect. Alternatively, the idea of breaking TRS by exerting circularly polarized ac-electric field and inducing QAH phase was suggested in light of the Floquet-Bloch theory 20,21 and has been realized recently 22,23 . Also, the experimental realization of the THM was recently made in ultracold atomic fermions in a periodically modulated honeycomb lattice 24 . However, it seems that realizing the THM in a simple two-dimensional honeycomb compounds in an equilibrium situation becomes at a glance an unrealistic task.
Here we show that the THM, original QAH model can be found in Fe-based honeycomb ferromagnetic insulators. With the help of strong Hund's coupling in Fe, electrons with one major spin-component (say down-spin) are fully polarized in occupied bands. Then electrons with other spin component (up-spin) form Haldane bands with finite Chern numbers, described by effective spinless fermions with complex n.n.n. hoppings. We find that a series of Fe-based honeycomb stoichiometric materials, AFe 2 (PO 4 ) 2 (AFPO, A=Ba,K,Cs,La), fall into a class of these materials described by the THM. Among them, compounds with A=K,Cs, and La exhibit a QAH effect, while BaFe 2 (PO 4 ) 2 does not show a QAH effect, because two Chern bands have opposite chirality.
BaFe 2 (PO 4 ) 2 (BFPO), a recently synthesized insulator, is the first example of two-dimensional Ising ferromagnetic oxides, where honeycomb layers are made of FeO 6 octahedra 25 Ferromagnetic transition occurs at T c ∼ 65 K, and Fe 2+ (d 6 ) high-spin moments on honeycomb lattices align along the layer-normal direction below T c . 26 . Interestingly, it also shows an intriguing re-entrant structural transition at T c from monoclinic P1 to rhombohedral R3 symmetries with most likely due to the coupling between ferromagnetic ordering and lattice structure via spin-orbit coupling (SOC). Signature of orbital angular momentum of ∼ 1µ B at Fe site was reported, implying significant role of atomic SOC in BFPO 27 . It was suggested that electronic correlation turns the system from a semimetal to a Mott insulator with Fe atomic orbital angular momentum 28 . We find that the atomic orbital momentum in BFPO signals possible Haldane bands via combined effects of Coulomb interactions and SOC under the ferromagnetic order. Two copies of Haldane Chern bands are identified; one set of Chern bands just above and the other set of Chern bands below the Fermi level, with opposite chirality Results.- Fig. 1 shows the evolution of Fe d-orbital in BFPO with respect to inclusions of relevant energy scales. In a 3d 6 configuration in the atomic limit, the dominant energy scale is the exchange splitting between different spin components introduced by the Hund's coupling. Density functional theory (DFT) calculation yields the energy difference of 5J H ∼ 3eV between the down and up spin states 28 , which is much larger than the cubic crystal field of ∼ 1 eV exerted by an oxygen octahedral cage surrounding Fe. Hence five electrons with down-spin occupy the bands well below the Fermi level, and only one electron with up spin is left in the t 2g minor spin states, which acts as spinless fermion. The system without SOC becomes half-metallic when no further on-site energy scales are included, as shown in Fig. 1. The presence of trigonal distortion in FeO 6 octahedra further splits the t 2g states into a 1g singlet and e g doublet; where θ = 2π/3. Note that, the t 2g triplet corresponds to an effective atomic orbital angular momentum l eff = 1 complex, a 1g and e g± to l eff n = 0 and ±1 eigenstates respectively, where n is the layer-normal direction. Finally, in the presence of ferromagnetic order parallel to n, as reported in BFPO, the spin degree of freedom is frozen and SOC takes the form of λl eff n S n , where the positive λ is the SOC magnitude of Fe d-orbital,l eff n is the effective atomic orbital angular momentum operator along n, and S n is magnitude of the FM order. Hence SOC under the FM order behaves as an atomic orbital Zeeman field to the e g± states making e g− state lower in energy than e g+ . Note that, without SOC these two atomic orbital states are degenerate, and the degeneracy is protected at Γ and K points even after hopping integrals are introduced to form the Bloch bands. The gap at Γ and K opens up when SOC is introduced.
The on-site Coulomb interaction enhances the gap between the e g+ and e g− Bloch states and gets bigger by the interaction parameter U . Due to the strong d-orbital Coulomb interaction in Fe, the system fully polarizes e g± orbitals, i.e., the both up spin e g− and e g+ orbitals are fully occupied and empty, respectively, as reported earlier in Ref. 28. Since the effective and real atomic orbital momenta are antiparallel to each other, the spin and the real orbital momenta at Fe add up to yield total magnetic moment of ∼ 5µ B , much larger than the size of d 6 high-spin moment S = 2µ B (assuming the g-factor ∼ 2). This is consistent to the value reported in experiment 27 , which confirms that the Hubbard U eff should be larger than 3eV, above which size of the total magnetic moment saturates close to the observed one.
Let us now construct an effective tight-binding model consisting of the e g+ orbital. Fig. 2 shows the schematic shape of e g+ orbital. First, n.n. hopping terms between the e g+ orbitals are real-valued due to the presence of n.n. inversion centers enforcing cancellation of complex phases. On the other hand, n.n.n. hopping channels can have complex values since they do not have any symmetry constraint. For example, let us consider one n.n.n. hopping channel between the e g+ orbitals at A-sublattice along the horizontal direction, as shown in lower figure in Fig. 2. Assuming only one off-diagonal hopping channel between the n.n.n. t 2g orbitals is active (between d yz and d xz orbitals, highlighted in blue and red colors in Fig.  2 respectively), in the e g+ subspace it yields complex hopping term t 2 e +2iθ , where the phase 2θ originates from e ±iθ assigned to d yz and d xz orbitals respectively. The presence of the n.n. inversion centers and additional three-fold symme- FIG. 2. (Color online) Illustration of a e g+ orbital at Fe site and depiction of next-nearest-neighbor (n.n.n.) complex hopping terms between the e g+ orbitals. In the hopping figure, dxz,yz orbitals, which contribute most to the horizontal n.n.n. hopping, are colored within each e g+ state. Note that, the n.n.n. hopping channels represented as dashed blue arrows can be transformed to each other by the inversion at the n.n. bond center and the threefold rotation symmetries.
try then generates all other n.n.n. channels; t 2 e +2iθ term for both A and B sublattices in a counterclockwise direction as shown in Fig. 2. For e g− orbitals, on the other hand, t 2 e −2iθ terms for counterclockwise n.n.n. hopping channels are obtained. The complex phase can deviate from 2θ depending on details of the t 2g hopping channels, but in general it does not vanish. Hence all of the conditions for realizing THM are fulfilled in BFPO, so that the system consists of two sets of THM (e g+ -and e g− -THM) with opposite chiralities to each other. As shown in Fig. 1, the e g− bands are fully occupied while e g+ bands are empty and mixed with a 1g band. Hence BFPO is in a trivial ferromagnetic Mott insulator phase, as the Fermi level lies in between two sets of Haldane bands with opposite Chern numbers.
The key ingredients for realizing the THM are spinpolarization due to strong Hund's coupling, and a separation of e g+ and e g− orbitals introduced by the broken TRS and SOC, as described above. Then, within the e g+ (or e g− ) subspace the phase factor of e ±2iθ emerges from the complex nature of the orbital wavefunction. To confirm our arguments above, DFT calculations incorporating Coulomb interactions and SOC are carried out. DFT+U 29 and Heyd-Scuseria-Ernzerhof (HSE) hybrid functional formalisms 30 are employed, which yield consistent results to each other when the value of U eff ≡ U − J parameter for Fe d-orbital in DFT+U computations is about 4 ∼ 5 eV 31 . Fig. 3(a) shows the DFT+U (U eff = 3 eV) bands for single-layer BFPO (SL-BFPO). First we comment that, while the bulk unit cell contains three Fe 2 (PO 4 ) 2 layers, the band splitting due to interlayer coupling is negligibly small as shown in Supplementary Material. This manifests the quasi-two-dimensional (2D) electronic structure of this compound, possibly one of the closest to the 2D limit among the other quasi-2D layered compounds ever synthesized. Hence in Fig. 3  shown. The spin is oriented along the n-direction, and all the down spin states are located below -2 eV, as shown in the projected density of states (PDOS) plot. About 1.5 eV below and 2 eV above the Fermi level, the up spin e g− and e g+ bands, respectively, are located. Chern numbers for the e g− bands are ±1 (for bulk ±3, with each layer contributing ±1).
While BFPO is a failed QAH though described by the THM, this system offers us an insight. One can obtain a halffilled e g± bands by adding one electron or hole per a formula unit, and this can be achieved by substituting Ba into alkali or lanthanide elements of similar ionic radii with Ba 2+ , say K 1+ and La 3+ for alkali and rare-earth ions, respectively 32 . Since BFPO is a layered compound with Ba layers residing between the Fe 2 (PO 4 ) 2 layers, replacing Ba into other cations can be done by intercalation of substitute cations or by using thin-film growth technique.
Next, the results for LFPO and CFPO are presented. Structural optimizations within R3 symmetry for both CFPO, KFPO, and LFPO with including van der Waals functionals were done to obtain the lattice parameters and internal coordinates, where the optimized structures are shown in Supplementary Materials. Note that, since results for CFPO and KFPO are almost identical to each other, except the c parameter, here we only show results from CFPO. Fig. 3(b) and (c) shows the LFPO bulk and zigzag-edge bands respectively, dominated by the e g+ states (U eff = 4 eV). The occupied and unoccupied bands have Chern number C = -1 and +1 per a Fe 2 (PO 4 ) 2 layer, respectively, with a well-defined bulk gap of ∼ 400 meV and showing a chiral edge state at one zigzag edge side as shown in Fig. 3(c). On the other hand, Fig. 3(d) and (e) shows the CFPO bands, showing e g− state character with opposite chirality. The bulk gap is about 220 meV, which is smaller than that of LFPO but still substantial. The hopping integrals for the e g± states in both compounds are obtained from the Wannier orbital calculations, where the values are shown in Table I. Note that they are well described by the THM. Unlike in 4d-or 5d honeycomb materials, such as α-RuCl 3 or A 2 IrO 3 (A=Li,Na) having the similar edgesharing octahedral structure, AFPO show almost negligible third-nearest neighbor hopping terms due to the spatially localized 3d orbitals 33 . t 3 term in LFPO is larger than that of CFPO due to the smaller in-plane lattice constant, and further enhancing t 3 with epitaxial strain may induce a phase transition from C = ±1 to ±2 phase as reported in α-RuCl 3 or A 2 IrO 3 33 . It should be commented that, although the size of the BFPO band gap in DFT+U and HSE results are well matched at U eff 5 eV, our results for all of AFPO systems in this work remain robust in a wide range of U eff value, between 1.8 and 7 eV. In CFPO (KFPO) and LFPO, a bipartite charge ordering which breaks the topological phase is not observed across the U eff parameter range we tested. These observations suggest the robustness of THM and the resulting topological phase in these systems.
Discussions.-Our result suggests a new strategy in realizing THM phase with sizable gap in condensed matter systems. Previously, there have been roughly two different approaches; depositing magnetic ions or interfacing magnetic systems onto graphene or other honeycomb lattices [34][35][36] , and planting magnetic ions in topological insulators 11,18,19 . We find that searching for systems with the ferromagnetic order in transition metal compounds with strong Hund's coupling and SOC is a promising way for the realization of THM. In addition to AFPO investigated in this study, there has been a report of possible quasi-2D Ising ferromagnets in several 3dtransition metal halides 37,38 . Since they have a t 2g orbital degree of freedom with spin splitting larger than the cubic crystal field, such system can show a similar atomic-orbital formation which may possibly lead to the formation of THM. Heterostructures of transition metal oxides with magnetic ions, such as a double perovskite Ba 2 FeReO 6 39 , can be another candidate.
Lastly we comment on the size of the gap in BFPO. In the paramagnetic high-temperature R3 phase, activation en-ergy estimation from the conductivity is about 0.2 eV 40 , while optical spectroscopy estimates the gap size to be about 1.5 eV 40 , which are significantly smaller than our results of ∼3 eV from the hybrid functional calculation. More reliable lowtemperature data should be measured for further studies, and we emphasize that our results is robust independent of the value of U and corresponding gap, as long as U eff is larger than 1.8eV. A further theoretical studies incorporating quantum fluctuation effects for these systems when a fraction of electron or hole is added would be interesting, which may reveal fractional Chern insulator phases in these systems 41 .
In summary, we propose a way to search for realistic materials described by the THM: effectively spinless fermion with complex n.n.n hopping integrals can be found in ferromagnetic insulators with strong Hund's coupling and finite SOC. We apply this idea to a series of Fe-based honeycomb ferromagnetic oxides. We show that AFPO series is represented by THM, and predict that honeycomb CFPO, KFPO, and LFPO are candidates for the original THM exhibiting a quantum Hall effect without Landau levels. Our study provides a platform for exploring correlated topological materials, including possible fractional Chern insulators via strong correlations.
Methods.-DFT computations in this work are done with Vienna ab-initio Simulation Package (VASP) 42 For the structural optimization, we employ the Vienna abinitio Simulation Package (VASP), which uses the projectoraugmented wave (PAW) basis set 42,43 . 520 eV of plane wave energy cutoff is used, and for k-point sampling 15×15×15 grid including Gamma point is employed for the rhombohedral primitive cell. A revised Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (PBEsol) 46 is used for structural optimization and total energy calculations, which yields the best agreement of calculated lattice parameters compared to local density approximation or other GGA functionals.
Structural optimizations using VASP are carried out in two stages. First, optimizations of cell parameters and internal coordinates are performed under the R3 symmetry constraint and in the absence of SOC. After that, the layer-normal clattice parameter in Fig. 4 is again optimized with using vdW-optB86b van der Waals functional 47 and with fixed in-plane a-lattice parameter. vdW-optB86b results yield 2% smaller and 0.5% larger value of c-parameter for KFPO and LFPO, respectively, compared to PBEsol-optimized c values. Force criterion of 1 meV /Å is used for internal coordinates and stress tensor optimizations.
Electronic structure calculations with including Coulomb interactions are both performed with using VASP and OPENMX 44 codes. In both codes, Dudarev's rotationally invariant DFT+U formalism 29 is employed with effective U eff ≡ U − J increased up to 6 eV for Fe d orbitals and 8 eV for La f -orbitals. For hybrid functional calculations, the exchange-correlation functional of Heyd-Scuseria-Ernzerhof 30 implemented in VASP is employed, with the mixing parameter α = 0.25 and the inverse effective screening length ω = 0.2Å −1 . For the HSE calculations, to reduce the computational cost, two-dimensional unit cells with vacuum of 20Å and 11×11×1 of k-grid sampling are employed for each system. Maximally localized Wannier orbital formalism 48,49 as implemented in OPENMX 50 is employed for computations of the e g± Wannier orbitals.   Table II presents the optimized structure for AFPO systems in the R3 symmetry, compared to the experimentally reported structures. Cell parameter optimizations for BFPO yields 0.9% smaller and 1.3% larger a and c parameters, respectively, compared to the experimentally observed ones at T = 1.8K. Optimized BFPO internal coordinates show fairly good agreement to the observed one at T = 293K. FeO 6 octahedral distortion is slightly larger compared to the roomtemperature structure. In KFPO and LFPO, the c parameter is increased and decreased by 7% compared to BFPO, respectively, reflecting the ionic character of the interlayer bonding mediated by the K 1+ , Ba 2+ , or La 3+ cations. Compared to KFPO, CFPO shows 10% larger c parameter due to the larger ionic radius of Cs 1+ compared to K 1+ , but without significant differences in the a parameter size and internal coordinates. KFPO and CFPO shows smaller octahedral distortion and Fe buckling than BFPO, with small mismatch of the in-plane lattice constant compared to the BFPO in-plane lattice constant (0.5% and 0.1% for KFPO and CFPO, respectively). LFPO, on the other hand, show 2% smaller in-plane lattice parameter, which is attributed to the completely filled three t 2g bonding orbitals as shown later in Fig. 5(g). The smaller value of c parameter in LFPO yields increased interlayer hopping terms compared to BFPO, CFPO, and KFPO, but not enough to break the quasi-two-dimensional picture.

III. SUPPLEMENTARY MATERIAL C: ELECTRONIC STRUCTURE RESULTS
Band structure and its evolution with respect to the inclusion of U in BFPO is presented in Ref. 28  comment some points which weren't mentioned in previous studies. First, the coupling between the different Fe 2 (PO 4 ) 2 is very weak in BFPO. Fig. 5(a) shows the bulk BFPO bands in a hexagonal unit cell with three Fe 2 (PO 4 ) 2 layers (grey lines) overlaid with a single-layer (SL) BFPO bands (black solid and dashed lines). Band splitting due to interlayer coupling is almost absent in the bulk bands, smaller then the linewidth to draw the SL bands in the plot. As shown later, this quasi-2D character is maintained in CFPO, KFPO, and BFPO, as mentioned in the previous section. Next, the three bands near the Fermi level in the absence of SOC and U , enclosed in a shaded box in Fig. 5(a), consist of t 2g bonding-orbitals centered at each n.n. bond center. This bonding-orbital formation is driven by a strong σ-like dd direct overlap between the two n.n. t 2g orbitals, where the overlap integral is about -0.45 eV and overcomes other t 2g hopping channels. The prevailing σ-overlap occurs due to the edge-sharing geometry of FeO 6 octahedra yielding rather short n.n. distance. Hence BFPO in the weakly correlated limit forms a bonding-orbital kagome lattice, as suggested in Ref. 51, but with significant n.n.n. and n.n.n.n. hopping terms. Degeneracies at Γ and K points are protected by the complex conjugation K and the product of inversion and K respectively, where K is a TRS operation followed by a global spin rotation reorienting the spin moments to the original direction. SOC breaks this pseudo-TRS and opens the gap at both points, as shown in Fig. 5(b).
As mentioned in the main text, the local atomic orbital picture survives in the strongly correlated regime. Hence a crossover happens in a intermediate strength of U (or equivalently U eff ). Fig. 5(c-e) shows the crossover, where the transition from a Chern to a trivial Mott insulator happens at U c eff = 1.3 eV in our calculations 45 . This critical value is smaller than the one reported in Ref. 28, U c − J = 2.5 -0.7 = 1.8 eV, due to the difference in the choice of correlated orbital projectors. The loss of the Chern insulator phase and the onset of the large orbital moment formation is a signature of the e g± -polarized atomic orbital picture. Just above the U c eff , the lowest unoccupied band carries a zero Chern number, however it crosses again with other t 2g bands and obtains a nontrivial Chern number in higher U eff value. Note that, in Ref. 45, the bulk Chern number in the weak-U regime was -3. While the sign difference originates from the direction of the magnetic moments, the difference in the Chern number magnitude merely comes from different definition. In Ref. 45, the total Chern number (per unit cell) is used, while in our paper it is the Chern number per a FePO 4 layer. The Chern number magnitude in Ref. 45 is consistent with ours (±1 per a FePO 4 layer) since the bulk primitive unit cell with rhombohedral symmetry employed in Ref. 45 implies the presence of three FePO 4 layers in the unit cell. Alternatively a conventional hexagonal unit cell including three FePO 4 layers can be used, and both choices of bulk unit cell yield ±3 (sign depending on the direction of the magnetic moment) of total Chern number in the weak-U eff regime. The DFT+U eff band gap agrees best with the HSE one in the range of U eff = 4.5 to 5 eV, and we chose 5 eV for presentation of the results. HSE calculation reproduces the formation of e g± -polarized band structure, as shown in Fig. 5(f). Compared to DFT+U eff , HSE bands show larger exchange splitting between the e g− and the fully occupied spin states, yielding well-separated e g− THM bands from the others. Also, the s-like bands at the bottom of conduction bands are pushed above in the HSE result. Now the bulk LFPO ( Fig. 5(g,h)) and CFPO (Fig. 5(i,j)) results are presented. Note that, KFPO bands are almost identical to those from CFPO and are not presented here. In LFPO, in the absence of SOC and U eff , the t 2g bonding orbital bands are fully occupied, as shown in Fig. 5(g). This induces the smaller in-plane lattice constant in LFPO compared to the others. Note that, the band splitting due to the interlayer coupling is visible, although not significant, due to the smaller interlayer distance here. Inclusion of U eff = 4 eV transforms the orbital landscape from the bonding-antibonding picture to the atomic one, and in Fig. 5(h) the half-filled e g+ THM bands are clearly seen. Interestingly, HSE result in Fig. 5(h) shows larger e g+ THM bandwidth with larger band gap of ∼ 0.8 eV at the K point, with much better separated e g+ THM bands compared to the DFT+U eff result. Similarly, CFPO HSE result in Fig. 5(j) shows much larger separation between the e g+ THM bands the lower ones, and the larger bandwidth and gap size of the THM bands. These enhancement of gap sizes and bandwidths, compared to the DFT+U eff results, might be attributed to the nonlocal correlation effect inherent in the hybrid functional approach but absent in DFT+U eff .