Order-by-disorder and spin-orbital liquids in a distorted Heisenberg-Kitaev model

Order-by-disorder and spin-orbital liquids in a distorted Heisenberg-Kitaev model

Eran Sela Raymond and Beverly Sackler School of Physics and Astronomy, Tel-Aviv University, Tel Aviv, 69978, Israel    Hong-Chen Jiang Department of Physics, University of California, Berkeley, California 94720, USA    Max H. Gerlach Institute for Theoretical Physics, University of Cologne, 50937 Cologne, Germany    Simon Trebst Institute for Theoretical Physics, University of Cologne, 50937 Cologne, Germany
July 6, 2019

The microscopic modeling of spin-orbit entangled Mott insulators such as the layered hexagonal Iridates NaIrO and LiIrO has spurred an interest in the physics of Heisenberg-Kitaev models. Here we explore the effect of lattice distortions on the formation of the collective spin-orbital states which include not only conventionally ordered phases but also gapped and gapless spin-orbital liquids. In particular, we demonstrate that in the presence of spatial anisotropies of the exchange couplings conventionally ordered states are formed through an order-by-disorder selection which is not only sensitive to the type of exchange anisotropy but also to the relative strength of the Heisenberg and Kitaev couplings. The spin-orbital liquid phases of the Kitaev limit – a gapless phase in the vicinity of spatially isotropic couplings and a gapped Z phase for a dominant spatial anisotropy of the exchange couplings – show vastly different sensitivities to the inclusion of a Heisenberg exchange. While the gapless phase is remarkably stable, the gapped Z phase quickly breaks down in what might be a rather unconventional phase transition driven by the simultaneous condensation of its elementary excitations.

75.10.Jm, 71.20.Be, 75.25.Dk, 75.30.Et

I Introduction

The intricate interplay of electronic correlations, spin-orbit coupling, and crystal-field effects in transition metal oxides has led to the discovery of an intriguing variety of quantum states of matter including Weyl semi-metals, axion insulators, or topological Mott insulators Witczak-Krempa et al. (2014). In the correlation dominated regime unusual local moments such as spin-orbit entangled degrees of freedom can form and whose collective behavior gives rise to unconventional types of magnetism including the formation of quadrupolar correlations or the emergence of so-called spin liquid states.Balents (2010) On the materials side a particularly prolific group of compounds are the Iridates, whose electronic state can be either weakly conducting or insulating. Common to all Iridates is that the Iridium ions typically occur in an Ir ionization state corresponding to a electronic configuration. For the insulating compounds a particularly intriguing scenario is the formation of a so-called Mott insulator Kim et al. (2008, 2009), in which a crystal field splitting of the orbitals into t and e orbitals and a subsequent spin-orbit entanglement leads to a Mott transition yielding a completely filled state and a half-filled doublet. The microscopic exchange between these spin-orbit entangled local moments has been argued Jackeli and Khaliullin (2009); Chaloupka, Jackeli, and Khaliullin (2010) to give rise to interactions which combine a spin-like contribution in form of an isotropic Heisenberg exchange with an orbital-like contribution in form of a highly anisotropic exchange whose easy axis depends on the spatial orientation of the exchange path. Such orbital exchange interactions are well known from the early work of Kugel and Khomskii Kugel and Khomskii (1972) on quantum compass models Nussinov and van den Brink () to induce a high level of exchange frustration, i.e. they inhibit an ordering transition of the local moments which cannot simultaneously align with all their nearest neighbors due to the competing orientations of the respective easy axis. This frustration mechanism is particularly effective in the so-called Kitaev model Kitaev (2006), a honeycomb compass model where the exchange easy axis points along the , , and -directions for the three different bond orientations in the honeycomb lattice, see Fig. 1(b). Its phase diagram parametrized in the relative coupling strength of the three types of exchanges exhibits two incarnations of spin liquid phases: an extended gapless spin liquid phase around the point of equally strong exchange interactions and gapped spin liquid phases if one of the three coupling strengths dominates, see Fig. 1(c) for a detailed phase diagram. On the materials side, the layered Iridates NaIrO and LiIrO, which form Mott insulators with the Iridium ions arranged on a hexagonal lattice as illustrated in Fig. 1(a), have recently attracted considerable attention as possible solid state incarnations Jackeli and Khaliullin (2009); Chaloupka, Jackeli, and Khaliullin (2010); Singh and Gegenwart (2010); Singh et al. (2012); Choi et al. (2012); Ye et al. (2012) of the Heisenberg-Kitaev model 111Generalizations of the Heisenberg-Kitaev model to lattice geometries beyond the hexagonal lattice have recently been considered in both two and three spatial dimensions Rousochatzakis et al. (); Kimchi and Vishwanath (2014); Kimchi, Analytis, and Vishwanath (); Lee et al. (2014a); Nasu et al. (2014); Lee et al. (2014b); Kin-Ho Lee et al. (); Hermanns and Trebst () motivated in part by the recent synthesis of three-dimensional honeycomb Iridates Modic et al. (); Takayama et al. ()..

Figure 1: (a) Crystal structure of the layered Iridates AIrO with A = (Na, Li). (b) Sketch of the microscopic interactions in the Kitaev honeycomb model. (c) Phase diagram of the quantum Kitaev model.

