Kitaev interactions between j=1/2 moments in honeycomb Na2IrO3 are large and ferromagnetic: insights from ab initio quantum chemistry calculations

Na$_2$IrO$_3$, a honeycomb 5$d^5$ oxide, has been recently identified as a potential realization of the Kitaev spin lattice. The basic feature of this spin model is that for each of the three metal-metal links emerging out of a metal site, the Kitaev interaction connects only spin components perpendicular to the plaquette defined by the magnetic ions and two bridging ligands. The fact that reciprocally orthogonal spin components are coupled along the three different links leads to strong frustration effects and nontrivial physics. While the experiments indicate zigzag antiferromagnetic order in Na$_2$IrO$_3$, the signs and relative strengths of the Kitaev and Heisenberg interactions are still under debate. Herein we report results of ab initio many-body electronic structure calculations and establish that the nearest-neighbor exchange is strongly anisotropic with a dominant ferromagnetic Kitaev part, whereas the Heisenberg contribution is significantly weaker and antiferromagnetic. The calculations further reveal a strong sensitivity to tiny structural details such as the bond angles. In addition to the large spin-orbit interactions, this strong dependence on distortions of the Ir$_2$O$_2$ plaquettes singles out the honeycomb 5$d^5$ oxides as a new playground for the realization of unconventional magnetic ground states and excitations in extended systems.


Introduction
The Heisenberg model of magnetic interactions, JS i · S j between spin moments at sites {i, j}, has been successfully used as an effective minimal model to describe the cooperative magnetic properties of both molecular and solid-state many-electron systems. A less conventional spin model -the Kitaev model [1] -has been recently proposed for honeycomb-lattice materials with 90 • metal-oxygen-metal bonds and strong spin-orbit interactions [2]. It has nontrivial topological phases with elementary excitations exhibiting Majorana statistics, which are relevant and much studied in the context of topological quantum computing [1,3,4,5,6,7]. Candidate materials proposed to host such physics are the honeycomb oxides Na 2 IrO 3 and Li 2 IrO 3 [2]. The magnetically active sites, the Ir 4+ species, display in these compounds a 5d 5 valence electron configuration, octahedral ligand coordination and bonding of nearest-neighbor (NN) Ir ions through two ligands [8,9]. In the simplest approximation, i.e., for sufficiently large t 2g -e g octahedral crystal-field splittings within the Ir 5d shell and degenerate Ir t 2g levels, the ground-state electron configuration at each Ir site is a t 5 2g effective j = 1/2 spin-orbit doublet [10,11,12,2]. The anisotropic, Kitaev type coupling then stems from the particular form the superexchange between the Ir j = 1/2 pseudospins takes for 90 • bond angles on the Ir-O 2 -Ir plaquette [2,13,14]. Recent measurements on Na 2 IrO 3 [8,9] indicate significant lattice distortions away from the idealized case of cubic IrO 6 octahedra and 90 • Ir-O-Ir bond angles for which the Kitaev-Heisenberg (KH) model was proposed [2,13]. Lower-symmetry crystal fields and distortions of the Ir-O-Ir bonds obviously give rise to finite Ir t 2g splittings [15,16] and more complex superexchange physics [17,18]. It has been actually shown that the interplay between "local" distortions of the O cage and longer-range crystal anisotropy is a key feature in 5d oxides [19,20,21,22,23] and the outcome of this competition is directly related to the precise nature of the magnetic ground state [23]. Moreover, the lower symmetry characterizing a given [Ir 2 O 10 ] unit of two edge-sharing octahedra allows in principle for nonzero anisotropic interaction terms beyond the Kitaev picture.
The inelastic neutron scattering data [8] and the magnetic ordering pattern [9,24] in Na 2 IrO 3 could in principle be explained in a minimal model by either i) a more conventional Heisenberg model with frustrated exchange couplings extending up to third NN's [8] or ii) a KH model with dominant antiferromagnetic (AF) Kitaev and smaller ferromagnetic (FM) Heisenberg NN terms [14]. The presence of strong Kitaev interactions has furthermore been suggested on the basis of recent resonant inelastic x-ray scattering experiments [25].
To clarify the signs and the strengths of the effective coupling constants in Na 2 IrO 3 , we here employ many-body ab initio techniques from wave-function-based quantum chemistry [26]. We establish that the situation is much more subtle than presumed so far. In model systems the KH Hamiltonian arises due to the destructive interference of different superexchange pathways that contribute equally to the effective intersite interaction [2,13]. This interference turns out to be rather fragile as even when we consider idealized structures with cubic IrO 6 octahedra, orthogonal Ir-O-Ir bonds and D 2h point-group symmetry of the Ir-Ir link, the computed low-energy magnetic spectrum does not support a pure KH model. A careful analysis shows that in the Kitaev reference frame [2,13] off-diagonal terms of the symmetric anisotropic exchangecoupling tensor are allowed by symmetry to be non-zero. The quantum chemistry (QC) calculations predict the latter are comparable in magnitude to the strength of the isotropic Heisenberg term. We find that the effective Kitaev coupling is, however, the dominant energy scale -depending on geometrical details in the range of 10-20 meV and FM -and that the NN Heisenberg J is AF and significantly weaker. For NN interaction parameters as derived in the QC study, we have further performed exact diagonalization (ED) calculations including additionally finite AF second and third order Ir-Ir Heisenberg couplings. These indicate the presence of zigzag AF order, in agreement with the experimentally observed spin texture [8,9,24].

