Hidden symmetries of the extended Kitaev-Heisenberg model: Implications for honeycomb lattice iridates {\boldsymbol{A}}_{\mathbf{2}}IrO{}_{\mathbf{3}}

Hidden symmetries of the extended Kitaev-Heisenberg model: Implications for honeycomb lattice iridates IrO

Jiří Chaloupka Central European Institute of Technology, Masaryk University, Kotlářská 2, 61137 Brno, Czech Republic    Giniyat Khaliullin Max Planck Institute for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany
July 12, 2019
Abstract

We have explored the hidden symmetries of a generic four-parameter nearest-neighbor spin model, allowed in honeycomb lattice compounds under trigonal compression. Our method utilizes a systematic algorithm to identify all dual transformations of the model that map the Hamiltonian on itself, changing the parameters and providing exact links between different points in its parameter space. We have found the complete set of points of hidden symmetry at which seemingly highly anisotropic model can be mapped back on the Heisenberg model and inherits therefore its properties such as the presence of gapless Goldstone modes. The procedure used to search for the hidden symmetries is quite general and may be extended to other bond-anisotropic spin models and other lattices, such as the triangular, kagome, hyper-honeycomb, or harmonic-honeycomb lattices. We apply our findings to the honeycomb lattice iridates NaIrO and LiIrO, and illustrate how they help to identify plausible values of the model parameters that are compatible with the available experimental data.

pacs:
75.10.Jm, 75.25.Dk, 75.30.Et

I Introduction

When relativistic spin-orbit coupling dominates over the exchange and orbital-lattice interactions, the orbital moment of an ion remains unquenched and a total angular momentum is formed. This was known to happen in compounds of late transition metal ions such as of cobalt (see, e.g., Ref. 1); however, the “cleanest” examples of spin-orbit coupled magnets emerged more recently: these are the iridium oxides SrIrO and NaIrO with perovskite and honeycomb lattice structures, correspondingly.

By construction, magnetic ordering in these systems necessarily involves interactions between orbital moments , in addition to a conventional Heisenberg exchange among the spin-part of total angular momentum [2]. Since the -moment, hosted by orbital in a crystal, is only an “effective” one [3], it need not be conserved during the electron hoppings, thus the -moment exchange interactions are generally not invariant [4]. Moreover, the orbital moments have a “shape” and hence the -interactions are anisotropic in real space, too, and thus strongly frustrated even on simple cubic lattices. Altogether, this results in nontrivial -Hamiltonians and orderings, including e.g. noncoplanar (multi-Q) states, “hidden” Goldstone modes, etc. [5; 6]. Via the spin-orbit coupling, these peculiar features of orbital physics are inherited by the “pseudospin-” wavefunctions and interactions [6; 7; 8; 9; 10; 11; 12; 13; 14]. In essence, frustrated nature and quantum behavior of -orbital moments [15; 16] are transferred to that of low-energy pseudospins .

Depending on the electron configuration of ions, the ground state pseudospin may take different values , and a variety of magnetic Hamiltonians with different symmetries and diverse behavior emerge in each case, because of different admixture of non-Heisenberg -interactions. Perhaps the most radical departure from a conventional magnetism is realized in compounds with apparently “nonmagnetic” ions, where a competition between spin-orbit and exchange interactions results in a nonmagnetic-magnetic quantum phase transition [17; 18; 19; 20].

The case of pseudospin iridates is of special interest. This is because SrIrO perovskite was found [22; 23; 24] to host cuprate-like magnetism, and honeycomb lattice iridates IrO (Na, Li) have been suggested [9] as a candidate material where the Kitaev model [25] physics might be realized. Following this proposal, a subsequent work [11] has introduced the minimal magnetic Hamiltonian for iridates IrO: the Kitaev-Heisenberg model (KH model) – a frustrated spin model with many attractive properties. Most importantly, its phase diagram contains a finite window of a quantum spin-liquid phase which emanates from the pure Kitaev point of the model with a known exact solution [25]. To reflect the later experimental findings in iridates, such as the zigzag (NaIrO [26; 27; 28]) and spiral (LiIrO [29]) type magnetic orderings, the initially proposed model was modified by including longer-range Heisenberg [30; 28] or anisotropic [31] interactions, extending the parameter range [32; 33; 34], by considering further anisotropic terms in the Hamiltonian [35; 36; 37; 38; 39; 40; 41; 42], or by including spatial anisotropy of the model parameters [43]. An alternative picture based on itinerant approach has been also suggested [44].

Despite the extensive efforts, no consensus concerning the minimal model for the honeycomb lattice iridates has thus far been reached. A reliable microscopic derivation of the exchange interactions is difficult and does not lead to a conclusive suggestion for the minimal Hamiltonian and its parameters. On the experimental side, the richest information about the underlying spin model would be provided by mapping momentum-resolved spin excitation spectrum. However, due to the lack of large enough monocrystals, the inelastic neutron scattering (INS) has been performed on powders only [28]. Another possible probe – resonant inelastic x-ray scattering (RIXS) – suffers from a small resolution at present. While it could be successfully applied in case of perovskite iridates [23], here the limitation comes from the much smaller energy scale of the excitations to be studied in detail by RIXS; however, the overall strength of magnetic interactions in NaIrO has been quantified [45; 46].

Nevertheless, the experimental data collected to date puts rather strong constraints on the possible models. First, the RIXS derived magnetic energies [45; 46] (of the order of 40 meV) are much higher than the ordering temperature () suggesting strong frustration. Second, the magnetic scattering intensity, measured by RIXS at zero momentum, , is as strong as elsewhere in the Brillouin zone, which implies a dominance of anisotropic, non-Heisenberg spin interactions. Third, the recent resonant x-ray scattering data [46] has revealed nearly ideal symmetry of the spin correlations in momentum space. Moreover, inelastic neutron scattering data [28] have indicated that a spin gap, if present, would be relatively small (less than ). All these observations taken together imply that the dominant pseudospin interactions in iridates are strongly frustrated, highly anisotropic in spin space, and yet highly symmetric in real space. By very construction, all these features are in fact the intrinsic properties of the KH model and its extended versions.