In this manuscript, we inspect the role of distortions, i.e. spatial anisotropies of the exchange coupling strength, on the collective spin-orbital state of the hexagonal Heisenberg-Kitaev model away from the exactly solvable Kitaev limit. Our motivation to do so has been twofold. First, early space group determinations of the layered Iridate NaIrO using powder x-ray diffraction scans Singh and Gegenwart (2010) hinted at space group C2/c, in which the hexagonal lattice formed by the Ir ions is slightly distorted along one of its three principal directions. However, more refined inelastic neutron scattering Choi et al. (2012) and single-crystal x-ray diffraction measurements Ye et al. (2012) later revealed that the correct space group of NaIrO is in fact space group C2/m and the hexagonal lattice formed by the Ir ions is an almost perfectly 120 symmetric honeycomb lattice. As we will show in this manuscript the collective spin-orbital states of these systems are nevertheless highly sensitive to small spatial anisotropies of the exchange couplings, which experimentally can be probed via external pressure measurements inducing small lattice distortions and concurrent exchange anisotropies. Second, we hoped to shed further light on the putative quantum critical point in the undistorted Heisenberg-Kitaev model Chaloupka, Jackeli, and Khaliullin (2010); Jiang et al. (2011); Schaffer, Bhattacharjee, and Kim (2012); Reuther, Thomale, and Trebst (2011) between a gapless spin-orbital liquid phase extending out of the Kitaev limit and a conventionally ordered “stripy” phase for the intermediate regime of roughly equally strong Heisenberg and Kitaev couplings. Our analysis shows that exchange coupling distortions are relevant perturbations in any field theoretical description of such a quantum critical point, which depending on their relative strength induce different types of conventionally ordered states in an order-by-disorder selection. This mechanism, which for an infinitesimally small distortion selects a subset of the six possible stripy spin-orbital orderings of the undistorted model, is at play for the entire stripy phase of the Heisenberg-Kitaev model in the intermediate coupling regime. In fact, the selection process turns out to be subtly sensitive not only on the sign of the distortion but also the relative coupling strength of Heisenberg and Kitaev exchange which leads to a total of four different stripy ordered phases in the phase diagram of the distorted Heisenberg-Kitaev model.

We will start our discussion by first considering the classical variant of the distorted Heisenberg-Kitaev model in section II. The phase diagram of the classical model already includes all of the conventionally ordered phases found in its quantum mechanical counterpart as well as its own variation of an order-by-disorder selection of ordered states in the presence of exchange coupling distortions. The entire phase diagram of the classical model as well as its finite-temperature behavior are discussed via extensive numerical simulations. We further consider in detail the classical limit of the Kitaev model, which in the absence of distortions is known to exhibit a classical spin liquid state with Coulomb gas correlations.Chandra, Ramola, and Dhar (2010) We show that the inclusion of exchange distortions leads to a break-down of these power-law correlations and a partial lifting of the residual entropy at zero-temperature, which is also reflected in characteristic signatures of the low-temperature specific heat behavior. We then turn to the quantum Heisenberg-Kitaev model in section III whose phase diagram we have determined via extensive numerical simulations relying on the density matrix renormalization group (DMRG) on finite two-dimensional clusters. The quantum order-by-disorder selection is discussed and found to be in perfect agreement with the numerical data. Finally we discuss the possibility of an exotic continuous quantum phase transition, where the Heisenberg exchange drives the system out of the gapped spin liquid phase of the distorted Kitaev model into a stripy ordered phase. Based on perturbative arguments we conjecture that this transition might be driven by the simultaneous condensation of the excitations of the spin liquid. We round off the manuscript with a summary and outlook in section IV.

Ii Classical Heisenberg-Kitaev model

Figure 2: Low-temperature phase diagram of the Heisenberg-Kitaev model under variation of the relative strength of the Heisenberg and Kitaev couplings and the distortion . For we additionally show histograms of the Néel magnetization , while the histograms for and display the stripy order parameter . In the undistorted case with the magnetization vector lies on one of the cubic axes, either in positive or negative direction, yielding a sixfold degeneracy. In the distorted models with depending on there is either a twofold degeneracy with the magnetization pointing in -direction or a fourfold degeneracy where the magnetization points in one of the - or -directions. The color coding is according to a normalization by highest density. Each histogram has been measured in a single parallel-tempering simulation of a system of size . At the model is symmetric and hence we find no preferred directions of ordering in the -histograms (right hand side).

We start our discussion of the distorted Heisenberg-Kitaev model by first considering its classical version. Its Hamiltonian is given by


where the spins are classical Heisenberg spins and the sums run over nearest neighbor bonds along the three principal directions of the honeycomb lattice labeled ,, and , see Fig. 1(b). The coupling constants parametrize the overall strength of the couplings along these three bonds, while the parameter parametrizes the relative strength of the Heisenberg and Kitaev exchange with corresponding to the Heisenberg limit and corresponding to the Kitaev limit. Note that the Heisenberg exchange is always antiferromagnetic, while the Kitaev exchange is always ferromagnetic. The choice of these coupling signs is motivated by the microscopic modeling Chaloupka, Jackeli, and Khaliullin (2010) of the layered Iridate compounds NaIrO and LiIrO. To be even more explicit, the Hamiltonian can be decomposed into three types of bond terms which read


The case of corresponds to spatially isotropic coupling strengths and the model reflects the rotational symmetry of the honeycomb lattice. We refer to this case as the undistorted Heisenberg-Kitaev model. To consider the effect of distortions, i.e. spatially anisotropic coupling strengths, we will vary the relative strength of the bond exchange while keeping the other two coupling strengths equal, i.e. . We further use the convention that the overall coupling strength is constant, i.e. , so that for varying we have .

ii.1 Phase diagram of the distorted HK model

