Spin–orbital locking, emergent pseudo-spin and magnetic order in honeycomb lattice iridates

The nature of the effective spin Hamiltonian and magnetic order in the honeycomb iridates is explored by considering a trigonal crystal field effect and spin–orbit (SO) coupling. Starting from a Hubbard model, an effective spin Hamiltonian is derived in terms of an emergent pseudo-spin-1/2 moment in the limit of large trigonal distortions and SO coupling. The present pseudo-spins arise from a spin–orbital locking and are different from the jeff = 1/2 moments that are obtained when the SO coupling dominates and trigonal distortions are neglected. The resulting spin Hamiltonian is anisotropic and frustrated by further neighbour interactions. Mean-field theory suggests a ground state with four-sublattice zigzag magnetic order in a parameter regime that can be relevant to the honeycomb iridate compound Na2IrO3, where a similar magnetic ground state has recently been observed. Various properties of the phase, the spin-wave spectrum and experimental consequences are discussed. The present approach contrasts with the recent proposals to understand iridate compounds starting from the strong SO coupling limit and neglecting non-cubic lattice distortions.


Introduction
The interplay between strong spin-orbit (SO) coupling and electron-electron interaction in correlated electron systems has recently been a subject of intensive study . In particular, 5d transition metal (e.g. iridium (Ir) or osmium (Os)) oxides are regarded as ideal playgrounds for observing such cooperative effects . Compared to 3d transition metal oxides, the repulsive Coulomb energy scale in these systems is reduced by the much larger extent of 5d orbitals, while the SO coupling is enhanced due to the high atomic number (Z = 77 for Ir and z = 76 for Os). Moreover, owing to the extended 5d orbitals, these systems are very sensitive to the crystal fields. As a result, the energy scales mentioned above often become comparable to each other, leading to a variety of competing phases. Precisely for this reason, one expects to see newer emergent quantum phases in such systems. Indeed, there have been several theoretical proposals, in the context of concrete experimental examples, for spin liquids [1,11,12,14,16,18,19], topological insulators [10,19,20], Weyl semimetals [21,22], novel magnetically ordered Mott insulators [2-6, 23, 24] and other related phases [8,9,23] in iridium and osmium oxides.
A typical situation in the iridates consists of Ir 4+ atoms sitting in the octahedral crystal field of a chalcogen, typically oxygen or sulfur [1][2][3][4]7]. This octahedral crystal field splits the five 5d orbitals of Ir into the doubly degenerate e g orbitals and the triply degenerate t 2g orbitals (each orbital has a further twofold spin degeneracy). The e g orbitals are higher in energy with the energy difference being approximately 3 eV. There are five electrons in the outermost 5d shell of Ir 4+ which occupy the low-lying t 2g orbitals, and the low-energy physics is effectively described by projecting out the empty e g orbitals [25]. A characteristic feature of most of the approaches used to understand these compounds is to treat the SO coupling as the strongest 3 b a Figure 1. The zigzag magnetic structure as found in [5]. The magnetic unit cell has four sites.
interaction at the atomic level, i.e. by considering the effect of extremely strong SO coupling for electrons occupying the t 2g orbitals. This decides the nature of the participating atomic orbitals in the low-energy effective theory. In this limit, the orbital angular momentum, projected to the t 2g manifold, carries an effective orbital angular momentum l eff = 1 with a negative SO coupling constant [2,3,15]. The projected SO coupling splits the t 2g manifold into the lower j eff = 3/2 quadruplet and the upper j eff = 1/2 doublet. Out of the five valence electrons, four fill up the quadruplet sector, leaving the doublet sector half-filled. Thus, in the limit of very strong SO coupling, the half-filled doublet sector emerges as the correct low-energy degrees of freedom. Considering the effect of Coulomb repulsion within a Hubbard model description and performing strong coupling expansion, various spin Hamiltonians for j eff = 1/2 are then derived within a strong-coupling perturbation theory [11,20,23].
In this paper, we, however, consider a different limit where the oxygen octahedra surrounding the Ir 4+ ions are highly distorted. While the above scenario of half-filled j eff = 1/2 orbitals is applicable to the undistorted case, as we shall see, it breaks down in the presence of strong distortions of the octahedra. In particular, we consider the effect of trigonal distortions, which may be relevant for some of the iridate systems including the much debated honeycomb lattice iridate, Na 2 IrO 3 . We show that, in this limit, a different 'doublet' of orbitals emerges as the low-energy degree of freedom. This doublet forms a pseudo-spin-1/2 that results from a kind of (physical) spin-orbital locking, so that the spin and orbital fluctuations are not separable (as discussed below). We emphasize that this pseudo-spin is different from the j eff = 1/2 doublet discussed above. The spin Hamiltonian for these pseudo-spins (equation (6)), on a honeycomb lattice, admits a four-sublattice zigzag (figure 1) pattern in a relevant parameter regime. Such a magnetic order has recently been observed in the experiments [5,28] on Na 2 IrO 3 and hence our theory may be applicable to this material.
Distortion of the octahedron surrounding the Ir 4+ generates a new energy scale associated with the change in the crystal field, which, as we shall see, competes with the SO coupling. Several kinds of distortion may occur, of which we consider the trigonal distortions of the octahedron where it is stretched/compressed along the body diagonal of the enclosing cube [25]. In the absence SO coupling, such a trigonal crystal field splits the t 2g manifold into e g (with two degenerate orbitals e 1g and e 2g ) and non-degenerate a 1g (again there is an added twofold spin degeneracy for each of these orbitals). The e g and a 1g levels are, respectively, occupied by three and two electrons in Ir 4+ . For large trigonal distortions, the splitting between them is big and the a 1g orbitals can be projected out. Now, if one adds SO coupling, the low-energy degrees of freedom are described by a subspace of the e g orbitals which form an emergent pseudospin-1/2 doublet out of |e 1g , ↓ and |e 2g ↑ states, where ↑ and ↓ represent the physical spin s z = 1/2, −1/2 (the spins are quantized along the axis of trigonal distortion). These pseudospin-1/2 are different from the j eff = 1/2 and j eff = 3/2 multiplets in the strong SO coupling limit as discussed above. Note the (physical) spin-orbital locking for the pseudo-spins, as alluded to above.
The two approaches, described in the last two paragraphs, for arriving at the low-energy manifold are mutually incompatible. This can be seen as follows: in the presence of sizeable trigonal distortions the j eff = 1/2 and j eff = 3/2 multiplets mix with each other and can no longer serve as good low-energy atomic orbitals. This dichotomy becomes quite evident in recent studies of Na 2 IrO 3 , where the Ir 4+ form a honeycomb lattice. Taking into account the strong SO coupling of Ir 4+ in Na 2 IrO 3 , a model for a topological insulator in the weak interaction limit [10] and a Heisenberg-Kitaev (HK) model for a possible spin liquid phase in the strong coupling limit [11] are proposed. These proposals prompted several experimental [4][5][6] and theoretical efforts [12,26] to understand the nature of the ground state in this material. Subsequently, it was found that Na 2 IrO 3 orders magnetically at low temperatures [4]. However, the magnetic moments form a 'zigzag' pattern (figure 1) which is not consistent with the ones that would be obtained by adding a weak interaction in a topological insulator (canted antiferromagnet) [10] or from a nearest-neighbour (NN) HK model (spin liquid or the socalled stripe antiferromagnet) [11]. While recent studies [27] show that a 'zigzag' order may be stabilized within the HK model by including substantial second and third neighbour antiferromagnetic interactions, it is hard to justify such large further neighbour exchanges without lattice distortions. (An alternative explanation that we do not pursue here is significant charge fluctuations, which would mean that the compound is close to metal-insulator transition. The resistivity data seem to support the fact that this compound is a good insulator [4].) If lattice distortions are responsible for the significant further neighbour exchanges, then there would be sizeable distortion of the oxygen octahedra, which in turn may invalidate the above j eff = 1/2 picture and thus the basic paradigm of the HK model, by mixing the j eff =1/2 and j eff = 3/2 subspaces. Recent finite-temperature numerical calculations on the HK model [26] also suggest possible inconsistencies with experiments on Na 2 IrO 3 . This necessitates the need for a different starting point to explain the magnetic properties of Na 2 IrO 3 .
The rest of this paper is organized as follows. In section 2, we derive the effective spin Hamiltonian in the limit of the large trigonal distortion and large SO coupling. This is done by taking the energy scale associated with trigonal distortion to infinity first, followed by that of the SO energy scale. This order of taking the limit gives a spin Hamiltonian in terms of emergent pseudo-spin-1/2, which is different from the HK model. This Hamiltonian is both anisotropic and frustrated. It also has further neighbour interactions, the effects of which are enhanced due to anisotropy that makes some of the NN bonds weaker. The origin of this anisotropy is trigonal distortion. We argue that this limit may be more applicable for the compound Na 2 IrO 3 . Having derived the spin Hamiltonian, we calculate the phase diagram and the 5 spin-wave spectrum within mean-field theory in section 3. We see that the 'zigzag' phase occurs in a relevant parameter regime. We also point out the experimental implications of our calculations in the context of Na 2 IrO 3 . Finally, we summarize the results in section 4. The details of various calculations are given in various appendices.