The KH model, supplemented by other symmetry allowed terms (see below), is therefore physically sound and plausible. However, there is a problem of its large parameter space (four parameters even within the nearest-neighbor model) resulting in complex phase diagrams, which makes the analysis of experimental data and the extraction of the model parameters a difficult task. In such cases, clarification of the underlying symmetry properties of the model is often of a great help. In general, the spin-orbital models in Mott insulators possess peculiar symmetries [6; 14] which are rooted in the bond-directional nature of orbitals. In this context, a special four-sublattice rotation [6] within spin space has proved itself as an extremely useful tool in the case of the original two-parameter KH model [11; 32; 47; 48; 49; 50]. It maps the Hamiltonian on itself but changes the Hamiltonian parameters, connecting thereby different points in the parameter space. Being an exact transformation, it transfers the complete knowledge about some point in the phase diagram, including the groundstate, excitation spectrum, response functions etc., to its partner. Based solely on this self-duality of the model, the entire phase diagram could be sketched and the deep relations between the phases understood. In addition, it also reveals points of hidden symmetry, where the system is exactly equivalent to a Heisenberg model for the rotated spins. Given its usefulness, it is highly desirable to find and analyze similar transformations for the extended versions of the KH model.

In this paper, we introduce a systematic method to derive dual transformations of bond-anisotropic spin Hamiltonians and demonstrate its results and their physical implications in the case of honeycomb iridates adopting the full nearest-neighbor model [36; 37; 39]. We find all the hidden -symmetry points of the model, the most peculiar one being characterized by a “vortex”-like pattern with a six-site unit cell, and demonstrate how the characteristics of the hidden Heisenberg magnet manifest themselves in the anisotropic situations. By identifying the points we characterize all the possible gapless Goldstone modes that may be encountered within the model. This is relevant in the context of real materials as the spin gap was found to be well below [28; 29], suggesting a connection to some of the points. Finally, using a self-duality of the model, we will provide a link between our fits of the earlier NaIrO data [32] and the recent experimental observation of the ordered moment direction [46]. We argue that this observation provides a direct access to the strength of the additional terms “extending” the KH model, and quantify the spin easy axis direction in terms of this “departure” from the pure KH model. This allows us to suggest plausible values of the model parameters that are compatible with the current data. While we focus here on the case of a honeycomb lattice as realized in NaIrO and more recently in RuCl [51], the method is general and expected to produce interesting results also in the context of the new structural families of iridates – recently synthesized hyper-honeycomb [52; 53] and harmonic-honeycomb lattices [54; 55], or the theoretically proposed hyperoctagon lattice [56].

The paper is organized as follows. Sec. II introduces the Hamiltonian and discusses its parameters. Sections III and IV introduce the method, derive and discuss the main results of the paper – the hidden symmetries of the model. Sec. V and Appendix B discuss the implications of the results for honeycomb iridates.

Ii Extended Kitaev-Heisenberg model

Figure 1: (Color online) (a) Top view of the honeycomb NaIrO plane, the definition of global axes, and the reference frame for the spin components. The - and -directions coincide with the crystallographic - and -axes. The three bond directions of the honeycomb lattice are labeled as , , and , its two sublattices are labeled by and . (b) Two edge-shared IrO octahedra of a -bond and the definition of the local spin axes , , (used in Eq. 1). (c) Simultaneous cyclic permutation of the Ir-Ir bond directions , , and the spin components , , when applying a rotation to the model.

We start by specifying the model Hamiltonian including all symmetry-allowed spin interactions on nearest-neighbor bonds. An ideal, undistorted, structure of the honeycomb NaIrO plane is shown in Fig. 1(a). We will utilize its rotational symmetry and the three sets of parallel mirror planes containing the shared edges of the IrO octahedra and cutting the Ir-Ir bonds into halves. The symmetry links the interactions for different bond directions while the presence of the mirror planes restricts the possible interactions for a given bond direction. A trigonal distortion (compression or elongation along the -axis) fully preserves these symmetries so that our Hamiltonian applies in that case as well. Furthermore, recent experiments [46] indicate a nearly ideal symmetry of the spin properties and hence suggest that additional terms, possibly induced by a monoclinic distortion present in NaIrO, can be neglected. Physically, this observation implies the robustness of the pseudospin wavefunctions against weak monoclinic distortions.

The bond Hamiltonian is most compactly expressed in a local, bond-dependent, reference frame for spins, presented in Fig. 1(b) for a -bond. Due to the mirror symmetry, the in-bond component is forbidden to interact with the and components [36]. Following the notation of Ref. 36 we arrange the allowed terms into the form

(1)

This four-parameter Hamiltonian extends the KH model (- and -terms) by the -term bringing further anisotropy among the diagonal components of the interaction, and the -term determining the only symmetry-allowed non-diagonal element in the exchange interaction tensor. Parameter would vanish for an isolated pair of undistorted octahedra; it becomes finite due to a trigonal distortion and/or due to the extended nature of orbitals in a crystal (“recognizing” the fact that the octahedra are canted relative to the crystal axis ).

To capture the symmetry, it is convenient to switch to cubic axes , introduced in Fig. 1(a) and pointing from an Ir ion to neighboring O ion positions in an ideal structure. The -bond Hamiltonian in the cubic reference frame, as derived in Ref. 37, reads then as:

(2)

with the correspondence and , often used below. For the other bond directions, the Hamiltonian is obtained by a cyclic permutation [see Fig. 1(c)], resulting in one-to-one correspondence between the three types of bonds and interactions, as required by symmetry. Physically, each type of bonds favors its own distinct “orbital setup” to optimize the hopping energy, and this is fingerprinted in pseudospin interactions via spin-orbit coupling. For completeness, the Appendix A shows the Hamiltonian in the global axes ; it has certain advantages moving the bond dependence from the operator forms to the coupling constants.