A summary of the low-temperature ordered states of this classical model is provided in the phase diagram of Fig. 2. The model exhibits a number of conventionally ordered states which we will discuss in the following.

We start by surveying the phases of the undistorted, symmetric model for , see the center row of Fig. 2. At we have an antiferromagnetic Heisenberg interaction stabilizing a Néel ordered phase with a staggered moment pointing along an arbitrary direction. Including a small (ferromagnetic) Kitaev interaction lowers the continuous O(3) symmetry of the Heisenberg model to a set of discrete symmetries including (i) time reversal symmetry, (ii) a spin rotation about the [111] spin axis along with lattice rotations about an arbitrary site, and (iii) an inversion symmetry around any plaquette or bond center. Yet the Néel order survives. Interestingly, the direction of the Néel staggered moment is determined by a classical order-by-disorder mechanism, which we will discuss in more detail in Section II.2. Upon further increasing the Kitaev exchange the system will eventually disfavor Néel order and undergo a first-order transition to an alternate ordered state exhibiting “stripy” order. To see the order of the resulting phase, fortunately, at after an appropriate change of spin variables the Hamiltonian reduces again to an O(3) symmetric model, albeit a ferromagnetic one Chaloupka, Jackeli, and Khaliullin (2010).

We briefly describe the four-sublattice basis transformation. Note that at the spin-spin interactions between , and spin components have equal magnitude but depending on the bond type two interactions are antiferromagnetic and one is ferromagnetic. This interaction can be transformed to a fully ferromagnetic one upon a relative -rotation of the two spins around the special axis. We denote the new spin variables by . Explicitly, to make this transformation on the full lattice we define a 16 site supercell with sites of types as depicted in Fig. 3. The new spin variables are obtained by a -rotation around , , or for sites of type 1, 2, and 3, respectively, and they are simply equal to on sites of type 0.

Figure 3: 16 site supercell used to transform the Hamiltonian at to a O(3) symmetric ferromagnetic one. Note that we use a smaller unit cell compared to Ref. 3.

After the four-sublattice basis transformation the Hamiltonian in the new spin variables reads Chaloupka, Jackeli, and Khaliullin (2010)


Thus we see that at the system has O(3) symmetry. The ground state is a ferromagnet in the variables. This translates to the stripy phases of the original spins; see Fig. 4. Similar to the Heisenberg point at , also at the direction of the ferromagnetic moment is arbitrary due to the O(3) symmetry. But any finite deviation from breaks the continuous symmetry down to a discrete one and we expect the ferromagnetic magnetization direction to be fixed at one of few discrete possibilities. As will be seen in section II.2 this happens by a classical order by disorder mechanism.

Figure 4: Illustration of the stripy , and phases where the arrows indicate the spin alignment along the , , and spin directions.

The Néel and stripy phases have direct analogs in the quantum case. The most interesting quantum phase occurring for , which is a spin liquid with gapless excitations in the form of emergent Majorana fermions, does not have an immediate classical analog. Instead, the system forms a classical spin liquid state – a so-called Coulomb gas, Chandra, Ramola, and Dhar (2010) which exists only in the Kitaev limit, i.e. , to which we will devote special attention in section II.4.

We now consider a finite amount of distortion . corresponds to strong dimers, while corresponds to dominating chains. As can be easily obtained by calculating the energies of the various ordered states discussed, the Néel ordered region splits up into one () in which spins are in the plane and another one () at which they point along the direction. Also in the stripy phases spins either point along for , or they lie in the plane for . Note, that from pure energetics the directions of the spins in the plane is not fixed. Also here the finite temperature order by disorder mechanism come to play; see section II.2.

The paragraph above relies on the following expressions for the energy per unit cell of the Néel and stripy phases with spins pointing along ,


Also the Néel-stripy phase transition lines can be found by equating energies. From


we obtain , giving the line boundary between Néel and stripy for . By comparing


we also obtain , giving the line boundary between Néel and stripy phases for . As a result there is a straight vertical line at marking the Néel-stripy transition in the low-temperature phase diagram of Fig. 2.

As outlined in section II.3 we numerically obtain the finite-temperature phase diagram shown in Fig. 5.

ii.2 Order by disorder and effective Ginzburg-Landau theory

At the magnetization points along an arbitrary direction due to the O(3) symmetry explicitly apparent in Eq. (3) (we refer to the variables in terms of which the Hamiltonian is ferromagnetic). At finite deviations from this symmetric point one expects the Kitaev anisotropic interactions to stabilize a discrete set of orientations of the magnetization. However, as Eq. (II.1) shows, on the mean field level all uniform ferromagnetic states in the O(3) order parameter manifold remain degenerate for . Similarly, the mean field energy in the stripy phases is still invariant under continuous rotations in this plane. Along the same lines, on the mean field level the order parameter in the Néel phase for is not determined.

As we will now see the Heisenberg-Kitaev model provides a simple example where Villain’s order by disorder mechanism comes into play and restricts the order parameter to lie in a subspace of the degenerate manifold. This mechanism requires finite temperatures, where entropic contributions to the free energy become effective. The formal procedure followed below is to integrate out the leading thermal fluctuations, and see that for certain directions of the ordered moment those fluctuations are softer and can further lower the free energy.

We shall consider explicitly the stripy region in terms of the variables. We introduce a slowly varying ferromagnetic order parameter field of unit length


and define gradients along the directions of the three bonds, , where , with unit vectors and . We set the length of these bonds to unity such that the hexagon area is and the area of the Brilloiuin zone is . Expanding the spin-spin interaction Eq. (3) up to second order in gradients we obtain the continuum Hamiltonian , with