Results and discussion
Multiconfiguration complete-active-space self-consistent-field (CASSCF) and multireference configuration-interaction (MRCI) calculations [26] were performed on embedded clusters made of two reference IrO 6 octahedra (for technical details, see Appendix A). All possible occupations were allowed within the set of t 2g orbitals at the two magnetically active Ir sites in the CASSCF calculations. The orbitals were optimized for an average of the lowest nine singlets and the nine triplet states. All these singlet and   Table 1.
The (x, y, z) coordinate frame used to express the KH Hamiltonian [2] is also drawn for one of the Ir-Ir links. For each Ir 2 O 2 plaquette in the actual C2/m structure [8], due to trigonal squashing of the IrO 6 octahedra (normal to the Ir honeycomb plane), the apical-like Ir-O bonds are not along the corresponding Kitaev axis.  triplet states entered the spin-orbit calculations, both at the CASSCF and MRCI levels.
In MRCI, single and double excitations from the Ir t 2g shells and the 2p orbitals of the bridging ligands were accounted for. A similar strategy of explicitly dealing only with selected groups of localized ligand orbitals was earlier adopted in QC studies on both 3d [27,28,29,30] and 5d [20,21,23] compounds, with results in good agreement with the experiment [20,21,27,28,29,30]. For a pair {i, j} of magnetic sites in systems in which the midpoint of the ij link displays inversion symmetry, the most general bilinear exchange Hamiltonian is whereS i ,S j are pseudospin operators (S = 1/2) [2,13] and the elements Γ αβ form a traceless symmetric second-rank tensor. It is convenient to choose the X axis along the Ir-Ir link and Z perpendicular to the plaquette defined by the two Ir ions and the two bridging ligands because in the C2/m crystal structure of Na 2 IrO 3 , for Ir-Ir bonds along b, see Fig. 1, the Ir-Ir axis is a C 2 axis with an orthogonal mirror plane [8,9], i.e., the symmetry of those [Ir 2 O 10 ] units is C 2h . With such a choice of the coordinate system only Γ Y Z = Γ ZY are finite and so ( The fact that Y and Z are not C 2 axes is related to the configuration of the four adjacent Ir sites -two of those are below and two above the XY plane, with no inversion center -and the trigonal squashing of the IrO 6 octahedra [8]. The KH Hamiltonian is however expressed in a (x, y, z) coordinate frame that has the (x, y) coordinates rotated by 45 • about the Z = z axis [2,13], as compared to (X, Y ), see Fig. 1, and Γ then becomes (see Appendix B for details). For a more transparent picture and better insight into the nature of the NN magnetic couplings, it is instructive to first consider two-octahedra clusters taken from an idealized crystalline model without trigonal distortions and with all adjacent Ir and Na sites modeled as identical point charges. In this case, the overall symmetry is D 2h and all off-diagonal couplings cancel by symmetry, in the (X, Y, Z) coordinate system with X along the Ir-Ir link. For an idealized [Ir 2 O 10 ] unit displaying D 2h symmetry C = 0 and the spin Hamiltonian reduces to . The off-diagonal xy coupling, last term in (4), is allowed by symmetry even for ideal octahedra at 90 • Ir-O-Ir bonding but has been neglected in earlier studies on Na 2 IrO 3 [2,13,14,31,16].
Results of spin-orbit calculations, both at the CASSCF (CAS+SOC) and MRCI (MRCI+SOC) levels, are listed for idealized [Ir 2 O 10 ] D 2h model clusters in Table 1. Such a cluster is highly charged, 12−. To ensure charge neutrality, we assigned to each of the 26 adjacent Na and Ir sites fictitious point charges of +12/26. In the simplest approximation, no farther embedding was used for these D 2h clusters. The ab initio calculations were performed for both i) regular IrO 6 octahedra and 90 • Ir-O-Ir bond angles and ii) distorted geometries with all ligands in the xy plane pushed closer to the Ir-Ir axis and therefore larger Ir-O-Ir bond angles but keeping the D 2h bond symmetry. The Ir-Ir distance d(Ir-Ir) and in the latter case the Ir-O-Ir angle were set to 3.133Å and 98.5 • , respectively, average values in the C2/m crystal structure reported in [8].
To determine the nature of each spin-orbit state we explicitly compute the dipole and quadrupole transition matrix elements among those four low-lying states describing the magnetic spectrum of two edge-sharing octahedra. A careful symmetry analysis reveals that the spin-orbit wave functions Ψ S , Ψ 1 , Ψ 2 and Ψ 3 defined in Table 1 transform according to the A g , B 2u , B 1u and A u irreducible representations, respectively. Standard selection rules and the nonzero dipole and quadrupole matrix elements in the QC outputs then clearly indicate which state is which. We also carried out the transformation of the spin-orbit wave functions from the usual {L 1 ,M L 1 ,L 2 ,M L 2 ,S,M S } basis in standard QC programs to the {S 1 ,S 2 ,M S 1 ,M S 2 } basis. This allows the study of Ψ 1 -Ψ 2 mixing when the point-group symmetry is reduced to C 2h , see below. Having the assignment of the states resolved, the Ψ S -Ψ 1 splitting provides J, the Ψ 2 -Ψ 3 splitting yields D, while Table 2. Energy splittings of the four lowest magnetic states and effective coupling parameters (meV) for two NN IrO 6 octahedra in the C2/m structure of Ref. [8]. The weight of (↑↓ + ↓↑)/ √ 2 and (↑↑ + ↓↓)/ √ 2 in Ψ ′ 1 and Ψ ′ 2 , respectively, is ≈ 98%, see text.