Few comments are in order concerning the model parameters. In general, calculation of exchange integrals in transition metal compounds with -- bonding geometry is an intricate task, because more hopping pathways are allowed as compared to a simpler case of -- bonding in perovskites (where theory [9] has correctly predicted the strength of dominant exchange constants). For instance, orbitals may also overlap directly, in addition to oxygen mediated hoppings; there is a large overlap between orbitals of and symmetries (forbidden in perovskites), etc., resulting in a number of competing ferromagnetic and antiferromagnetic contributions which are difficult to evaluate, in particular in compounds with small Mott and/or charge-transfer excitation gaps. The uncertainties in interaction parameters and further affect the theoretical estimates.

Initial consideration [6] of the pseudospin one-half exchange interactions in -bonding geometry resulted in (hitting a “hidden” point by chance) in the cubic limit; later work [8; 9; 11] using different approximations has changed this estimate both in terms of the signs and values of and , illustrating the difficulties described above. It was also found that the non-diagonal element allowed in cubic symmetry may take sizable values [37; 36; 39; 40]. Further, is expected to become as large as the other parameters if trigonal splitting of the orbital level, caused by a compression along -axis, becomes comparable to spin-orbit coupling ; also, the trigonal field suppresses the parameter . These trends are easy to understand: large trigonal field suppresses the in-plane components of orbital moment and , leaving the axial component the only unquenched one; thus the pseudospin one-half Hamiltonian, written most conveniently in global axes in this limit, may not contain anything but and type terms: , identical for all bonds. This is what has indeed been found by explicit calculations [6; 35] in the limit of . This implies and in this limit (see also Appendix A), while and may take any values depending on the microscopic details. Although this limit is not very realistic for Ir ion with large spin-orbit constant [3; 57], we may expect sizable values of both and in NaIrO where seems to exceed [58; 59]. The role of and terms should further increase in other compounds based on pseudospin one-half Co, Ru, and Rh ions with smaller .

In general, the high-energy behavior of spins and orbitals in transition metal compounds is well captured by the Kugel-Khomskii models [4] and their descendants [6]. However, the low-energy physics and ultimate magnetic “fixed-point” are heavily influenced by many “unpleasant” details originating from orbital-lattice coupling and distortions, unavoidable in real materials. In perovskites, the Kugel-Khomskii energy scale is given by independent on spin-orbit coupling; however, this leading term drops out for pseudospins one-half in the edge-shared, -bonding geometry [6; 9], so the “high-energy” scale is set up by the subleading terms. In iridates, the hope [9] is that the Kitaev-type coupling is the leading one among these subleading terms. Since this coupling is itself a correction to , this expectation may or may not hold in reality.

To summarize up to now: in real materials even with an ideal symmetry, all the four exchange parameters may play a significant role. This motivates us to regard the Hamiltonian (1,2) as an effective model with arbitrary parameters, and look for some general symmetry arguments that may help to identify plausible parameter windows in the analysis of experimental data.

Iii Systematic construction of dual transformations

Having fixed the model Hamiltonian, we are ready to explore its dual transformations. By a dual transformation we mean a prescription for site-dependent rotations in the spin space, , which transforms a spin Hamiltonian into a formally new Hamiltonian . We are interested in self-dual transformations of that map the model onto itself, preserving its all symmetry properties. That is, the rotated partner (i) has the same four terms albeit with different parameters , and (ii) it respects the rotation rules encoded in Fig. 1(c), hence preserving the original distribution of the three types of bond-dependent interactions on a lattice.

Starting with the Hamiltonian expressed as where are matrices, we obtain with . For a self-dual transformation, the matrices are identical to , but the parameters are replaced by , and the one-to-one correspondence between the bond directions and interactions remains intact. These two points in the parameter space are linked by the transformation and knowing the solution at one of the points, we may “rotate” it to the other one.

In this section we give an algorithm to find the self-dual transformations for the extended KH model that map it onto itself. We have found a single self-dual transformation operating in full parameter space of the model; we will show it shortly below and return to it later when discussing experimental data.

However, studying the hidden symmetries of the model, we have identified a number of restricted self-dual transformations that operate only in some regions of the parameter space, where constants are all finite but obey certain relations, or some of them are simply zero. Our primary interest is in the special class of such transformations of the type , which convert the Heisenberg model into the full model and vice versa. These transformations, to be discussed in the next section, reveal points of hidden symmetry – by inverting the transformation the anisotropic model with the parameters can be exactly mapped back to the Heisenberg model with the exchange constant .

iii.1 Algorithm

A systematic search for the dual transformations seems to be an intricate task. Fortunately, it can be easily performed by computer on a finite cluster of the lattice using the following simple algorithm. We give it specifically for the case of a self-dual transformation:

A) As a first step, we choose two rotation matrices , on a selected bond . They have to preserve the -form given by (1), which leaves us with only a few choices, each having only one free angular parameter.

B) Next, we randomly choose nonzero values of the initial parameters and use the relation together with the symmetry to determine the new Hamiltonian matrices for the three bond directions.

C) Knowing all the bond Hamiltonians, we may now determine further rotation matrices by utilizing relations of the type and proceeding neighbor-by-neighbor. To fully determine the rotation matrices, about two thirds of the bonds need to be used.

D) The bonds of the remaining third are used to check consistency, the Hamiltonian matrix determined by using the rotation matrices belonging to the bond has to be identical to that determined in step B. If the total difference on all the remaining bonds equals zero, we have just constructed a self-dual transformation. By scanning through the entire interval of the free parameter introduced in step A, we find all the self-dual transformations.

The above procedure may be easily adapted to find the dual transformations such as . In this case, in step A of the algorithm, we use the symmetry of the Heisenberg model and choose as an identity matrix. The choice of the second matrix is restricted by the requirement that is of the form.