For simplicity we focus on the case .

We now consider the partition function of the continuum model Eq. (8),


We proceed by describing the magnetization in terms of fluctuations corresponding to two Goldstone modes and around a uniform magnetization ,


Here , and the set of unit vectors forms an orthonormal basis. This allows to rewrite the partition function as


hence introducing an effective Hamiltonian of by integrating over the fluctuations,


In appendix A we compute explicitly by expanding up to quadratic order in the fluctuations . Up to a constant and up to quadratic order in , we obtain the symmetry allowed anisotropic term


This is the main result of this section. Its negative sign restricts the magnetization in the stripy phase to lie along one of the cubic axes. This term is quadratic in , implying the same conclusion for both sides of the point in the phase diagram at . Similarly, in the stripy phases, by the same argument the magnetization is restricted to either the or cubic axes. On the classical level the Néel ordered phase has an equivalent description as the ferromagnet, and our order by disorder calculation implies that the Néel order parameter is restricted to point along one of the cubic axes.

ii.3 Numerical results

Figure 5: Finite-temperature phase diagrams for (a) , (b) and (c) . We estimate the temperature of transition to the ordered phase by the intersection point of order parameter Binder cumulant plots for and for averaged over multiple pairs of lattice sizes . See Fig. 6 for example data that went into this calculation. The dashed lines at indicate the parametrizations for which the Heisenberg-Kitaev model is symmetric and as a consequence of the Mermin-Wagner theorem is not expected to display finite-temperature transitions in good agreement with our numerical analysis. The degenerate dimer-covering states at also do not undergo any phase transition for .
Figure 6: Binder cumulant curves of the order parameter evaluated over for various system sizes close to their crossing, which gives an estimate of the transition temperature for several example parameter sets. Shown are the cumulant of for the transitions to (a) the stripy phase and (b) the stripy phase as well as (c) the cumulant of for the transition to the Néel phase. Symbols with error bars are single temperature data, while continuous lines are interpolated by multiple histogram reweighting.

Our analysis of the classical Heisenberg-Kitaev model is complemented by an extensive finite-temperature Monte Carlo study. In our simulations the classical spins are situated on the vertices of hexagon-shaped clusters with periodic boundary conditions, which realize the symmetry of the honeycomb lattice and allow to observe unbiasedly all possible orientations in the stripy phases; see Fig. 7. A cluster with a side length of plaquettes contains sites.

We apply the standard Metropolis algorithmNewman and Barkema (1999); Janke (2008) with two different types of proposed moves: In one lattice sweep we first perform local updates of each individual spin, where the new orientation is chosen from an angular region around the old orientation, which has been tuned in such a way during thermalization that acceptance ratios of are maintained at all temperatures. In a second stage we then propose “bond-flip” moves. In one of these moves we choose a random pair of nearest-neighbor sites together with their associated bond-direction . Then for the spins at both sites we reverse the sign of the spin-component linked via that bond in the Kitaev interaction: and , whereas the other components are not modified. While the bond-flip update would not be ergodic on its own, in combination with the single-spin update it greatly accelerates simulation dynamics in the stripy phases, vastly facilitating equilibration.

To further improve ergodicity we combine these canonical updates with a parallel-tempering scheme.Geyer (1991); Hukushima and Nemoto (1996) Here we simulate multiple replicas of the spin system concurrently at different temperatures and exchange configurations between them in a controlled manner that satisfies detailed balance. In this way short autocorrelation times at high temperatures can be exploited to easily overcome free energy barriers at low temperatures, and we can reach all relevant regions of phase space in a single simulation regardless of initial conditions.

We measure two vector order parameters to distinguish between different antiferromagnetic spin-alignments:


Here and stand for the two sublattices of the honeycomb lattice, while the four honeycomb sublattices formed by the sites of the different types of the supercell of Fig. 3 are denoted by , , and . Fig. 7 shows how these sublattices are assigned in our finite lattices. corresponds to perfect Néel order, while is realized for perfect stripy order. The preferred orientations of the magnetization vectors and reflect which ordering directions are possible in the different Néel and stripy phases. In Eq. (II.3) we have chosen an asymmetric definition of the order parameter , where one of the sublattice magnetizations is counted negative and three are counted positive. With this definition is simultaneously an order parameter for the stripy , and phases on the same lattice.222Note that Refs. Price and Perkins, 2012, 2013 use an alternative definition of the stripy-order parameter, which is specified on a choice of sublattices different from ours. By measuring histograms of the components of and we were able to verify the analytical arguments of section II.2. We obtain planar representations of and by mapping the three Cartesian basis vectors to the complex plane as in , and and show the resulting histograms as insets in the phase diagram of Fig. 2. Both the carefully chosen shape of the finite lattices and the parallel-tempering algorithm are essential tools allowing us to fully explore configuration space in our simulations as reflected in these histograms.

Recently Price and Perkins Price and Perkins (2012, 2013) studied the undistorted, -symmetric classical Heisenberg-Kitaev model at finite temperature. Following their analysis we study the Binder cumulants of the absolute valued order parameters


in order to pinpoint the precise temperature of the transitions into the ordered phases. At criticality their values depend only weakly on the system size. Hence the intersection point of or curves evaluated for different gives a good estimate of the critical temperature.

