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

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–7]. Candidate materials proposed to host such physics are the honeycomb oxides Na₂IrO₃ and Li₂IrO₃ [2]. The magnetically active sites, the Ir⁴⁺ species, display in these compounds a 5d⁵ 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 (GS) electron configuration at each Ir site is a t⁵_2g effective j = 1/2 spin–orbit doublet [2, 10–12]. 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₂–Ir plaquette [2, 13, 14].

Recent measurements on Na₂IrO₃ [8, 9] indicate significant lattice distortions away from the idealized case of cubic IrO₆ 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–23] and the outcome of this competition is directly related to the precise nature of the magnetic GS [23]. Moreover, the lower symmetry characterizing a given [Ir₂O₁₀] 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₂IrO₃ 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 NNs [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₂IrO₃, we here employ many-body ab initio techniques from wave-function-based quantum chemistry (QC) [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₆ 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 exchange-coupling tensor are allowed by symmetry to be nonzero. The 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].

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₆ 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 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–30] and 5d [20, 21, 23] compounds, with results in good agreement with the experiment [20, 21, 27–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

$$\begin{equation} \mathcal{H}_{ij} = J_0\tilde{\boldsymbol{S}}_i\cdot \tilde{\boldsymbol{S}}_j + \displaystyle\sum\limits_{\alpha,\beta\in\{X,Y,Z\}} \Gamma_{\alpha\beta} {\tilde{S}}_i^{\alpha} {\tilde{S}}_j^{\beta}, \end{equation} \tag{ 1 }$$

where $\tilde {\boldsymbol {S}}_i$, $\tilde {\boldsymbol {S}}_j$ are pseudospin operators ($\tilde {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₂IrO₃, for Ir–Ir bonds along b, see figure 1, the Ir–Ir axis is a C₂-axis with an orthogonal mirror plane  [8, 9], i.e. the symmetry of those [Ir₂O₁₀] units is C_2h. With such a choice of the coordinate system only Γ_{Y Z} = Γ_ZY are finite and so

$$\begin{equation} \boldsymbol{\Gamma}=\left( \begin{array}{@{}ccc@{}} A & 0 & 0\\ 0 & B & C\\ 0 & C & -A-B\\ \end{array} \right)_{\!\!\{X,Y,Z\}}. \end{equation} \tag{ 2 }$$

The fact that Y and Z are not C₂ 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₆ 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 figure 1, and Γ then becomes

$$\begin{equation} \boldsymbol{\Gamma}=\left( \begin{array}{@{}ccc@{}} (A+B)/2 & (A-B)/2 & -C/\sqrt{2}\\ (A-B)/2 & (A+B)/2 & C/\sqrt{2}\\ -C/\sqrt{2} & C/\sqrt{2} & -A-B \end{array} \right)_{\!\!\{x,y,z\}} \end{equation} \tag{ 3 }$$

(see appendix B for details).

Zoom In Zoom Out Reset image size

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₂O₁₀] unit displaying D_2h symmetry C = 0 and the spin Hamiltonian reduces to

$$\begin{equation} \mathcal{H}_{ij}^{{D_{2{\rm h}}}} = J\,\tilde{\boldsymbol{S}}_i\cdot \tilde{\boldsymbol{S}}_j + K\,{\tilde{S}}_i^z {\tilde{S}}_j^z + D\left( {\tilde{S}}_i^x {\tilde{S}}_j^y + {\tilde{S}}_i^y {\tilde{S}}_j^x \right), \end{equation} \tag{ 4 }$$

where $K=-\frac {3}{2}(A+B)$, J = J₀ − K/3 and $D = \frac {1}{2}(A-B)$. 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₂IrO₃ [2, 13, 14, 16, 31].

Results of spin–orbit calculations, both at the CASSCF (CAS + SOC) and MRCI (MRCI + SOC) levels, are listed for idealized [Ir₂O₁₀] 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₆ 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].

Table 1. Energy splittings of the four lowest magnetic states and effective exchange parameters (meV) for two edge-sharing IrO₆ octahedra taken out of the idealized structural model with D_2h bond symmetry. In a pure KH model, Ψ₂ and Ψ₃ are degenerate. d(Ir–Ir) = 3.133 Å, see text.

Method	CAS+SOC	MRCI+SOC
$\measuredangle ({\mathrm {Ir}}$–O–Ir)=90°:
$\Psi _{\mathrm {S}}= (\uparrow \downarrow -\downarrow \uparrow )/\sqrt 2$	0.0	0.0
$\Psi _{\mathrm {2}}= (\uparrow \uparrow +\downarrow \downarrow )/\sqrt 2$	0.7	0.4
$\Psi _{\mathrm {3}}= (\uparrow \uparrow -\downarrow \downarrow )/\sqrt 2$	0.7	1.1
$\Psi _{\mathrm {1}}= (\uparrow \downarrow +\downarrow \uparrow )/\sqrt 2$	1.0	1.1
(J,K,D)	(1.0,−0.6,0.0)	(1.1,−0.7,−0.7)
$\measuredangle$(Ir–O–Ir)=98.5° :
$\Psi _{\mathrm {2}}= (\uparrow \uparrow +\downarrow \downarrow )/\sqrt 2$	0.0	0.0
$\Psi _{\mathrm {3}}= (\uparrow \uparrow -\downarrow \downarrow )/\sqrt 2$	1.2	2.1
$\Psi _{\mathrm {S}}= (\uparrow \downarrow -\downarrow \uparrow )/\sqrt 2$	4.3	5.1
$\Psi _{\mathrm {1}}= (\uparrow \downarrow +\downarrow \uparrow )/\sqrt 2$	3.9	6.5
(J,K,D)	(−0.4,−6.6,−1.2)	(1.4,−10.9,−2.1)

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, Ψ₁, Ψ₂ and Ψ₃ 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₁,M_L1,L₂,M_L2,S,M_S} basis in standard QC programs to the $\{{\tilde {S}}_1,{\tilde {S}}_2,{\tilde M}_{S_1},{\tilde M}_{S_2}\}$ basis. This allows the study of Ψ₁–Ψ₂ mixing when the point-group symmetry is reduced to C_2h, see below. Having the assignment of the states resolved, the Ψ_S–Ψ₁ splitting provides J, the Ψ₂–Ψ₃ splitting yields D, while the difference between the energy of Ψ₁ and the average of the E₂(Ψ₂) and E₃(Ψ₃) 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₂O₂ plaquette is not a negligible factor, as presently assumed in simplified superexchange models for Kitaev physics in honeycomb Na₂IrO₃ [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 [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₂IrO₃, displaying different Ir–O–Ir bond angles and slightly different Ir–O and Ir–Ir distances [8]. While the [Ir₂O₁₀] 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 $\{\tilde {S}_1,\tilde {S}_2,\tilde {M}_{S_1},\tilde {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).

Table 2. Energy splittings of the four lowest magnetic states and effective coupling parameters (meV) for two NN IrO₆ octahedra in the C2/m structure of [8]. The weight of $(\uparrow \downarrow + \downarrow \uparrow )/\sqrt {2}$ and $(\uparrow \uparrow + \downarrow \downarrow )/\sqrt {2}$ in Ψ'₁ and Ψ'₂, respectively, is ≈98%, see text.

Method	CAS+SOC	MRCI+SOC
$\measuredangle$(Ir–O–Ir)=99.45°, d (Ir1–Ir2)=3.138 Å(×1)a:
Ψ'2	0.0	0.0
$\Psi _{\mathrm {3}}= (\uparrow \uparrow - \downarrow \downarrow )/\sqrt 2$	0.2	0.5
$\Psi _{\mathrm {S}}= (\uparrow \downarrow - \downarrow \uparrow )/\sqrt 2$	4.4	5.5
Ψ'1	6.3	10.5
(J,K,D)	(1.9,−12.4,−0.2)	(5.0,−20.5,−0.5)
$\measuredangle$(Ir–O–Ir)=97.97°, d (Ir2–Ir3)=3.130 Å(×2)b:
Ψ'2	0.0	0.0
$\Psi _{\mathrm {3}}= (\uparrow \uparrow - \downarrow \downarrow )/\sqrt 2$	0.3	1.2
$\Psi _{\mathrm {S}}= (\uparrow \downarrow - \downarrow \uparrow )/\sqrt 2$	4.6	6.7
Ψ'1	5.8	8.2
(J,K,D)	(1.2,−11.3,−0.3)	(1.5,−15.2,−1.2)

^ad (Ir–O_1,2) = 2.056 Å. ^bd (Ir–O₁) = 2.065 Å, d (Ir–O₂) = 2.083 Å.

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 [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₆ 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₂IrO₃ 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, 16, 31, 35]. We therefore performed ED calculations for a KH model supplemented with second and third NN couplings J₂ and J₃ (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₂ and J₃ 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 figure 2, shows that the zigzag phase is indeed stable in the region of J₂,J₃ ≳ 2 meV. We note that positive J₂ and J₃ values of 4–5 meV would be consistent with the experimentally observed Curie–Weiss temperature ≈ − 125 K [35] using $\theta _{\rm{CW}}=-\tilde {S}(\tilde {S}+1)(J+2J_2+J_3+K/3)/k_{\mathrm { B}}$ [8]. Thus we propose that an extended spin Hamiltonian based on the NN 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₂IrO₃.

Zoom In Zoom Out Reset image size

To discuss in more detail the generic phase diagram in figure 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₁–J₂–J₃ model [36–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₂ and J₃. These findings are relevant in the context of recent experimental data that indicate a qualitatively different AF GS for the related compound Li₂IrO₃ [39].

In summary, for the honeycomb iridate Na₂IrO₃, the ab initio QC calculations show that in a reference system with X along the Ir–Ir link and Z perpendicular on the Ir₂O₂ 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–20 meV, while the NN Heisenberg superexchange and the off-diagonal xy coupling are significantly weaker. The QC data additionally 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₂IrO₃, Li₂IrO₃ and Li₂RhO₃ 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 GSs in different 5d⁵ or 4d⁵ 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–42] but clearly 5d systems with stronger spin–orbit couplings may provide new playgrounds in this research area too.

Note added in proof. After submission of our paper, results of ED calculations including off-diagonal D type terms were also reported in [43].

We thank N A Bogdanov, G Khaliullin, D I Khomskii, H Gretarsson, Y-J Kim, H Stoll and P Fulde for insightful discussions. LH acknowledges financial support from the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG).

For the computation of the intersite spin couplings, two reference NN IrO₆ 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⁴⁺ 5d⁵ NN's with closed-shell Pt⁴⁺ 5d⁶ species. This is an usual procedure in QC studies on transition-metal systems, see, e.g.  [20–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 QC 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⁴⁺ 5d⁶ 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 nine-root 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 Ψ₁, Ψ₂ or Ψ₃, respectively, defined as in table 1, should not change energy. This provides a quick alternative way of identifying the low-lying spin–orbit states.

In the main body of this article, we made use of two different local reference frames, namely {x_b,y_b,z_b} and { $\mathbf {X}_b,\mathbf {Y}_b,\mathbf {Z}_b\} =\{\frac {\mathbf {x}_b+\mathbf {y}_b}{\sqrt {2}},\frac {-\mathbf {x}_b+\mathbf {y}_b}{\sqrt {2}},\mathbf {z}_b\}$, for each of the three bond types b = 1 − 3, see figure B.1. Choosing the global frame to be {x,y,z} ≡ {x₁,y₁,z₁}, the local frames {X_b,Y_b,Z_b} are expressed as

$$\begin{eqnarray} \{\mathbf{X}_1,\mathbf{Y}_1,\mathbf{Z}_1\}&=&\left\{\frac{\mathbf{x}+\mathbf{y}}{\sqrt{2}}, \frac{-\mathbf{x}+\mathbf{y}}{\sqrt{2}}, \mathbf{z}\right\},\nonumber\\ \{\mathbf{X}_2,\mathbf{Y}_2,\mathbf{Z}_2\}&=&\left\{\frac{\mathbf{x}-\mathbf{z}}{\sqrt{2}}, -\frac{\mathbf{x}+\mathbf{z}}{\sqrt{2}}, \mathbf{y}\right\},\\ \{\mathbf{X}_3,\mathbf{Y}_3,\mathbf{Z}_3\}&=&\left\{ \frac{\mathbf{y}+\mathbf{z}}{\sqrt{2}}, \frac{-\mathbf{y}+\mathbf{z}}{\sqrt{2}}, \mathbf{x}\right\}.\nonumber \end{eqnarray} \tag{ B.1 }$$

We note that the system is invariant under a two-fold rotation around X₁ and that the directions of X₂ and X₃ are chosen to map to each other under this rotation (otherwise the sign of the exchange coupling parameter C is arbitrary).

Zoom In Zoom Out Reset image size

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 {X_b,Y_b,Z_b}:

$$\begin{eqnarray} &&\fl \mathcal{H}_{\langle ij\rangle \in b} = J_{0,b} \boldsymbol{\tilde{S}}_i\cdot\boldsymbol{\tilde{S}}_j + A_b {\tilde{S}}_i^{X_b}\tilde{S}_j^{X_b} + B_b \tilde{S}_i^{Y_b}\tilde{S}_j^{Y_b} -(A_b+B_b) \tilde{S}_i^{Z_b}\tilde{S}_j^{Z_b} + C_b (\tilde{S}_i^{Y_b}\tilde{S}_j^{Z_b}+\tilde{S}_i^{Z_b}\tilde{S}_j^{Y_b})\nonumber\\ &&=J_b\boldsymbol{\tilde{S}}_i\cdot\boldsymbol{\tilde{S}}_j + K_b {\tilde{S}}_i^{Z_b}{\tilde{S}}_j^{Z_b} + D_b ({\tilde{S}}_i^{X_b}{\tilde{S}}_j^{X_b}-{\tilde{S}}_i^{Y_b}{\tilde{S}}_j^{Y_b}) + C_b ({\tilde{S}}_i^{Y_b}{\tilde{S}}_j^{Z_b}+{\tilde{S}}_i^{Z_b}{\tilde{S}}_j^{Y_b}), \end{eqnarray} \tag{ B.2 }$$

where $J_b=J_{0,b}+\frac {1}{2}(A_b+B_b)$, $K_b=-\frac {3}{2}(A_b+B_b)$ and $D_b=\frac {1}{2}(A_b-B_b)$. We can rewrite these terms in the global frame {x,y,z} using (B.1)

$$\begin{eqnarray} &&\fl \mathcal{H}_{\langle ij\rangle \in b=1} = J_1 \boldsymbol{\tilde{S}}_i\cdot\boldsymbol{\tilde{S}}_j + K_1 {\tilde{S}}_i^{z}{\tilde{S}}_j^{z} + D_1 ( {\tilde{S}}_i^{x}{\tilde{S}}_j^{y}+{\tilde{S}}_i^{y}{\tilde{S}}_j^{x})+ \frac{C_1}{\sqrt{2}} ({\tilde{S}}_i^{y}{\tilde{S}}_j^{z}+{\tilde{S}}_i^{z}{\tilde{S}}_j^{y}-{\tilde{S}}_i^{x}{\tilde{S}}_j^{z}-{\tilde{S}}_i^{z}{\tilde{S}}_j^{x})\,, \nonumber\\ &&\fl \mathcal{H}_{\langle ij\rangle \in b=2} = J_2 \boldsymbol{\tilde{S}}_i\cdot\boldsymbol{\tilde{S}}_j + K_2 {\tilde{S}}_i^{y}{\tilde{S}}_j^{y} + D_2(-{\tilde{S}}_i^{x}{\tilde{S}}_j^{z}-{\tilde{S}}_i^{z}{\tilde{S}}_j^{x})+\frac{C_2}{\sqrt{2}} (-{\tilde{S}}_i^{x}{\tilde{S}}_j^{y}-{\tilde{S}}_i^{y}{\tilde{S}}_j^{x}-{\tilde{S}}_i^{z}{\tilde{S}}_j^{y}-{\tilde{S}}_i^{y}{\tilde{S}}_j^{z}),\nonumber\\ &&\fl \mathcal{H}_{\langle ij\rangle \in b=3} = J_3 \boldsymbol{\tilde{S}}_i\cdot\boldsymbol{\tilde{S}}_j + K_3 {\tilde{S}}_i^{x}{\tilde{S}}_j^{x} + D_3 ({\tilde{S}}_i^{y}{\tilde{S}}_j^{z}+{\tilde{S}}_i^{z}{\tilde{S}}_j^{y})+{ \frac{C_3}{\sqrt{2}} (-{\tilde{S}}_i^{x}{\tilde{S}}_j^{y}-{\tilde{S}}_i^{y}{\tilde{S}}_j^{x}+{\tilde{S}}_i^{x}{\tilde{S}}_j^{z}+{\tilde{S}}_i^{z}{\tilde{S}}_j^{x})}. \end{eqnarray} \tag{ B.3 }$$

For Na₂IrO₃ bonds 2 and 3 are equivalent and therefore J₂ = J₃, K₂ = K₃, D₂ = D₃, C₂ = C₃.

In C_2h symmetry, D_b = 0 and for a given bond b the eigenvalues of the Hamiltonian defined by equation (B.2) are

$$\begin{equation} E_{\mathrm{S}}=-\frac{3J_{0,b}}{4}, \end{equation} \tag{ B.4a }$$

$$\begin{equation} E_1^{\prime} =\frac{J_{0,b}+A_b+\sqrt{(A_b+2B_b)^2+4C_b^2}}{4}, \end{equation} \tag{ B.4b }$$

$$\begin{equation} E_2^{\prime} =\frac{J_{0,b}+A_b-\sqrt{(A_b+2B_b)^2+4C_b^2}}{4}, \end{equation} \tag{ B.4c }$$

$$\begin{equation} E_3 =\frac{J_{0,b}-2A_b}{4}\,. \end{equation} \tag{ B.4d }$$

These expressions become more complicated for point-group symmetries lower than C_2h. The analysis of the spin–orbit wave functions in the transformed $\{{\tilde {S}}_1,{\tilde {S}}_2,{\tilde M}_{S_1},{\tilde M}_{S_2}\}$ basis shows, however, that the mixing of the D_2h Ψ_i eigenvectors for distorted clusters and lower symmetries of the [Ir₂O₁₀] blocks (see table 2) is negligible in the C2/m structure determined by Choi et al [8].

For D_2h symmetry (see table 1), C_b = 0 and the eigenvalues in the effective two-site, four-state problem are E_S = −3J_0,b/4, E₁ = (J_0,b + 2A_b + 2B_b)/4, E₂ = (J_0,b − 2B_b)/4 and E₃ = (J_0,b − 2A_b)/4.