By inspecting the rotation matrices of the cluster, we can identify the particular unit cell of the transformation. Note, that even if our cluster is smaller than this unit cell, we do not miss the corresponding transformation, so that the method is completely systematic [60].

iii.2 Self-duality of the extended Kitaev-Heisenberg model

The systematic procedure described above has identified only a single self-dual transformation . This is not surprising given the complexity of the model. The corresponding parameter transformation may be written in a matrix form

(3)

For convenience, we also give the transformation of the parameters entering the Hamiltonian (2):

(4)

In terms of the spins, the transformation, labeled for future reference as , is simply a global -rotation about the -axis defined in Fig. 1(a). The individual , , and components transform according to

(5)

at every site. By applying the transformation twice, we get an identity and the matrices in (3,4) are thus self-inverse. Despite its apparent triviality, this transformation will play an essential role when discussing the real materials, see Sec. V below.

Iv Points of hidden SU(2) symmetry

In this paragraph we find and characterize all the points of hidden symmetry present in the extended KH model. At these special points in the parameter space, the anisotropic model can be mapped back to a Heisenberg ferromagnet or antiferromagnet. The points of the original KH model have been identified [11] by virtue of the four-sublattice transformation introduced in Ref. 6. The corresponding ordering patterns on the honeycomb lattice are of stripy and zigzag type. A similar symmetry analysis of the KH model was performed for other relevant lattices [47].

The extended KH model of course inherits the points of the KH model and contains several new ones in addition. They are identified by dual transformations of the type which is less general than . Because of this, we obtain a relatively rich set of dual transformations characterized by two-, four-, and six-sublattice structure of the rotations. In terms of parameters, all the non-trivial points are listed in Table 1. We now proceed with the detailed description of the corresponding transformations.

Table 1: Parameter values for the points in units of the exchange constant of the hidden Heisenberg model.

iv.1 Summary of the SU(2) points and the corresponding rotations on the sublattices

We first give a summary of the transformations as represented by rotations in the real space. Each of them generates an infinite number of orderings, since the ordered moment direction in the underlying Heisenberg model can be chosen arbitrarily. Figure 2 shows a few important examples.

Figure 2: (Color online) (a,b) Unit cells for the four- and six-sublattice transformations. (c,d) Stripy and zigzag patterns related to the FM and AF order of a hidden Heisenberg magnet via the four-sublattice transformation . The spins take the -axis direction. (e,f) “Vortex”-like patterns generated by the six-sublattice transformation . The spins are lying in the lattice plane in the case presented. The colors of the arrows in panels (d) and (f) indicate the sublattices of the hidden AF order. (g) Brillouin zones of the honeycomb (inner hexagon) and the completed triangular lattice (outer hexagon). The characteristic vectors of the four-sublattice transformation and of the six-sublattice transformation are shown in red and blue, respectively. (h) Bragg spots of the patterns in panels (c-f). The dot size is proportional to .

The simplest transformation is -rotation about -axis at one of the two sublattices of the honeycomb lattice:

(6)

Its physical relevance is small due to the dominance of and the complete absence of (corresponding to the case of strong trigonal field splitting, as explained above). As a curiosity, if we choose the spins to lie in the honeycomb plane, converts the FM pattern to AF and vice versa. We may thus have an AF/FM ordered pattern but the hidden nature revealing itself e.g. in the spin dynamics is that of Heisenberg FM/AF, respectively.

The next transformation has a four sublattice structure depicted in Fig. 2(a) with -rotations about cubic , , and axes applied at sublattices , , and , respectively, and no rotation involved at sublattice . Written explicitly:

(7)

This transformation, introduced earlier in Ref. 6, is a self-dual transformation of the original two-parameter KH model and has been already heavily used in this context. Applying the transformation to an ordered Heisenberg FM/AF with the moments pointing along the -axis, we get the stripy/zigzag order shown in Fig. 2(c,d).

Perhaps the most surprising point of the model is linked to the six-sublattice transformation . Its rotations are most conveniently described in the cubic coordinates. On the lattice sites , , and [see Fig. 2(b)] they correspond to cyclic permutations among the spin components. On the lattice sites , , and the rotations correspond to anti-cyclic permutations which have to be followed by a spin-inversion. Altogether the transformation can be written as

(8)

It is easy to see, that for and , these rotations lead to the isotropic Heisenberg Hamiltonian. As an example, we consider the -bond - of Fig. 2(b). By exchanging and at site , the non-diagonal -term in (2) becomes diagonal and the inversion ensures its proper sign. Sample patterns generated by and showing a “vortex”-like structure are presented in Fig. 2(e,f). The peculiarity of the points is now best demonstrated: the Hamiltonian is completely anisotropic containing and terms only, the ordered spins form a very unusual pattern, yet the hidden nature of the system is exactly that of the Heisenberg FM or AF, including e.g. the presence of gapless Goldstone modes.

Apart from revealing a hidden point of the present model, the transformation has a remarkable property that deserves a special attention. Namely, applying to the Kitaev Hamiltonian, we notice that it re-distributes three types of Ising-interactions on a honeycomb lattice such that at each hexagon a Kekulé type pattern is formed [61]. We thus arrive at the so-called Kekulé-Kitaev model [62]. In other words, the Kitaev and Kekulé-Kitaev models are exact dual partners linked via the transformation. This observation should be helpful in studying both models, in particular of their extended versions including a Heisenberg term [62; 63].

Two more transformations providing points are obtained as the combinations and . They share the sublattice structure with and , respectively. Adopting the extended KH model, the former one is probably the point closest to the real situation in NaIrO as will be discussed in Sec. V.

iv.2 Implications for the phase diagram