Interestingly, Price and Perkins found that for , the entrance to the ordered phases (Néel or stripy) from the high-temperature paramagnetic phase undergoes two consecutive phase transitions, via a small sliver of a critical Kosterlitz-Thouless phase. In this intermediate phase, the effective model is a 6-state clock model, corresponding to the 6 possible stripy or Néel phases, where an effective U(1) symmetry emerges. However, for the distorted model there are only 2 or 4 degenerate stripy or Néel phases. In this case the intermediate U(1) symmetric phase is not expected José et al. (1977) and no evidence of it is found in our numerical analysis.

We apply standard multiple histogram reweighting techniquesFerrenberg and Swendsen (1989); Chodera et al. (2007) to the temperature-sorted observable time series in combination with numerical minimization routinesPress et al. (2007) to find the intersection points for systems of different sizes up to . Statistical uncertainties are estimated by performing the entire analysis on jackknife resampled data sets.Efron (1982) Plots of Binder cumulants close to their crossing points are given in Fig. 6 for several parameter sets. We average over the results for different values of to estimate the transition temperatures shown in Fig. 5.

Figure 7: example of the finite lattices used in our Monte Carlo simulations, where opposing boundaries are periodic as indicated. The numbers show the division into the sublattices used in the definition of the order parameter in Eq. (II.3), which matches the supercell in Fig. 3.
Figure 8: A zero-temperature configuration of the symmetric Kitaev model is a generic dimer covering state, to which we can associate a divergence-free field.

For the symmetric case our approach resolves the lower of the two transition temperatures of the analysis of Ref. Price and Perkins, 2013. For the distorted model, only a single transition is expected as argued above. We associate this transition with the crossing point of the Binder cumulant curves of different sizes as shown in the examples in Fig. 6.

ii.4 Emergent magnetostatics in the Kitaev limit

Before concluding our discussion of the classical Heisenberg-Kitaev model, we will briefly discuss the physics of the Kitaev limit . While its quantum mechanical counterpart is well known as a paradigmatic, exactly solvable spin model harboring various spin liquid ground states, the classical Kitaev model certainly deserves some attention as well. In its undistorted form it is one of the simplest, analytically tractable classical spin models that evades a thermal phase transition and harbors a classical spin liquid state, which at zero temperature exhibits an extensive degeneracy and pair correlations decaying with a characteristic power-law Chandra, Ramola, and Dhar (2010). These zero-temperature features can be traced back to an effective description in terms of emergent magnetostatics – an example of a so-called Coulomb gas Henley (2010). We will briefly review the arguments showing the origin of this emergent spin liquid in the classical Kitaev model in the following with a more detailed and self-consistent account being given in appendix B. We then discuss the effect of finite distortions, which lead to a (partial) lifting of the zero-temperature degeneracy and a break-down of the Coulomb correlations. However, characteristic remnants of the Coulomb description remain as signatures in the low-temperature specific heat as we detail in the subsequent subsection.

As noted earlier the undistorted classical Kitaev model incorporates a high level of exchange frustration with each spin being subject to competing magnetic exchanges that equally favor alignment along one of the three orthogonal axes of a classical Heisenberg spin. As one approaches the zero-temperature limit of this model it is easy to see Chandra, Ramola, and Dhar (2010) that the total energy of the system can be minimized by spin configurations where spins align in a pairwise fashion along one of the three easy axes of the magnetic exchange, i.e. the one favored by the bond between the two spins forming a pair. An example of such a spin configuration is illustrated in Fig. 8(a). Since every spin is part of precisely one such aligned pair, we can identify each pair of aligned spins with a ‘dimer’. As a consequence, any such energy minimizing spin configuration can be mapped to a hardcore dimer covering of the honeycomb lattice as illustrated in Fig. 8(b) where every site (spin) is part of precisely one dimer. This mapping allows two immediate conclusions. First, it is well known since the early work of Wannier Wannier (1950), Kasteleyn Kasteleyn (1963), and Elser Elser (1984) that the number of dimer coverings on the hexagonal lattice grows exponentially in the system size and as thus we can immediately estimate the zero-temperature degeneracy of the spin model. Second, it has long been appreciated Henley (2010) that the hard-core dimer constraint on a bipartite lattice allows a mapping of any dimer covering to a divergence-free field configuration, which is schematically illustrated in Fig. 8(c). It is precisely this description of the zero-temperature spin configurations in terms of a divergence-free magnetic field that allows to draw the connection to an emergent Coulomb gas description. The latter is well known to give rise to power-law correlations, which translated back to the original spin model are pair correlations of the form

For a detailed and self-consistent description of the Coulomb gas formulation of the zero-temperature classical Kitaev model we refer the reader to appendix B.

When introducing distortions of the exchange couplings the extensive degeneracy of zero-temperature states is immediately lifted. For two spin configurations are singled out where spins align along the -direction again in a pairwise fashion – with both states being mapped to an identical dimer covering as illustrated on the left-hand side in Fig. 9. As a consequence, the spin liquid physics disappears entirely and the system undergoes a conventional symmetry breaking thermal phase transition into one of the two states. For a different picture emerges. While the extensive zero-temperature degeneracy is still lifted, the system retains a subextensive degeneracy down to zero temperature where the spins align in pair-wise fashion along the zig-zag chains spanned by the and -bonds as illustrated on the right-hand side in Fig. 9. The consequence of this lifting again is the loss of Coulomb correlations, but the system still evades a conventional ordering transition down to zero temperature with characteristic features arising for instance in the specific heat as discussed in the next subsection.

Specific heat and zero modes

One characteristic feature of the extensive manifold of zero-temperature spin configurations is that it gives rise to certain soft fluctuations called zero modes. Following the pioneering work of Chalker et al. Chalker, Holdsworth, and Shender (1992), we show in the remainder of this section that these zero modes reduce the specific heat in its limit in a universal way – a characteristic signature that as we show can easily be tracked by numerical simulations of the classical spin model.