Method
the difference between the energy of Ψ 1 and the average of the E 2 (Ψ 2 ) and E 3 (Ψ 3 ) eigenvalues equals −K/2, see Appendix B. The QC data in Table 1 indicate AF J's, FM K's and off-diagonal anisotropic couplings comparable in strength to the isotropic J interaction. Interestingly, the ab initio MRCI calculations indicate a much stronger K for nonorthogonal Ir-O-Ir bonds. This shows that deviations from rectangular geometry on the Ir 2 O 2 plaquette is not a negligible factor, as presently assumed in simplified superexchange models for Kitaev physics in honeycomb Na 2 IrO 3 [2,13,14]. The effect of the MRCI treatment is also stronger for nonorthogonal Ir-O-Ir bonds: the CAS+SOC K and D coupling parameters are enlarged by more than 50% by including O 2p to Ir 5d charge-transfer effects, Ir t 2g to e g excitations and additional correlations accounted for in MRCI while J changes sign. This strong enhancement of the FM K for nonorthogonal Ir-O-Ir bonds and the tiny effect of the MRCI treatment on the FM K for rectangular geometry also disagrees with predictions of present approximate superexchange models that indicate the t 2g to e g excitations and hopping as the dominant superexchange mechanism, giving rise to a large AF K [14].
Relative energies and the resulting effective coupling constants are next given in Table 2 for the experimentally determined C2/m crystal structure of Ref. [8]. For this set of calculations we used effective embedding potentials as described in Appendix A. There are two inequivalent Ir-Ir links in Na 2 IrO 3 , displaying different Ir-O-Ir bond angles and slightly different Ir-O and Ir-Ir distances [8]. While the [Ir 2 O 10 ] block with larger Ir-O-Ir bond angles (upper part in Table 2) displays C 2h symmetry, for the other unit of edge-sharing octahedra the point-group symmetry is even further reduced to C i (lower part in Table 2). The expressions of the spin-orbit wave functions in the transformed {S 1 ,S 2 ,M S 1 ,M S 2 } basis show, however, that the mixing of the Ψ i terms as expressed in the idealized D 2h geometry is negligible in the C2/m structure. Therefore the ab initio data is mapped also in this case on the effective model described by (4).
As for the idealized D 2h configuration, the MRCI+SOC results indicate large FM Kitaev couplings, weaker AF Heisenberg superexchange and sizable D anisotropic interactions. The latter are not included in the plain KH model [2,13,14,31] while the signs of K and J that we compute are different from those proposed in the recent model-Hamiltonian analysis of Ref. [14]. We note that in agreement with our findings, relatively large FM Kitaev couplings K have been earlier predicted by Kimchi and You [31] from the analysis of the phase diagram obtained by ED on modest size clusters and by Foyevtsova et al. [16] on the basis of an effective superexchange model fed with electronic-structure parameters obtained from density-functional calculations for the same C2/m structure [8]. However, the NN Heisenberg J is also FM in the latter work, different from the small AF values we find in the MRCI calculations. We also find that on each hexagonal Ir 6 unit the two Ir-Ir links along the b-axis have effective coupling constants significantly different from the set of parameters associated with the other four Ir-Ir "bonds" due to subtly different oxygen distortions. Together these findings stress the importance of lattice distortions and symmetry issues and lay the foundation for rigorous ab initio investigations of unusually large anisotropic interactions such as the Kitaev exchange in strongly spin-orbit coupled 5d oxides ‡.
It is known experimentally that Na 2 IrO 3 displays zigzag AF order at low T 's [8,9,24]. It has been also argued that the longer-range magnetic interactions, up to the second and third Ir coordination shells, are sizable and AF [8,31,35,16]. We therefore performed ED calculations for a KH model supplemented with second and third NN couplings J 2 and J 3 (see Appendix B), on a 24-site cluster with periodic boundary conditions as used in earlier studies [13,14]. We disregarded the presence of two structurally and magnetically different sets of Ir-Ir links and on the basis of the QC results of Table 2, used on all bonds J, K and D coupling constants of 3, −17.5 and −1 meV, respectively (approximately averaged over all bonds). For a given set of J 2 and J 3 values the dominant order is determined according to the wave number Q = Q max giving a maximum value of the static structure factor S(Q). The resulting phase diagram, see Fig. 2, shows that the zigzag phase is indeed stable in the region of J 2 ,J 3 2 meV. We note that positive J 2 and J 3 values of 4-5 meV would be consistent with the experimentally observed Curie-Weiss temperature ≈ −125 K [35] using θ CW = −S(S + 1)(J + 2J 2 + J 3 + K/3)/k B [8]. Thus we propose that an extended spin Hamiltonian based on the nearest-neighbor anisotropic exchange terms found from the ab initio QC calculations supplemented by further-neighbor exchange integrals could provide a realistic starting point to explain the magnetism of Na 2 IrO 3 . ‡ Detailed QC studies of the anisotropic terms have been so far confined to 3d oxides, where the spin-orbit interaction is just a small perturbation and the anisotropic coupling parameters are orders of magnitude weaker than the isotropic Heisenberg exchange [32,33,34].
To discuss in more detail the generic phase diagram in Fig. 2, we note that in between the ordered phases we find a narrow spin liquid phase characterized by Q max ≃ 0 and very weak spin-spin correlations. We checked that this phase is adiabatically connected to the Kitaev liquid phase appearing for larger |K| [13,14]. The canted "stripe I" phase corresponds to the canted "phase III" found in earlier investigations of the isotropic J 1 -J 2 -J 3 model [36,37,38]. Interestingly, a finite D enlarges the extent of both the spin liquid and the stripe I phases, showing that the planar spin fluctuations effectively amplify frustration in the model. The ED results show that the phase diagram is also very sensitive to the longer-range exchange couplings J 2 and J 3 . These findings are relevant in the context of recent experimental data that indicate a qualitatively different AF ground state for the related compound Li 2 IrO 3 [39].