After examining the nature of the individual points, we want to visualize now their positions in the parameter space, get a sketch of the phase diagram, and infer the relations between the individual phases. The result can be compared with the published phase diagrams of Refs. 37 and 38, obtained by classical analysis and partly complemented by exact diagonalization. For this reason, we adopt the representation of the parameter space introduced in Ref. 37. The overall energy scale irrelevant for the phase diagram is removed and , , are parametrized using “spherical” angles and via , , and , keeping as a separate parameter of the phase portrait.

Figure 3: (Color online) (a) Depiction of the points using the parametrization of Ref. 37, , , . The distance from the center of the circle corresponds to going from (center) through (dashed circle) to (solid circle). The polar angle is . Filled squares show the points with , open squares those with nonzero values given on the right along with the transformation label. The color of the points indicates their hidden FM (blue) or AF (red) nature. The green square (circle) shows the parameter values specified in Sec. V when discussing NaIrO (LiIrO). (b) The same as in panel (a) but with .

Shown in Fig. 3 is the complete set of points of the extended KH model. The outer rings correspond to the original KH model and contain the trivial points and the two well-known hidden points of the KH model characterized by a stripy and zigzag pattern. Still within the plane is the “vortex” point associated with a “vortex”-like pattern. The corresponding phases determined by these points can be observed in Figs. 2 and 3 of Ref. 37, with the point lying in their phase.

Three more points , , and characterized by a nonzero value of are shown as projected onto the plane. For the case presented in Fig. 3(a), they are of AF character, one of them appears for positive (point ) and two for negative (points and ) values of . The point (given by ) of hidden AF nature can possess FM pattern as discussed in the previous paragraph. The region between the true FM Heisenberg point and the point in the phase diagram obtained classically is therefore filled by the FM phase extending as increases [see panels (c) and (e) of Fig. 2 of Ref. 38]. However, the (hidden) nature of this phase changes from FM to AF which should manifest itself e.g. on the character of the magnon dispersion. Similarly, the presence of the points () and () of zigzag and “vortex” character, respectively, explains the enlarged region of the corresponding phases in the classical phase diagram for [see the panels (a) and (d) of Fig. 2 of Ref. 38]. We also observe an intimate relation between the zigzag phase emanating from the () point and that connected to the zigzag point of the original KH model (given by ). Due to the additional rotation, their ordered moment directions are related by -rotation about the global -axis. This point will be further discussed in Sec. V. Finally, similar conclusion as for the case presented in Fig. 3(a) can be drawn for the case shown in Fig. 3(b). The points are related by inversion with respect to the center of the circle and the opposite FM/AF nature.

In summary, we have illustrated that the gross features of the phase diagram of the extended, four-parameter KH model can be deduced solely by inspecting the nature of the points of hidden symmetry and their location in the parameter space.

iv.3 Spin excitation spectra

We proceed further by inspecting the spin excitation spectra at the points associated with and transformations, and see how they are related to those of the simple Heisenberg model. To this end, the dual transformations have to be expressed in Fourier space and relations between the Fourier components of the dual partners have to be established. The situation is somewhat complicated by the two-sublattice structure of the honeycomb lattice, requiring us to introduce an additional index [see the labels and in Fig. 1(a) for the convention used below].

In both cases, it is convenient to use the cubic axes . The four-sublattice transformation has three characteristic vectors and touching the Brillouin zone boundary in the middle of its edges [see Fig. 2(g)]. The rotation matrices have a simple diagonal form, reflecting only the sign changes of the respective components

(9)

The six-sublattice transformation written in Fourier representation has a full matrix structure

(10)

with the factor and the matrices

(11)

where . The characteristic vectors and shown in Fig. 2(g) again touch the boundary of the Brillouin zone, now in its corners. The dual transformation takes a general form (here for sublattice ) which translates into

(12)

i.e., the Fourier components get shifted by the characteristic vectors. As a side result, the above relation gives the Bragg spots derived from the Bragg spots of Heisenberg FM/AF () and presented in Fig. 2(h).

To study the spin excitations, we employ the spin susceptibility tensor defined as

(13)

It is evaluated at the points by first decomposing into the and -sublattice contributions via

(14)

applying the dual transformation in the Fourier form of Eq. 12 to get back to the underlying Heisenberg model, and using the spin susceptibility for the Heisenberg model obtained within linear-spin-wave (LSW) approximation. In the case of , this brings simple -shifts by , , and for the individual components. For , the corresponding expressions are somewhat more involved containing a non-shifted contribution and shifted contributions combining pairs of the characteristic vectors , , and . Without going into details, the presence of both shifted and non-shifted parts can be easily understood based on Eq. 10.

Figure 4: (Color online) (a) LSW dispersion of the Heisenberg FM (blue) and AF (red) on the honeycomb lattice. The width of the lines indicates the trace of the spin susceptibility tensor, , calculated in the LSW approximation. (b) The same for the stripy and zigzag state presented in Fig. 2(c,d). Energy is scaled by of the hidden Heisenberg magnet. (c) The same for the “vortex”-like patterns presented in Fig. 2(e,f).

Presented in Fig. 4 are the traces of the spin susceptibility tensor of the Heisenberg model and the extended KH model at the two hidden points under consideration. For completeness, we demonstrate both hidden FM and AF case characterized by quadratic and linearly dispersing Goldstone modes, respectively. The situation is more transparent for the four-sublattice patterns – stripy (hidden FM) and zigzag (hidden AF) – since the spinwave dispersions are just shifted with the = points replacing the Goldstone points and of the Heisenberg case. In our example, we have chosen -axis as the ordered moment direction. For the magnons, which are in fact deviations of the ordered moment in and directions, only and shifts are active, selecting four out of the six -points in total. The remaining two are the Bragg spots reached from and by shifts active for the ordered spin component. The Bragg spots and the Goldstone points are thus complementary in this case. The spin excitations associated with the six-sublattice patterns are significantly more complicated. They contain both shifted Goldstone modes [in Fig. 4(c) such a mode appears at =-point coinciding with ] and Goldstone modes at the characteristic momenta and of the underlying Heisenberg model. In the latter case just the intensity of the modes has been transferred by the dual transformation, making e.g. the linear Goldstone mode at the most intense one in the hidden AF case.