To start our discussion of the analytical arguments we consider fluctuations around a given dimer covering or spin configuration, respectively. Each spin belonging to a dimer on a -bond gives rise to possible fluctuations in the two directions orthogonal to . For example for a spin belonging to a dimer and pointing along we write


The fluctuations in the and directions influence also the component due to the unit constraint .

Figure 9: Left: The preferable dimer covering state for dominating . Right: Typical dimer covering states for

Let denote the set of dimerized bonds. For , and assuming for simplicity that this is a type bond, the Kitaev spin-spin interaction reads


We see that up to quadratic order, fluctuations do not interact across dimerized bonds (no coupling terms for ). On the other hand, for a non-dimerized bond [see Fig. 8(b)] the Kitaev interaction reads


Thus, expanding the Hamiltonian in to quadratic order, the fluctuation corrections consist of decoupled terms which live on the non-dimerized bonds and read



Figure 10: Low temperature behavior of the specific heat per spin in the classical Kitaev model with different distortions. Shown are Monte Carlo results obtained at temperatures demonstrating that in the limit of one finds , and . The data has been obtained for systems of side length .

Interestingly, for ,


This implies the existence of a zero mode: does not appear in . This zero mode has been identified Baskaran, Sen, and Shankar (2008) to be a sliding degree of freedom of the dimer covering states. For low enough temperatures fluctuations become small and the partition function becomes


For any quadratic eigenmode , with energy , the contribution to the specific heat then becomes


independent of the coefficient . However, in our system we have to further consider the contributions of the zero modes. For those modes we need to go to quartic order, i.e. , for which the contribution to the specific heat can be estimated to be


again independent of the coefficient . In a standard state without zero modes (such as a ferromagnetic state) we would have two quadratic modes ( and ) per spin. This would give the zero temperature value of the specific heat per spin


However, in the Coulomb phase of the classical Kitaev model, we have only one zero mode for each quadratic mode, hence


We now consider the effect of a finite distortion, i.e. , which splits the degeneracy of the various dimer covering states. For , namely , the dimer covering with only -dimers has the lowest energy; see Fig. 9.

At the same time, fluctuations around this state are described by Eq. (19), which can be written as


For the two coefficients in this equation are positive, leaving no zero modes. Hence


For the dimers cover or bonds in the ground state; see Fig. 9. Now consider fluctuations around these 1D covering states. The Hamiltonian for the fluctuations is the same as Eq. (19), but now there are two types of non-dimerized bonds. For with or , has the form of Eq. (21), implying a zero mode. But for with , the Hamiltonian has the form of Eq. (27), implying no zero mode. As a result the specific heat per spin becomes


Our Monte Carlo calculations, summarized in Fig. 10, nicely reproduce these fractions and are thus able to pinpoint the different constraints on the dimer covering states underlying the Coulomb gas.

Iii Quantum Heisenberg-Kitaev model

Figure 11: Phase diagram of the distorted quantum Heisenberg-Kitaev model.

We now turn to a discussion of the quantum version of the distorted Heisenberg-Kitaev model, i.e. we again consider the Hamiltonian


where the spins are now quantum mechanical spin-1/2 degrees of freedom. (In our convention are represented by Pauli matrices .) The exchange parameter again interpolates between the antiferromagnetic Heisenberg model and the ferromagnetic Kitaev model and the distortion of the exchange couplings is parametrized by with the simultaneous conditions that all three spin exchange couplings add up to a constant, i.e. , and . The case then corresponds to the undistorted situation where the spin exchange along all three bonds has equal magnitude, i.e. . The limit corresponds to decoupled dimers on the -bonds, while the opposite limit of corresponds to decoupled zig-zag chains along the - and -bonds.

When exploring the -parameter space we find that the above model not only harbors quantum analogues of all classically ordered states, but exhibits a number of additional genuinely quantum states including a valence-bond solid and two spin-orbital liquid phases, which both extend well beyond the well-studied Kitaev limit of the quantum model. In fact, one of the more interesting features of the extended phase diagram of the quantum Heisenberg-Kitaev model is the possible occurrence of unconventional continuous phase transitions between these gapped and gapless spin-orbital liquid phases and conventionally ordered states.

In the following, we will first discuss the general quantum phase diagram of the distorted Heisenberg-Kitaev model and the numerical simulations underlying its determination and then focus our discussion on the possibly interesting quantum critical behavior associated with the phase transition out of one of the spin-orbital liquid phases.

iii.1 Phase diagram of the quantum model

The phase diagram of the quantum Heisenberg-Kitaev model in the presence of exchange distortions is summarized in Fig. 11. Similar to the classical model we find an extended Néel ordered phase around the Heisenberg limit which upon distorting the exchange interactions undergoes a quantum order-by-disorder transition locking the spin orientation in the ordered phases to the ( or ) direction for (), respectively. For the system undergoes a transition into a valence bond solid (VBS), which adiabatically connects to the limit of isolated dimer singlets on the -bonds in the limit (and ).