The effective Hamiltonian
In the cubic environment the t 2g orbitals are degenerate when there is no SO coupling. Trigonal distortion due to compression or expansion along one of the four C 3 axes of IrO 6 octahedra lifts this degeneracy. Although it is possible that the axes of trigonal distortions are different in different octahedra [20], we find that uniform distortions are more consistent with the experiments (see below) on Na 2 IrO 3 . Hence we consider uniform trigonal distortion.

The trigonal Hamiltonian
Let us denote the axis of this uniform trigonal distortion of the octahedron by the unit vectorn = 1 √ 3 [n 1 , n 2 , n 3 ], where n α = ±1. Since there are two directions to each of the four trigonal axes we may choose a 'gauge' to specifyn. This is done by taking n 1 n 2 n 3 = +1. The Hamiltonian for trigonal distortion, when projected in the t 2g sector, gives [20] (in our chosen gauge) x y ] and tri is the energy scale for trigonal distortion. The eigenstates are (ω = e i2π/3 ) The trigonal distortion splits the t 2g sector into the doubly degenerate e g and the non-degenerate a 1g with energies tri /3 and −2 tri /3, respectively. A description based on the Hubbard model for the e g orbitals may be systematically derived starting from the t 2g orbitals. This is done in appendix A. This has the following general form: where e † i Mσ is the electron creation operator in the e g orbital (M = 1, 2) with spin σ (=↑, ↓); t i M; j M are the effective hopping amplitudes within the subspace and U is the effective on-site Coulomb's repulsion. We note that Hund's coupling (which arises from the orbital dependence of the Coulomb repulsion) for the t 2g orbitals only renormalizes U in this restricted subspace. (See appendix A for details.)