Conclusions
In sum, for the honeycomb iridate Na 2 IrO 3 , the ab initio quantum chemistry calculations show that in a reference system with X along the Ir-Ir link and Z perpendicular on the Ir 2 O 2 plaquette the X-Y anisotropy is significant and gives rise in the rotated (x, y, z) Kitaev-type frame [2,13] to off-diagonal anisotropic terms beyond the plain Kitaev-Heisenberg model. Nevertheless, the calculations predict that the largest energy scale is the Kitaev interaction, 10 to 20 meV, while the NN Heisenberg superexchange and the off-diagonal xy coupling are significantly weaker. The quantum chemistry data additionaly establish that the Kitaev term is FM. Further, all NN couplings are highly sensitive to subtle distortions involving the O ions. This makes the material dependence along the Na 2 IrO 3 , Li 2 IrO 3 and Li 2 RhO 3 series an interesting topic for future investigations. Large variations of the effective couplings as function of bond lengths and bond angles and a variable degree of "inequivalence" of those sets of parameters for structurally distinct Ir-Ir links in the honeycomb layer may in principle give rise to very different types of magnetic ground states in different 5d 5 or 4d 5 honeycomb compounds. Strong exchange anisotropy is also a topic of active research in the field of molecular magnetism and single-molecule magnets. The focus there has so far been on 3d and 4d compounds [40,41,42] but clearly 5d systems with stronger spin-orbit couplings may provide new playgrounds in this research area too. After submission of our paper, results of ED calculations including off-diagonal D type terms were also reported in Ref. [43].