For the quantum model exhibits an symmetry that is again rooted in the observation that for this ratio of the Heisenberg and Kitaev couplings the model can be mapped via the four-sublattice basis transformation illustrated in Fig. 3 to a ferromagnetic Heisenberg model. In fact, such a mapping exists for all values of the distortion , i.e. the quantum model exhibits an entire symmetric line for . In the four-sublattice rotated basis the ground state of the quantum model is a simple ferromagnet for , which transformed back into the original basis becomes a ‘stripy ferromagnet’ akin to the illustrations in Fig. 4. In the undistorted case the ground state is sixfold degenerate with the six possible stripy states of Fig. 4 having equal weight in the ground state. This picture changes immediately upon moving away from the line and distorting the exchange couplings. Again a quantum order-by-disorder transition (detailed in appendix C) selects a subset of these six stripy states with four different phases emerging around the undistorted point in the middle of our phase diagram in Fig. 11. In complete analogy to the classical model, a subset of two stripy FM states locking the spins into the -direction is selected for as well as for . For the other two quadrants and the opposite subset of four stripy FM states with the spins locking into either the or -directions are selected by the quantum order-by-disorder mechanism, see appendix C for details.

Arguably the most interesting phases in our phase diagram are the two spin liquid phases emerging for dominating Kitaev couplings. For the undistorted Heisenberg-Kitaev model it was previously established Chaloupka, Jackeli, and Khaliullin (2010); Jiang et al. (2011) that the stripy FM phase gives way to a gapless spin liquid phase for , i.e. Kitaev couplings which are about 8 times larger than the isotropic Heisenberg exchange. This gapless spin liquid phase remains stable when introducing an exchange distortion and is found to occupy a rather extended regime in the -parameter space as illustrated in Fig. 11. For the pure Kitaev model it is well known Kitaev (2006) that the gapless spin liquid can be gapped out into a topological spin liquid if one introduces an exchange distortion that renders one of the three coupling exchanges dominant, i.e. in our notation, see Fig. 1(c). Upon including a Heisenberg exchange this gapped phase must remain stable for a finite parameter regime – however, since the gap itself is rather small the regime occupied by this topological spin liquid in our -parameter space reduces to a small sliver as illustrated in Fig. 11. We come back to a more detailed discussion of the emergence of this topological phase as well as the nature of the quantum phase transition out of this phase into the stripy phase in the next subsection.

Our approach to map out the phase diagram of the quantum Heisenberg-Kitaev model as discussed above is based on various numerical techniques, in particular exact diagonalization (ED) studies and density-matrix renormalization group (DMRG)White (1992) calculations for small, but highly symmetric clusters with up to sites for the DMRG (ED) calculations, respectively. In order to minimize finite-size effects we employed periodic boundary conditions and chose the clusters such that they preserve the SU(2) symmetry of the four-sublattice basis transformation introduced in Section II.1. We used clusters of and sites. For the DMRG calculations we typically kept up to states in the DMRG block and performed multiple sweeps to converge the observables with the typical truncation error becoming of the order of or smaller. The location of the phase boundaries in the phase diagram (see Fig. 11) are determined by the peak position of the second derivatives of the ground state energy density, i.e. and . A similar approach has previously been used to successfully map out the phase diagram of the (undistorted) Heisenberg-Kitaev model in a magnetic field Jiang et al. (2011). Data for these derivatives along representative cuts in the -parameter space are shown in Figs. 12 and 13. A very sharp peak in the second derivative – corresponding to a jump of the first derivative of the ground state energy density, i.e., and – is taken as a signature for a first-order transition and marked by the red solid lines in the phase diagram, while a relative shallow peak in the second derivative data is interpreted as possibly indicating continuous phase transitions.

Figure 12: (Color online) Second derivative of the ground state energy density as a function of for two different values of the distortion .
Figure 13: (Color online) Second derivative of the ground state energy density as a function of the distortion for two different values of the coupling .

To further identify the nature of different phases and compare with the classical Heisenberg-Kitaev model, we calculate a ‘bond magnetization’, i.e. the expectation value of the bond operator where denotes the component of spin, and is the unit vector along an -bond. As illustrated in Fig. 14 this bond magnetization is a very useful tool to track the quantum order-by-disorder selection in the distorted stripy phases. For example, in the stripy- phase for and , the -bond magnetization is positive since points in the same direction in -bond, while and are negative because are antiparallel along the and bonds (not shown). In addition to the stripy phase, this bond operator can also be used to study the phase transition between different phases, which will increase or decrease rapidly across the phase boundary. As an example, we plot the bond operator with in Fig. 14(b), in which the dotted lines are the phase boundaries determined by , and consistent with the ones determined by the second derivative of the ground state energy density.

Figure 14: (Color online) Bond magnetization as a function of for two different values of the distortion . Here .

Finally, we want to shortly comment on the quantum order-by-disorder mechanism playing out in the distorted stripy phase. As mentioned earlier, for precisely the system exhibits an additional symmetry and its ground state can be characterized by a conventional ferromagnetic order parameter in terms of the four-sublattice transformed spin variables introduced in Section II.1. For small deviations from the symmetry of the model is reduced to a discrete one. But as we saw in Section II.1 when discussing the classical model, one can quickly see that on the mean field level the actual direction of the ferromagnetic order is not fixed upon introducing a distortion. In fact, as we have shown in Section II.1 thermal fluctuations are ultimately responsible for the eventual ordering along cubic axes in the classical model. An analogous argument applies to the quantum model where quantum fluctuations will favor a locking of the spin orientation along the cubic axes of the model for finite distortions at zero temperature. This effect has been first commented on in Refs. Khaliullin, 2001 and Chaloupka, Jackeli, and Khaliullin, 2010 and is discussed in detail in appendix C.

iii.2 Phase transition out of the Abelian topological phase