The projected spin-orbit (SO) coupling
The SO coupling, when projected in the e g subspace, yields a block diagonal form (see appendix B for details): where s i is the spin operator at the site i, λ ≈ 500 meV is the SO coupling parameter and τ z = +1(−1) refers to the e 1g (e 2g ) orbital. Thus the projected SO interaction acts as a 'Zeeman coupling' where the direction of the 'magnetic field' is along the trigonal axis or opposite to it [13]. Thus it is natural to choose the direction of spin quantization along the axis of trigonal distortion. This then gives the active atomic orbitals after incorporating the SO coupling. These active orbitals are Krammer's doublet |e 1g , ↓ and |e 2g , ↑ as shown in figure 2.

The spin Hamiltonian
Hence the low-energy physics may be described by considering only the above atomic orbitals. The starting point for the calculations is projection of the Hubbard model (equation (3)) in the space spanned by Krammers's doublet |e 1g , ↓ and |e 2g , ↑ . The bandwidth of this projected model is narrow and the effect of the Hubbard repulsion is important. Indeed, it can easily render the system insulating. To capture the magnetic order in this Mott insulator, we do a strong-coupling expansion int/U to obtain an effective 'pseudo-spin' model in terms of the pseudo-spin-1/2 operators, where, ρ α (α = x, y, z) are the Pauli matrices and a, b = (e g1 ; ↓), (e g2 ; ↑). The 'pseudo-spin' Hamiltonian has the following form up to the quadratic order: Here i j , i j and i j refer to summation over the first, second and third NNs, respectively.