A similar analysis of the spin excitations as presented here for and points can be performed for the remaining points. Due to the nature of the relevant transformations, no other characteristic vectors appear. Therefore, , and their counterparts shifted by the vectors and entering the transformations and constitute the entire set of the wavevectors of the Goldstone modes that can be observed within the extended KH model.

V Application to the real materials

The aim of this work was to study the basic symmetry properties of the extended KH model – a promising spin Hamiltonian for the magnetism of honeycomb iridates. Below we illustrate how this knowledge, taken together with the experimental data, helps to locate the plausible windows in otherwise very large parameter space even for this nearest-neighbor (NN) model. We will show that, despite having only a single result, the search for full self-dual transformations of the extended KH model provides us with a surprisingly useful tool in the context of NaIrO. This utility of emerges due to the recent observation of the magnetic moment direction [46] which, as we see shortly, imposes an important constraint on the model parameters. This is because, in general, the data on magnetic easy axes in a crystal, along with the magnon gaps and torque magnetometry data, provides a direct information on the symmetry and strength of the anisotropy terms in spin Hamiltonians, and the case of NaIrO is of course not any special in this sense.

To begin with, we recall that NaIrO shows so-called zigzag order, where the spins on and bonds are parallel and form ferromagnetic chains that run along direction and couple antiferromagnetically along the -axis. This relatively simple collinear magnetic structure has been first explained [30; 28] as due to -NN and -NN Heisenberg couplings (which are often relevant in compounds with -bonding geometry – well known example is quasi-one dimensional cuprates). This model emphasizes a geometrical frustration which is realized at large values of and resolved by the symmetry breaking zigzag formation.

However, as argued in the Introduction, more recent data [45; 46] suggests that the origin of frustrations is largely related to the non-Heisenberg-type interactions which are bond-dependent and hence highly frustrated even on the level of NN-models. A minimal NN-model of this sort is the KH model, which has been shown [32] to host zigzag order in its phase diagram indeed. We follow this way of reasoning and explore below the extended version of the KH model as the basic NN-model for iridates. On the way, we will also see the point where the data may require the presence of additional terms too, suggesting that the both “zigzag theories” above are the part of a full story.

In Ref. 32, the available experimental data on NaIrO has been fitted using the two-parameter KH model, regarding it as a phenomenological spin Hamiltonian with arbitrary parameters. For and , the model was found consistent with experiments in terms of the type of magnetic ordering, temperature dependence of static magnetic susceptibility and the low energy spin-excitation spectrum being compared to powder INS. Later, RIXS experiments [45; 46] confirmed the presence of a high energy branch of spin excitations, with an even better agreement obtained if the LSW calculation of Ref. 32 is replaced by a more suitable exact diagonalization [65].

However, the recent data [46] on the moment direction came about as an unexpected surprise, challenging at first glance the above coherent description of NaIrO. The point is that within the original two-parameter KH model, the zigzag order is characterized by the spins pointing towards one of the oxygen ions [see Fig. 1(a)]. This expectation is generic and guaranteed by the “order-from-disorder” physics [66] which typically selects one of the high-symmetry cubic axes as the easy one, when a spin Hamiltonian contains the compass-type or Kitaev-type bond-dependent anisotropy [16; 11; 40; 67], independent on parameter values. The resonant magnetic x-ray diffraction data [46] shows instead that the magnetic easy axis is in fact far away from any of the Ir-O bond directions: it is oriented “nowhere” slightly below a midpoint between the two, and , oxygen ions in Figs. 1(a) and 5(a). This is a clear indication of the significance of the and terms in the spin Hamiltonian [68].

To reconcile all the data at hand using now full four-parameter model, we first notice that the above two easy axis directions – the one observed in NaIrO and the one expected from the KH model as used in Ref. 32 – are roughly related to each other simply by a -rotation about the -axis. This observation gives an immediate hint how to obtain a starting parameter point when fitting the current data set for NaIrO within the extended, four-parameter KH model in an appealingly easy way, and resolve the above apparent problem with the moment direction.

As discussed in Sec. IV.2, the -associated zigzag phase of the KH model with is related to the zigzag phase connected to the point () of Fig. 3(a) via . Remarkably, due to the nature of – a global -rotation of the magnetic moments about -axis – all the aforementioned consistent results [32] of the two-parameter KH model are fully preserved if we apply Eq. 3 to the parameters and of Ref. 32 given above, the only change is the spin easy axis being rotated to the proper direction as in experiment. The corresponding set of parameters obtained via (3) is: , , , . We would like to emphasize that these numbers should not be taken literally; rather, they fix the signs of the parameters involved and put an upper limit for and , as we explain below.

Figure 5: (a) Pseudospin angle relative to the -plane (see inset) as a function of the parameter . Dashed lines show the “magic” angle and its complement , determined by the -axis and -plane, respectively, as sketched in the inset. (b) The phase diagram as a function of long-range couplings and anisotropy parameter . Starting with the “bare”, -derived values of , , , and , we have scaled and simultaneously to vary . To stay within the zigzag phase at the smaller values of , one needs to have finite couplings. Otherwise, the NN-only extended KH model with negative switches to the incommensurate and “stripy” [11; 69; 32] ground states. The inset shows the exchange bonds and .
Figure 6: (a) Map of the -dependent classical energy (per site, in units of meV) obtained by Luttinger-Tisza method [71; 72] for the “bare”, i.e. -derived parameters , , , relevant for NaIrO. The hexagon indicates the first Brillouin zone. (b) The same for the parameters , relevant to LiIrO. (c) Length of the ordering vector for varying , keeping the other parameter values unchanged. The dashed (solid) line was calculated using the above parameters relevant to NaIrO (LiIrO). Points and show the values used in panels (a) and (b), respectively. (d) LSW dispersions for the parameters used in panel (a) (left) and panel (b) (right). In the latter case, we have taken instead of to stay in the zigzag phase at the border to the spiral phase.