For , , the system decouples into -dimers with Hamiltonian given in Eq. (II). The Heisenberg term has the singlet state as the ground state with an excited triplet , whereas the Kitaev ferromagnetic term has a degenerate pair of ground states and a second degenerate pair of excited states . The energies of these states are , , and . For , one can formulate an effective interaction between the doublet degrees of freedom localized on links and represented by effective spins . Thus for a given link for the state and for the state . Following Kitaev, Kitaev (2006) those spins can be located on the links of a square lattice; see Fig. 15. For small , the dimer-dimer interaction can be represented as an effective interaction between the ’s. For the Kitaev model () the leading interaction is generated at forth order in , and is a 4-spin interaction equivalent to the toric code model. Explicitly for Kitaev (2006)


with , and the plaquette operator , where runs over all hexagonal plaquettes of the honeycomb lattice, which become either plaquettes or stars on the square lattice of the toric code.

Figure 15: Mapping from the honeycomb lattice to the toric code lattice.

In the presence of the Heisenberg term, we find an interaction already at first order in which reads simply


with . Here are nearest neighbors in the square lattice of the toric code. One immediately sees that this term stabilizes a Néel order of the effective spins , which is equivalent to the stripy- phase in Fig. 4. Therefore the phase transition between the topological phase and stripy phase emanates from the right-top corner of the phase diagram. By comparing the energy scales of the interaction , in Eq. (31), stabilizing the topological phase and the Ising interaction , in Eq. (32), one immediately sees that the transition line approaches the right-top point as


This high power of () is consistent with the very small area occupied by the gapped topological phase in our phase diagram in Fig. 11.

Possibility of condensation of excitations

We propose a simple model to understand the quantum phase transition between the gapped topological phase and stripy phase. This model contains just the two competing interactions which stabilize either phase: the toric code Hamiltonian Eq. (31) and the Ising Hamiltonian Eq. (32). In order to introduce standard notation for the toric code model, after permuting the spin indices , and then performing a rotation along for spins living on vertical bonds in the toric code square lattice defined in Fig. 15, the model becomes


where the coupling have been separated into star and plaquette operators with couplings and , shown in Fig. 16. In our case . In the last term always belongs to an horizontal bond and to a vertical bond and are nearest neighbors; see Fig. 16.

Figure 16: represents star operators in Eq. (34), plaquette operators, and Ising coupling on nearest neighbors.

As a function of there must be a quantum phase transition between the gapped topological phase to the Ising ordered phase at with spontaneously broken local Ising symmetry


One can write this model in terms of the excitations of the gapped topological phase: (i) electric excitations living on stars with eigenvalue of


and (ii) magnetic excitations living on plaquettes with eigenvalue of


In the physical Hilbert space both and excitations occur in pairs. Such pairs are created, respectively, by


where () is an arbitrary path along the lattice (dual lattice) connecting stars , (plaquettes , ) where the two excitations are created. One can check that for or , and otherwise, and the ’s satisfy similar relations. Independent of the choice of contours, and commute if the corresponding contours cross an even number of times and anticommute otherwise.

The Hamiltonian is simply


Here each horizontal edge is shared by two stars and and each vertical edge shares two plaquettes and . The reference star and reference plaquette are arbitrary, and can be thought of as being located at infinity (with open boundary conditions).

Clearly in the Néel phase there is a finite expectation value of




and their relative sign is opposite for . In the topological phase all excitations are gapped and uncorrelated. Thus a natural question is how the - and - excitations condense. Typically, excitations condense at a phase transition as their kinetic energy exceeds the mass gap. From the effective model Eq. (39), we see that to first order in individual - and - excitations can not hop thus their excitation energy is and , respectively. On the other hand their bound state does acquire kinetic energy of order . It can hop along the direction hence lowering the gap to . This suggests an interesting type of quantum phase transition consisting of a condensation of the bound states for large enough , which is unusual due to the fermionic nature of those composite particles. It is interesting to explore this possibility on a quantitative level in the future.

iii.3 One-dimensional limit of the Heisenberg-Kitaev model

In the limit of , corresponding to the bottom in the phase diagram in Fig. 11, the system decomposes into decoupled Heisenberg-Kitaev chains. The physics of such chains has previously been partially explored, in particular with regard to its energy dynamics Steinigeweg and Brenig (). Here we will apply one dimensional (1D) field theoretical methods to analytically construct the 1D phase diagram of such Heisenberg-Kitaev chains, and to gain insight into the 2D case by studying the limit of weakly coupled chains.

iii.3.1 Phase diagram

Our phase diagram of the 1D Heisenberg-Kitaev (HK) model is shown in Fig. 17(a). It contains three exactly solvable points: (i) For the model is the antiferromagnetic Heisenberg chain which is described by a conformal field theory (CFT) with central charge . (ii) At , the model written in terms of the spin variables is the ferromagnetic Heisenberg chain, which has dynamical critical exponent . (iii) At , the system is also critical and can be described by a CFT with corresponding to gapless Majorana chains Kitaev (2006). Below we describe the phases in between these three exactly solvable points.

Figure 17: (a) Phase diagram of the 1D Heisenberg Kitaev model. (b) Phase diagram of the XXZ model.

It is convenient to express the 1D HK Hamiltonian as the sum of the well studied XXZ model,


and a perturbation. Indeed, our model reads




The signs correspond to alternating chains. The well known phase diagram of the XXZ model is summarized in Fig. 17(b).

We begin by analyzing the small limit. At the perturbation to the XXZ chain vanishes and . Our system lies inside the gapless Luttinger liquid phase of the XXZ model which extends in the range . This phase is described by a Luttinger liquid theory Giamarchi (2004) which is characterized by Luttinger parameter and velocity , given exactly by