7
The different exchange couplings are given in terms of the underlying parameters of the Hubbard model as where α = 1, 2, 3 denotes that i j are the first, second or third NNs, respectively and the last expression defines J (0α) i j . T (0α) Before moving on to the details of the spin Hamiltonian, we note that, on projecting to the subspace of |e 1g , ↓ and |e 2g , ↑ , the spin and orbitals are no longer independent. Instead, at every site there is a pseudo-spin-1/2 degree of freedom where the spin is locked to the orbital wave function. This we refer to as spin-orbital locking.

Application to Na 2 IrO 3
We now apply the above results to the case of Na 2 IrO 3 . The early x-ray diffraction experiments [4] suggested a monoclinic C2/C structure for the compound and distorted IrO 6 octahedra. However, more recent experiments see a better match for the x-ray diffraction data with the space group C2/m [28,29]. They also unambiguously confirm the presence of uniform trigonal distortion of the IrO 6 octahedra. However, the magnitude of such a distortion is not clear at present. Further, experimental measurements suggest: (i) the magnetic transition occurs at T N = 15 K, while the Curie-Weiss temperature is about CW ≈ −116 K. This indicates the presence of frustration. (ii) The high-temperature magnetic susceptibility is anisotropic; the inplane and out-of-plane susceptibilities are different. This may be due to a trigonal distortion of the IrO 6 octahedra [4]. (ii) The magnetic specific heat is suppressed at low temperatures [4]. (iv) A recent resonant x-ray scattering experiment [5] suggests that the magnetic order is collinear and has a four-site unit cell. (v) The magnetic moments have a large projection on the a-axis of the monoclinic crystal [5]. (vi) A combination of these experimental findings and density functional theory (DFT) calculations strongly suggest a 'zigzag' pattern for the magnetic moments, as shown in figure 1 in the ground state [5], which has since been verified independently by two groups using neutron scattering [28,29].
Taking these phenomenological suggestions, we try to apply the above calculations to the case of Na 2 IrO 3 . At the outset, we must note that, in the above derivation of the spin Hamiltonian, we have assumed the trigonal distortion to be the largest energy scale, followed by the SO coupling. While this extreme limit of projecting out the a 1g orbitals is most likely not true for Na 2 IrO 3 , we expect the real ground state to be adiabatically connected to this limit. With this in mind, we now consider the case of Na 2 IrO 3 . Clearly, the exchanges (equation (7)) depend on both the direction of the bond and the direction of the trigonal distortion. So it is important to ask about the direction of the latter. Comparing the crystallographic axes of Na 2 IrO 3 , we find that the direction [1, 1, 1] is perpendicular to the honeycomb plane, while the other three directions make an acute angle to it. In the monoclinic C2/m structure, uniform trigonal distortion in these four directions may not cost the same energy. In experiments [5], the moments are seen to point along the a-axis of the monoclinic crystal which is parallel to the honeycomb plane. This, along with the fact that the magnetic moment in our model is in the direction ofn (explained below), seems to suggest thatn = 1 √ 3 [−1, −1, 1] is chosen in the compound (see figure 3). In the absence of a better theoretical understanding of the direction of the trigonal distortion, we take this as an input from the experiments.
To identity different hopping paths (both direct and indirect), we consider various overlaps (see appendix C) and find that, while J (3z) = 0, J (1z) = J (2z) = 0 are approximately (spatially) isotropic and antiferromagnetic. For the exchanges of the Heisenberg terms, both J (2) and J (3) are antiferromagnetic and isotropic (both of them result from indirect hopping mediated by the Na s-orbitals and are expected to be comparable). For the NN Heisenberg exchanges, the couplings are antiferromagnetic, but much more spatially anisotropic. We find that for the chosen direction of the trigonal distortion, the coupling along one of the NN exchanges (J (1) ) (viz b 1 in figure 3) is different from the other two neighbours (J (1) ) (b 2 and b 3 in figure 3).