For the representative parameters given above, the pseudospin makes a “magic” angle of from the -plane, as follows from construction. This is slightly lower than observed [46; 70]. Now, we inspect how the angle varies as a function of the anisotropy parameters. The result is illustrated in Fig. 5(a) and shows that the exact value of heavily influences the “departure” from the KH model quantified by the parameter . (The corresponding Eqs. 32, 33, and 35 for the spin direction are derived in the Appendix B by minimizing the classical energy). The “magic” angle appears at – as obtained for the above parameters. Yet, as observed in Fig. 5(a), by increasing for example by only, we already find and get closer to the regime. Therefore, more detailed measurements and fits of the ordered spin direction are highly desirable to get the actual values of the parameters and relative to the Kitaev term . Doing so, it is crucial to take into account the fact that the pseudospin direction and magnetic moment direction are not the same in general; while they coincide in the cubic limit, a sizable trigonal-field splitting might be present in NaIrO [45; 59]. It is thus important to quantify this splitting by independent measurements.

At this point, longer-range couplings become a part of the full spin model for iridates, for the following reason. As a partner of the zigzag phase of Ref. 32, the present NN-model with large is well in its zigzag ordered state. But this is not so at smaller values of (e.g. for ), which are required to get the spin angles or above, see Fig. 5(a). Incorporating moderate and couplings into the model, we can however stabilize the zigzag phase, see Fig. 5(b), and hence obtain the ordered spin angles above the “magic” one. The values of of the order of are indeed suggested by ab-initio calculations [39]. This shows again the key importance of the experimental data on moment directions for quantifying the balance between the two zigzag-supporting mechanisms discussed above: based on geometrical frustration, and on frustration driven by the non-Heisenberg nature of interactions in spin-orbit coupled magnets. Recent observations [46] of a pronounced spin-space anisotropy on one hand, and an “intermediate” spin direction that requires finite values on the other hand, suggest that both mechanisms are at play in NaIrO.

Altogether, the present analysis using the symmetry properties of the model, taking into account the recent data on moment direction [46], as well as considering the role of the couplings suggests a plausible window in the parameter space of an effective spin model for NaIrO: , with positive (AF) Heisenberg couplings and . The leading anisotropy terms and are both negative, while a smaller term may in principle take any sign. This parameter window is globally consistent with experimental observations on NaIrO we are aware of to date, and may be used as a guide in future analysis, in particular once -resolved spin response becomes available, and the ordered pseudospin and magnetic moment directions (they differ in general) are obtained and confirmed by independent measurements.

Even though this general result still leaves quite a freedom, it is of great help by fixing the signs of most relevant couplings and their hierarchy. This is the main outcome of the present theory in the context of real materials. Further, we note that the Kitaev coupling can be deduced from overall magnetic energy scale, and the spin and magnetic moment directions should determine the parameter hence . From a careful analysis of the zigzag stability condition, magnon gaps and dispersions, paramagnetic susceptibility data, etc., one should be able to quantify all the model parameters including , , and .

Considering this result in the context of microscopic theories, we notice first that the signs of and above are consistent with the original calculations of these parameters for honeycomb iridates [9; 11] as well as with the later studies [36; 37; 38; 39; 40]. Next, we may conclude that a contribution from hopping that favors pseudospin interaction with [6] is not significant in iridates; this is also consistent with the recent calculations [73; 36]. Further, the present symmetry analysis resolves an apparent conflict with the theoretical [9; 11] and the positive that follows from the best data-fit using the KH model [32]: in fact, the pure KH model with and the extended one with and sizable terms are -dual partners (the latter one being physical).

More surprisingly, relatively large anisotropy term is required to “turn” the moment direction well away from the pure KH model position. A positive implication of this observation is that this term makes it much easier to stabilize the zigzag order (the pure KH model with large would require large long-range couplings otherwise). In a view of the discussion in Sec. II, this suggests a presence of sizable trigonal field effects in NaIrO. Eventually, an unusual – out of any crystal symmetry axis – orientation of pseudospins [46] should originate from a competition among the several anisotropy terms , , and of different symmetry and physical origin.

To conclude our discussion of NaIrO: it seems that the extended KH model, likely further “extended” by moderate longer-range couplings, is indeed a good candidate model for this compound. Even though these extensions (to be still quantified by future experiments) reduce the chances for “pure” Kitaev-model physics in iridates, the model itself is highly interesting due to its rich internal structure and hidden symmetries that we have uncovered in this work.

Motivated by the above, we further consider the case of LiIrO. Since the data is limited here, the discussion will be brief and suggestive only. Due to the smaller Curie-Weiss temperature and more “ferromagnetic” behavior of its spin susceptibility [74; 75], this compound was located closer to point of the KH model [32]. Even though the parameters and given in Ref. 32 correspond to the zigzag phase while LiIrO shows a spiral magnetic ordering [29], these parameters can be used to get a hint of the direction in the parameter space to consider. We therefore transform the above parameters using (3) to obtain: , , , . Representing the parameters for Na and Li compounds obtained via transformation (3) in Fig. 3, we see that both are close to the point (), with Li being closer, as expected. To approach the spiral state observed in LiIrO, we first note that, in first approximation, in both cases and that Li compound is characterized by a much smaller . For simplicity, we set , assume to roughly preserve the overall energy scale, and reduce the parameter associated with the trigonal distortion, which is expected to be much smaller in LiIrO with the bond angles being closer to . The Luttinger-Tisza [71; 72] maps of the classical energy for the Na and Li case presented in Fig. 6(a,b) confirm the zigzag and incommensurate magnetic ordering, respectively. For the incommensurate ordering wavevector is obtained as [see Fig. 6(b,c)], this would predict a magnetic Bragg peak in powder neutron diffraction experiments at a -value that could be consistent with experiments on powder LiIrO [29]. Finally, Fig. 6(d) compares the spin excitations obtained using LSW approximation. In the case of NaIrO, the dispersion is identical to that presented in Fig. 3 of Ref. 32, possessing a low and high-energy branches. As the parameter is reduced, these two branches gradually merge, leading to a steeper dispersion compared to NaIrO, which might be consistent with powder inelastic neutron scattering experiments on powder LiIrO [29]. The predicted dispersion is illustrated in Fig. 6(d) for a point on the boundary between the zigzag and the spiral phase. A further minor reduction of to enter the spiral phase and get the proper ordering vector should not affect this result dramatically, apart from the changes at low energies forming an “hour-glass” shape characteristic to spiral magnets (see, e.g., Ref. 76).