Appendix A: Computational details in the QC calculations
For the computation of the intersite spin couplings, two reference NN IrO 6 octahedral units are considered. Since it is important to describe the finite charge distribution at sites in the immediate neighborhood [44,45], the closest 22 Na neighbors and the other four, adjacent octahedra are also explicitly included in the actual cluster. However, to make the analysis of the low-lying magnetic states tractable, we cut off the spin couplings with the adjacent 5d ions by replacing those open-shell Ir 4+ 5d 5 NN's with closed-shell Pt 4+ 5d 6 species. This is an usual procedure in quantum chemistry studies on transition-metal systems, see, e.g., Refs. [46,20,21,22,23,45,47], and here allows a straightforward mapping of the ab initio data onto the effective spin model. It has been shown in earlier work that this way of modeling the NN 5d ions does not affect the size of the t 2g splittings at the central, reference site [20,21] and that the computed Heisenberg couplings agree well with estimates derived from experiment [20,21]. The surrounding solid-state matrix is described as a finite array of point charges fitted to reproduce the crystal Madelung field in the cluster region. Most of the ab initio calculations were carried out with the molpro quantum chemistry package [48]. Test calculations in finite magnetic fields were also performed with the orca program [49].
We used energy-consistent relativistic pseudopotentials with quadruple-zeta basis sets for the valence shells of the two reference Ir ions [50], all-electron quintuple-zeta basis sets for the bridging ligands [51] and triple-zeta basis functions for the other O's of the two reference octahedra [51]. For the NN Pt 4+ 5d 6 species triple-zeta basis sets were applied [50] while the other O ions not shared with the central octahedra were modeled with minimal atomic-natural-orbital basis functions [52]. All occupied shells at the NN Na + sites were incorporated in the large-core pseudopotentials and each of the Na 3s orbitals was described with a single basis function [53]. For the central Ir ions and the two bridging ligands we also employed polarization functions, two Ir f and four d O functions [50,51]. To separate the metal 5d and O 2p valence orbitals into different groups, we used the orbital localization module available in molpro. The MRCI calculations were carried out for each spin multiplicity, singlet or triplet, as nineroot calculations. Only the four low-lying spin-orbit states are relevant for the analysis of the NN magnetic interactions. The higher-lying spin-orbit states imply an excitation energy of at least 0.6 eV. This gap concerns the j = 1/2 to j = 3/2 transitions [15].
The calculations in finite magnetic fields performed with the orca package were used to crosscheck the assignment of the lowest four magnetic states made on the basis of the dipole/quadrupole transition matrix elements, symmetry analysis and selection rules. These test calculations were carried out only for idealized D 2h structural models such as those of Table 1. It can be shown that when B Oz, B Oy or B Ox, state Ψ 1 , Ψ 2 or Ψ 3 , respectively, defined as in Table 1, should not change energy. This provides a quick alternative way of identifying the low-lying spin-orbit states.

Appendix B: Lattice spin model
In the main body of this article, we made use of two different local reference frames, for each of the three bond types b = 1−3, see Fig. 3. Choosing the global frame to be {x, y, z} ≡ {x 1 , y 1 , z 1 }, the local We note that the system is invariant under a two-fold rotation around X 1 and that the directions of X 2 and X 3 are chosen to map to each other under this rotation (otherwise the sign of the exchange coupling parameter C is arbitrary).
To be explicit, let us write down all NN interaction terms for each of the three types of bonds b in the corresponding reference frame where . We can rewrite these terms in the global frame {x, y, z} using (5): For Na 2 IrO 3 bonds 2 and 3 are equivalent and therefore J 2 = J 3 , K 2 = K 3 , D 2 = D 3 , C 2 = C 3 . In C 2h symmetry, D b = 0 and for a given bond b the eigenvalues of the Hamiltonian defined by Eqn. 6 are These expressions become more complicated for point-group symmetries lower than C 2h . The analysis of the spin-orbit wavefunctions in the transformed {S 1 ,S 2 ,M S 1 ,M S 2 } basis shows, however, that the mixing of the D 2h Ψ i eigenvectors for distorted clusters and lower symmetries of the [Ir 2 O 10 ] blocks (see Table 2) is negligible in the C2/m structure determined by Choi et al. [8].
For D 2h symmetry (see Table 1