Mean-field theory and magnetic order
We now consider the mean-field phase diagram for the above anisotropic spin Hamiltonian. For J (1) being the largest energy scale, the classical ground state for the model can be calculated within mean-field theory as a function of x 0 =J (1) /J (1) and y 0 = J (2) /J (1) (we have taken J (2) = J (3) ). A representative mean-field phase diagram is shown in figure 4. It shows a J (1) ; y 0 = J (2) J (1) , where J (1) (J (1) ) are related to the strong(weak) NN exchange and J (2) is the second and third neighbour exchanges (see equation (7)). We take the Ising anisotropy to be 5% of J (1) . Note that, due to Ising anisotropies, one has zigzag order at y 0 = 0. region of the parameter space where the zigzag order is stabilized 5 . The effect of the Ising anisotropies J (1z) and J (2z) is to pin the magnetic ordering along the z-direction of the pseudospin quantization which is also the direction of the trigonal distortionn. They also gap out any Goldstone mode that arises from the ordering of the pseudo-spins. The latter results in the exponential suppression of the specific heat at low temperatures. The other competing phase with a collinear order is the regular two-sublattice Neel phase.
The nature of the ground states may be understood from the following arguments. In the presence ofn in the [−1, −1, 1] direction, the NN exchange coupling becomes anisotropic. When it is strong in one direction (J (1) ) and weak in two other directions (J (1) ), for the bonds where the NN coupling becomes weak, the effects of the small second and third neighbour interactions become significant. Since the latter interactions are antiferromagnetic, they prefer anti-parallel alignment of the spins. As there are more second and third neighbours, their cumulative effect can be much stronger. This naturally leads to the zigzag state. The NN antiferromagnetic interactions on the weaker bonds compete with the antiferromagnetic second and third neighbour interactions and frustrate the magnet. This suppresses the magnetic ordering temperature far below the Curie-Weiss temperature.

The spectrum for spin-orbital waves
The low-energy excitations about this magnetically ordered zigzag state are gapped spin-orbital waves. Signatures of such excitations may be seen in future resonant x-ray scattering experiments. It is important to note that these 'pseudo-spin' waves actually contain both orbital and the spin components due to spin-orbital locking.
We calculate the dispersion of such spin-orbital waves to quadratic order using the wellknown Holstein-Primakoff methods. The details are given in appendix D. A representative  panel (b)). The values used for the parameters are the same as that used for the calculation of the mean-field phase diagram ( figure 4). We note that, as expected, the spectrum is gapped.
spin-wave spectrum in the zigzag phase is shown in figures 5(a) and (b). The spectrum is gapped and the bottom of the spin-wave dispersion has some characteristic momentum dependence.

