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

The microscopic modeling of spin-orbit entangled $j=1/2$ Mott insulators such as the layered hexagonal iridates ${\mathrm{Na}}_2{\mathrm{IrO}}_3$ and ${\mathrm{Li}}_2{\mathrm{IrO}}_3$ 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 that include not only conventionally ordered phases but also gapped and gapless spin-orbital liquids. In particular, we demonstrate that in the presence of distortions, i.e., 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_2$ 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_2$ phase quickly breaks down in what might be a rather unconventional phase transition driven by the simultaneous condensation of its elementary excitations.

Figure 1

(a) Crystal structure of the layered iridates $A_2{\mathrm{IrO}}_3$ with $A=\text{Na}$, Li. (b) Sketch of the microscopic interactions in the Kitaev honeycomb model. (c) Phase diagram of the quantum Kitaev model.

Figure 2

Low-temperature phase diagram of the Heisenberg-Kitaev model under variation of the relative strength $\alpha$ of the Heisenberg and Kitaev couplings and the distortion $J_z$. For $0<\alpha <1/3$, we additionally show histograms of the Néel magnetization ${\mathbf{m}}_{\text{N}}$, while the histograms for $1/3<\alpha <1/2$ and $1/2<\alpha <1$ display the stripy order parameter ${\mathbf{m}}_{\text{S}}$. In the undistorted case with $J_z=1$, 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 $J_z\gtrless 1$ depending on $\alpha$, there is either a twofold degeneracy with the magnetization pointing in $\pm z$ directions or a fourfold degeneracy where the magnetization points in one of the $\pm x$ or $\pm y$ 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 $L=32$. At $\alpha =1/2$ the model is O(3) symmetric and hence we find no preferred directions of ordering in the ${\mathbf{m}}_{\text{S}}$-histograms (right-hand side).

Figure 3

16 site supercell used to transform the Hamiltonian at $\alpha =1/2$ to a O(3) symmetric ferromagnetic one. Note that we use a smaller unit cell compared to Ref. [(3)].

Figure 4

Illustration of the stripy $x$, $y$, and $z$ phases where the arrows indicate the spin alignment along the $x$, $y$, and $z$ spin directions.

Figure 5

Finite-temperature phase diagrams for (a) $J_z>1$, (b) $J_z=1$, and (c) $J_z<1$. We estimate the temperature of transition to the ordered phase by the intersection point of order parameter Binder cumulant plots $Q_2^{\text{N}}\left(T\right)$ for $\alpha <1/3$ and $Q_2^{\text{S}}\left(T\right)$ for $\alpha >1/3$ averaged over multiple pairs of lattice sizes $L$. See Fig. 7, for example, data that went into this calculation. The dashed lines at $\alpha =0,1/2$ indicate the parameterizations for which the Heisenberg-Kitaev model is O(3) 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 $\alpha =1$ also do not undergo any phase transition for $T>0$.

Figure 6

$L=4$ 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 ${\mathbf{m}}_{\text{S}}$ in Eq. (14), which matches the supercell in Fig. 3.

Figure 7

Binder cumulant curves of the order parameter evaluated over $T$ for various system sizes $L$ close to their crossing, which gives an estimate of the transition temperature for several example parameter sets. Shown are the cumulant of ${\mathbf{m}}_{\text{S}}$ for the transitions to (a) the stripy $z$ and (b) the stripy $xyz$ phases as well as (c) the cumulant of ${\mathbf{m}}_{\text{N}}$ for the transition to the Néel $z$ phase. Symbols with error bars are single temperature data, while continuous lines are interpolated by multiple histogram reweighting.

Figure 8

A zero-temperature configuration of the $C_3$ symmetric Kitaev model is a generic dimer covering state, to which we can associate a divergence-free field.

Figure 9

(Left) The preferable dimer covering state for dominating $J_z$. (Right) Typical dimer covering states for $J_x=J_y>J_z$.

Figure 10

Low-temperature behavior of the specific heat per spin $C_v\left(T\right)$ in the classical Kitaev model with different distortions. Shown are Monte Carlo results obtained at temperatures $T\ge 1/2000$ demonstrating that in the limit of $T\to 0$, one finds $C_v^{J_z=1}\to 3/4$, $C_v^{J_z>1}\to 1$ and $C_v^{J_z<1}\to 7/8$. The data have been obtained for systems of side length $L=16$.

Figure 11

Ground-state phase diagram of the quantum Heisenberg-Kitaev model with regard to the relative strength $\alpha$ of the Heisenberg and Kitaev couplings and the distortion $J_z$. The line $J_z=1$ corresponds to the undistorted model, which for increasing coupling parameter $\alpha$ shows a sequence of a Néel ordered antiferromagnet, a stripy-ordered ferromagnet, and a gapless spin liquid around the Kitaev limit. For a finite distortion, the magnetically ordered states undergo a quantum order-by-disorder transition, which locks the spin orientation in the ordered phases to the $z$ ($x$ or $y$) direction as indicated in the diagram. The quantum phase transitions between these various ordered states are all first-order as indicated by the red solid lines in the phase diagram. For large distortions $J_z\gtrsim 1.35$, the Néel antiferromagnet is destabilized in favor of a valence-bond solid (VBS) state. Arguably, the most interesting phases in this phase diagram are an extended gapless spin liquid around the undistorted Kitaev limit as well as a topologically ordered spin liquid for $J_z>3/2$ around the Kitaev limit. The possibly continuous phase transitions out of the gapped spin liquid phase is discussed in the main text.

Figure 12

Second derivative of the ground-state energy density as a function of $\alpha$ for two different values of the distortion $J_z$.

Figure 13

Ground-state energy density as a function of the distortion $J_z$ for two different values of the coupling $\alpha$. The kink at $J_z=1$ clearly indicates a first-order transition.

Figure 14

Bond magnetization $B_a^{\gamma }$ as a function of $\alpha$ for two different values of the distortion $J_z$. Here, $\gamma =a=x,y,z$.

Figure 15

Mapping from the honeycomb lattice to the toric code lattice.

Figure 16

$A$ represents star operators in Eq. (34), $B$ plaquette operators, and $J$ Ising coupling on nearest neighbors.

Figure 17

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

Figure 18

Entanglement scaling for the 1D Heisenberg-Kitaev chain for various coupling strength $\alpha$. (a) Data for a periodic $L=100$ site chain in the gapped phase for $\alpha <1/2$. (b) Data for a periodic $L=100$ site chain in the gapless phase for $\alpha >1/2$. (c) Data for an open $L=256$ site chain in the gapless phase for $\alpha >1/2$ and $\alpha =1$. Boundary effects of the open chain result in an odd-even staggering.

Figure 19

Quantum fluctuation contribution to the ground-state energy per unit cell as a function of azimuthal angle along the equator as parameterizing the magnetization $\widehat{\mathbf{e}}$, with $\alpha =0.45$ $J_z=1.5$, or alternatively for $\alpha =0.6$ $J_z=0.5$ showing minima along cubic axes. Numerical evaluation of the $k$ sum in Eq. (C6) has been done for a honeycomb lattice with $10\times 10$ unit cells with periodic boundary conditions.