Vi Conclusions

To summarize, we have analyzed non-trivial symmetries of the extended Kitaev-Heisenberg model on the honeycomb lattice. As a main result, we have identified the complete set of points in the parameter space where this bond-anisotropic model can be transformed to a simple Heisenberg model and is therefore characterized by hidden symmetry. Such a dual transformation can be performed using a particular choice of sublattice rotations of the spins, specific for each of the points. The sublattice structure of the transformations creates a number of ordering patterns which together with the location of the hidden points in the parameter space give a good overview of the global phase diagram of the model. In terms of the spin excitations, the hidden symmetry manifests itself by the presence of Goldstone modes inherited from the symmetric Heisenberg FM/AF on the honeycomb lattice. Their characteristic vectors and even the full spin excitation spectra are easily obtained by an explicit transformation of the FM/AF case.

One of the special transformations linked to the hidden points reveals at the same time an exact duality between the Kitaev and Kekulé-Kitaev models; this result should be useful in theoretical studies of these and related models.

We emphasize that, adopting the extended KH model, all the above results are necessary consequences of its symmetry which is in turn dictated by the underlying symmetry of the lattice.

Having the results of the general symmetry analysis at hand, we were able to find the region of the parameter space that is consistent with the observed properties of the honeycomb lattice iridates NaIrO and LiIrO. Further, a relation between the ordered moment direction and the model parameters is derived, which may help to quantify these parameters from future experiments.

Finally, our method to systematically explore the hidden symmetries is general and can be applied to other bond-anisotropic models as well. In the context of the iridate materials, the symmetry analysis of the extended KH model on hyper-honeycomb and harmonic honeycomb lattices is of a great interest.

Acknowledgements.
We would like to thank B.J. Kim for sharing with us the experimental data [46] which motivated this study, R. Coldea, G. Jackeli, I. Kimchi, N.B. Perkins, S. Trebst, and S.E. Sebastian for helpful discussions and comments. JC acknowledges support by ERDF under project CEITEC (CZ.1.05/1.1.00/02.0068), EC 7 Framework Programme (286154/SYLICA), and Czech Science Foundation (GAČR) under project no. 15-14523Y.

Appendix A XYZ form of the Hamiltonian

The Hamiltonian expressed in terms of the spin components , , and , corresponding to the reference frame in Fig. 1(a), takes the form

(15)

Here the symmetry of the model is embodied in the factors and , where the angles are determined by the bond directions: for the , , and bonds, respectively. In terms of the original parameters , the exchange constants entering (A) read as

(16)
(17)
(18)
(19)

Note that it is the and terms which bring about the bond-directionality of the interactions, and hence they naturally support -symmetry breaking orderings such as zigzag in the present model. Physically, these terms arise from the exchange processes that involve the in-plane components of orbital momentum and which “know” the bond directions, like the orbitals do in the Kugel-Khomskii models.

It is also noticed that the and terms change the -component of total angular momentum by and , correspondingly. This is because the -orbital angular momentum is not a conserved quantity in a crystal, and this commonly shows up in effective spin Hamiltonians due to the spin-orbit coupling.

A strong trigonal field splits -level such that the lowest Kramers doublet (pseudospin) wavefunctions become simple products of and spin , states, correspondingly; i.e., there will be a one-to-one correspondence between the real spin and pseudospin directions. Since the total spin is conserved during the hoppings, pseudospin is then conserved, too. Thus, the spin non-conserving terms and must vanish in this limit, which implies and simultaneously. Physically, a strong compression along trigonal axis dictates that this axis becomes the “easy” (or “hard”) one for moments. Since this limit is not realized in iridates, we will not use the -form of Hamiltonian in this paper; however, it might be useful for pseudospin one-half Co, Rh, and Ru compounds where the spin-orbit and crystal field effects may strongly compete.

Appendix B Analysis of the classical energy and moment direction

In this appendix we show the expressions used in the classical energy analysis. We first give the Hamiltonian in its momentum space form utilized within the Luttinger-Tisza method. By minimizing the classical energy in the zigzag phase we then find the ordered moment direction.

Transforming the spin operators via and similarly for the -sublattice, we cast the Hamiltonian to the form

(20)

where the -vectors cover the first Brillouin zone of the triangular lattice of . The simplest expressions for the matrices of the momentum-space Hamiltonian (20) are obtained using the cubic axes , , . Complementing the interactions in Eq. 2 by long-range and , we arrive at

(21)

and

(22)

Here the momentum-dependent factors read as

(23)
(24)
(25)
(26)

In the Luttinger-Tisza method [71; 72], the matrices for running through the Brillouin zone are diagonalized. The -vector and the eigenvector corresponding to the minimum eigenvalue then determine the ordering resulting on a classical level and the minimum eigenvalue itself gives the classical energy per site. This approach relaxes the spin-length constraint which should be checked afterward.

Next, we evaluate the classical energy for the zigzag state with the ordering vector and an arbitrary ordered moment direction given by a unit vector . The corresponding zigzag pattern is captured by . Using (20), we get for the classical energy per site:

(27)

with the matrix

(28)

or equivalently