Experimental implications
Apart from the already discussed exponential suppression of low-temperature magnetic specific heat, the above calculation predicts an interesting feature in the magnetic susceptibility. The relation between the magnetic moment and the pseudo-spins is where µ B is the Bohr magneton. This follows from the twin facts that, in e g subspace, the angular momentum transverse ton is quenched and the spins are locked to the orbitals with the axis of quantization beingn in our pseudo-spin sector (see appendix E). Thus, the magnetization is sensitive to the z-component of the pseudo-spin (the direction of which is shown in figure 1). Indeed the magnetization has the largest projection along the a-axis of the monoclinic crystal. This was seen in experiments [5] and was the motivation for choosing the [− . Equation (10) suggests that the magnetic susceptibility is highly anisotropic and depends on the cosine of the angle between the direction of the magnetic field andn. Indeed, signatures of such anisotropy have already been seen in experiments [4]. We emphasize that within this picture, the in-plane susceptibility also varies with the direction of the magnetic field. So the ratio χ ⊥ /χ can be less than or greater than 1.
The current experiments [4] do not tell the in-plane direction of the magnetic field, and hence, we cannot comment on the ratio at present. However, the above picture is strictly based on atomic orbitals. One generally expects that there is also hybridization of the Ir d-orbitals with the oxygen p-orbitals. Such a hybridization will contribute a non-zero isotropic component to the susceptibility [24]. Also, as remarked earlier, in the actual compound, the SO coupling scaling may not be very small compared to the trigonal distortion limit scale. An additional perturbation coming from the mixing with the a 1g orbitals will also contribute to decreasing the anisotropy of the susceptibility.

Summary and conclusion
In this paper, we have studied the effect of trigonal distortion and SO coupling and applied it to the case of the honeycomb lattice compound Na 2 IrO 3 . We find that, in the limit of large trigonal distortion and SO coupling, a pseudo-spin-1/2 degree of freedom emerges. The lowenergy Hamiltonian, in terms of this pseudo-spin, gives a 'zigzag' magnetic order as seen in recent experiments on Na 2 IrO 3 . We have also calculated the low-energy spin-wave spectrum and elucidated various properties of the compound that have been observed in experiments. The pseudo-spin couples the physical spin and the orbitals in a non-trivial manner, signatures of this may be seen in future inelastic x-ray resonance experiments probing the low-energy excitations.
While very recent experiments [28,29] clearly indicate the presence of trigonal distortions, their magnitude is as yet not confirmed. On the other hand, the only available numerical estimate of the energy scale for trigonal distortion comes from the DFT calculations by Jin et al [13] (based on C2/C structure). It suggests that tri ≈ 600 meV. While it is not clear whether such a large value is in conformity with the experiments, at present the detection of trigonal distortion in experiments is highly encouraging from the perspective of the present calculations.
In the light of the above calculations, it is tempting to predict the case of Li 2 IrO 3 where recent experiments suggest a more isotropic honeycomb lattice [6,30]. A possibility is that sizeable trigonal distortion is also present in Li 2 IrO 3 (so that the above discussion holds), but the axis is perpendicular to the plane. What may be the fallouts in such a case? Our present analysis would then suggest that the antiferromagnetic exchanges are isotropic and equally strong for the three NNs. This would develop two-sublattice Neel order in the pseudo-spins with the magnetic moments being perpendicular to the plane. Also the further neighbour exchanges are rather weak (compared to Na 2 IrO 3 ) and hence frustration is quite small. Indeed recent experiments see ordering very close to the Curie-Weiss temperature, the later being calculated from the hightemperature magnetic susceptibility data [6,30]. However, the present experiments do not rule out the possibility of small or no trigonal distortions in Li 2 IrO 3 , in which case the limit of the HK model [11,26] may be appropriate.
show that there are contributions from both direct and indirect hoppings for the first, second and third NNs. Projecting them into the e g orbitals we obtain the effective hopping amplitudes, which are then used in equation (3). As to the Coulomb repulsion term, we find that it has the following form:Ũ where U = U 0 − 2J H /3. This form is then used in equation (3). The reason for this special form ofŨ M 1 M 2 lies in the fact that the e g orbitals have equal weight of the three t 2g orbitals (see the wave functions in equation (2)).  We shall make an approximation here. We shall leave out the directional dependence of the magnitudes on the direction. The argument is that the essential directional dependence due to the trigonal distortion has been taken care of by the parameter 1 . When the DFT [13] resultsareusedtofindthetight-bindingparameters(seefootnote5),itisfoundthat(theyuse

C.3. Third NN
The third NNs are listed in figure C.1(c). The hopping to the third NN is mediated by the Na atoms. Again we shall neglect the directional dependence and take these to be in the magnitudes of the hopping amplitudes. The result is summarized in table C.3.
Tight-bindingfittotheDFTresults( [13];seefootnote5)indeedshowsthatthishopping energy scale is of the order of for the other direction. Since there are four sites per unit cell (see figure 1(a) of the main text) the quadratic Hamiltonian is a 8 × 8 matrix given by where H cl is the classical part dealt with in the previous section. The spin-wave Hamiltonian has the following form: Here † k = [a † k,1 , a † k,2 , a † k,3 , a † k,4 , a −k,1 , a −k,2 , a −k,3 , a −k,4 ] (the subscript 1, 2, 3, 4 refers to the four sites in the unit cell as shown in figure 1(a) of the main text) and where N cell is the number of unit cells and