Ordered states in the Kitaev-Heisenberg model: From 1D chains to 2D honeycomb.

We study the ground state of the 1D Kitaev-Heisenberg (KH) model using the density-matrix renormalization group and Lanczos exact diagonalization methods. We obtain a rich ground-state phase diagram as a function of the ratio between Heisenberg (J = cosϕ) and Kitaev (K = sinϕ) interactions. Depending on the ratio, the system exhibits four long-range ordered states: ferromagnetic-z, ferromagnetic-xy, staggered-xy, Néel-z, and two liquid states: Tomonaga-Luttinger liquid and spiral-xy. The two Kitaev points [Formula: see text] and [Formula: see text] are singular. The ϕ-dependent phase diagram is similar to that for the 2D honeycomb-lattice KH model. Remarkably, all the ordered states of the honeycomb-lattice KH model can be interpreted in terms of the coupled KH chains. We also discuss the magnetic structure of the K-intercalated RuCl3, a potential Kitaev material, in the framework of the 1D KH model. Furthermore, we demonstrate that the low-lying excitations of the 1D KH Hamiltonian can be explained within the combination of the known six-vertex model and spin-wave theory.

Introduction

It is well-known that magnetic ordering (e.g. the Néel state) is induced by the spontaneous breaking of some sort of symmetry at low temperatures. If one can prevent such a conventional magnetic ordering, an exotic type of quantum-disordered state — the so-called quantum spin liquid (QSL) — may arise. In the spin liquid phase, a macroscopic number of quasi-degenerate low-energy states compete with each other. Since Anderson’s proposal of a resonating valence bond state as a theoretical realization of a QSL in 1973[1], geometrical frustration[2], a situation where the network of magnetic interactions is incompatible with the spatial symmetry, had been considered to be almost the key to achieve such situation[3]. But in 2006, the advent of the Kitaev model[4] shifted the focus on another way to realize a QSL — frustration due to bond-dependent exchange interactions[5]. This model has nontrivial topological phases with elementary excitations exhibiting a Majorana algebra. After that, the microscopic origin of such Kitaev-type interactions in the d5 transition metal compounds was worked out[6,7]. In the last decade, there has been a growing number of researches on the Kitaev materials[8,9].

The first Kitaev material candidates are 5d5 honeycomb iridates A2IrO3 (A = Na or Li)[7,10]. The Ir4+ ions centered in the edge-sharing IrO6 octahedra form a planar hexagonal network. In each Ir4+ ion, the two e levels are split off by the crystal field and five electrons in the 5d shell are put into the t2 orbitals with a finite effective angular momentum l = 1. Besides, the strong spin-orbit coupling (SOC) for 5d electrons splits up the manifold into a fully occupied quartet band and a half-filled doublet band. The latter half-filled band gives an effective pseudospin model, in association with the opening of a Mott gap. In the edge-sharing IrO6 octahedra, there are two Ir-O-Ir exchange paths with a ~90° bond angle. The anisotropic Kitaev-type coupling then arises from the particular superexchange interactions between the pseudospins. A more realistic spin model to describe the magnetic properties is the Kitaev-Heisenberg (KH) model, that accounts for the residual Heisenberg-type couplings. It is known that the Kitaev QSL is most likely preempted by a certain level of the Heisenberg interaction. In fact, the above iridates exhibit long-range magnetic order at low temperatures[11-14], possibly assisted by off-diagonal and/or long-range interactions.

At present, one of the most promising materials close to the Kitaev QSL is the 4d5 ruthenium trichloride α-RuCl3 in its honeycomb crystal phase[15-21]. For this compound the effective description, as a Ru3+ 4d analogue to the iridates, seems to still give a good description of the magnetic properties, although some degree of admixture between and in the states may arise from a weaker SOC for 4d electrons than in the Ir 5d orbitals. Neutron and Raman scattering studies gave evidence for fractionalized excitations typical of the Kitaev QSL[22-24], and both theoretical[25] and experimental investigations[26-31] indicated, for this material, the existence of a transition into a QSL in the presence of external magnetic field. Furthermore, a recent X-ray diffraction measurement suggested the Kitaev QSL is induced by a small applied pressure of about 0.7 GPa[32]. A possibility of pressure-induced QSL has been also reported for the other Kitaev materials, honeycomb iridates β-Li2IrO3[33] and γ-Li2IrO3[34]. In both materials the low temperature phase is magnetically ordered without field and under ambient pressure. Namely, the Kitaev QSL is derived through a melting of the magnetic ordering by field or pressure. Therefore, it is crucial to deeply understand the characteristics of the magnetically ordered states.

Very recently, a temperature-dependent electron energy loss spectroscopy (EELS) measurement was performed in a K-intercalated RuCl3, denoted as K0.5RuCl3 (Knupfer M., private communications). The intercalated K+ ions provide charge carriers, however, a sharp gap was observed at ~0.9 eV, instead of the charge gap E = 1.1–1.2 eV for the undoped α-RuCl3[35-37]. This indicates an insulating feature of K0.5RuCl3 and differs from the pseudogap behavior seen for charge localization in disordered metals. This has been interpreted as half of the pseudospins being replaced by nonmagnetic d6 ions. Therefore, this insulating state can be pictured as a formation of superlattice by charge disproportionation (charge ordering). Different possible charge ordering patterns are shown in Fig. 1. Of particular interest in the present context is the zigzag chain-type pattern exhibited in Fig. 1(c), where chains of nonmagnetic ions are separated by chains of magnetic ions with KH interactions.

Figure 1

Possible charge ordering patterns of K0.5RuCl3; (a) triangular, (b) dimer and (c) zigzag. (d) Lattice structure of the 1D Heisenberg-Kitaev model. The labels ‘x’ and ‘y’ indicate the x-bond and y-bond, respectively (see text).

The above motivations lead us to studying the 1D KH model not only for fundamental theoretical reasons, but also as a possible minimal spin model to describe the magnetic properties of the K-intercalated RuCl3. So far, there are few studies on the 1D KH model[38-46]. Using the density-matrix renormalization group (DMRG) method, we calculate total spin, real-space spin-spin correlation functions, static spin structure factor, central charge, and various order parameters in the ground state. Based on the results, we obtain the ground-state phase diagram, including four long-range ordered and two liquid phases, as a function of the ratio between Heisenberg and Kitaev interactions. Moreover, the relevance of the phase diagram to that of the honeycomb-lattice KH model is discussed. It is striking that all the magnetically ordered states of the honeycomb-lattice KH model can be interpreted in terms of coupled 1D KH chains. Furthermore, we calculate the dynamical spin structure factors via the Lanczos exact diagonalization (ED) technique. The basic low-lying excitations are considered by use of the known six-vertex model and the spin-wave theory (SWT). The present investigation can thus contribute to an elucidation of the fundamental properties of the K-intercalated RuCl3 as well as a better insight in the understanding of the physics of the honeycomb-lattice KH model.

Model and Numerical Methods

Pattern of charge ordering in the K-intercalated RuCl3

The insulating feature of K0.5RuCl3 with a sharp peak at ~0.9 eV in the EELS spectrum can be explained by charge ordering associated with a superlattice formation rather than a random distribution of K+ ions (Knupfer M., private communications). The superlattice structures are simply constructed by replacing half of the pseudospins with nonmagnetic sites in the original hexagonal cluster. In practice, there are three possibilities for the charge ordering pattern, they shown in Fig. 1(a–c). Recently, for exfoliated α-RuCl3, a first-order structural phase transition at ~150 K between low-temperature C2/m and high-temperature P31 structures was reported[47]. At high temperature, the emergence of charge ordering, originating from anisotropy in the charge distribution along Ru-Cl-Ru hopping pathways, was observed. It suggests a strong tendency to the bond-directional anisotropy along the Ru-Ru axes. In this sense a realization of the zigzag-type charge ordering [Fig. 1(c)] may be most likely. Besides, the zigzag-type charge ordering can gain the largest exchange energy (see below). To fix the detailed charge distribution, experimental observation and/or microscopic analysis including structural distortion by the K-intercalation are required in future.

1D Kitaev-Heisenberg Hamiltonian

At present, it is widely believed that the magnetic properties of the undoped α-RuCl3 are well described by the KH model on a honeycomb lattice. If we assume the zigzag-type charge ordering in Fig. 1(c), the zigzag chains are well separated by the nonmagnetic ions. Then, each chain is considered to be a 1D KH model, which is equivalent to a system obtained by removing the z-bonds from the honeycomb-lattice KH model. The Hamiltonian of 1D KH model is written asfor a system with L sites, where is a spin - operator at site i, and the Kitaev and exchange couplings are defined as K = sinϕ and J = cosϕ via a phase parameter ϕ. Throughout the paper, we take as the energy unit. The system has two kinds of neighboring links and they appear alternately along the chain. Hence, the structural unit cell contains two lattice sites. Hereafter, we call the links (2i − 1, 2i) and (i, 2i + 1) “x-link” and “y-link”, respectively. By rewriting Eq. (1) aswe can easily notice that a XXZ Heisenberg chain containing exchange () and Ising () terms is disturbed by sign-alternating double-spin-flip () fluctuations[46].

Numerical methods

We employ the density-matrix renormalization group (DMRG) method, which is one of the most powerful numerical techniques for studying various quasi 1D quantum systems[48]. Open boundary conditions are applied unless stated otherwise. This enables us to calculate ground-state and low-lying excited-state energies, as well as static quantities, quite accurately for very large finite-size systems. We are thus allowed to carry out an accurate finite-size-scaling analysis to obtain energies and quantities in the thermodynamic limit L → ∞. We hence study chains with several lengths up to L = 200 sites for a given ϕ. Since, in hindsight, the system (1) exhibits only commensurate phases in the ground state and the largest magnetic unit cell contains four lattice sites, its size is taken as L = 4n (n: integer). For each calculation, we keep up to m = 1200 density-matrix eigenstates in the renormalization procedure and extrapolate the calculated quantities to the limit m → ∞ if needed. Since the SU(2) symmetry is broken in system (1) and total S is no longer a good quantum number except at ϕ = 0, π, one may have some difficulty in obtaining accurate results in comparison to usual DMRG calculations. Nevertheless, in this way, the maximum truncation error, i.e., the discarded weight, is less than 1 × 10−10 while the maximum error in the ground-state is less than 1 × 10−8.

For the dynamical properties calculation, we used the Lanczos exact diagonalization (ED) method. To examine the low-energy excitations for each phase, we calculate the dynamical spin structure factor, defined aswhere γ is z or −(+), and E are the v-th eingenstate and the eigenenergy of the system, respectively (v = 0 corresponds to the ground state). Under periodic boundary conditions, the spin operators can be precisely defined bywhere r is the position of site i and the sum runs over either i even or i odd sites. They provide the same results. The momentum is taken as () since the lattice unit cell includes two sites and the number of unit cells is in a system with L sites. We calculate both spectral functions S±(q, ω) and S(q, ω) as they are different due to broken SU(2) symmetry except at ϕ = 0 and π. We study chains with L = 24, namely, 12 unit cells, by the Lanczos ED method. As shown below, system (1) contains only commensurate phases with unit cell containing one, two, or four sites. Therefore, a quantitative discussion of the low-lying excitations is possible even within the L = 24 chain.

Results

Ground-state properties

Quantum phase transitions

Lets first look at the ground-state energy and total spin with respect to ϕ in order to capture the overall appearance of quantum phase transitions. The results, using a periodic 24-site KH chain, are shown in Fig. 2. We clearly see discontinuities in the first derivative of the ground-state energy at four values: , π, , and , which indicate the first-order phase transitions. The second derivative is, in a precise sense, continuous except for the above four ϕ points, nonetheless, there exists a distinguishable peak around , which may corresponds to a second-order (or continuous) phase transition. Furthermore, as shown below, we find another phase transitions at ϕ = 0. Therefore, we suggest that the simple 1D model (1) exhibits a variety of phases including six quantum phase transitions. It will be confirmed by studying various corresponding order parameters or spin-spin correlation functions. We also confirm that the ground-state energy of the 1D KH chain is always lower than that given by the dimer-type charge ordering [see Fig. 2(a)].

Figure 2

(a) Ground-state energy E0, (b) the second derivative of E0 with respect to ϕ, and (c) total spin S as a function of ϕ, obtaind with a 24-site periodic Kitaev-Heisenberg chain. Dotted line in (a) indicates ground-state energy for the dimer-type charge ordering.

Ferromagnetic-xy phase

Since both K and J are ferromagnetic (FM), a long-range FM ordered state is naively expected in the range . At the spin isotropic point ϕ = π, the spins can align along any arbitrary spatial direction due to SU(2) spin rotation invariance and the total spin takes the maximum value [see Fig. 2(c)]. Away from the isotropic point, the SU(2) symmetry is broken to U(1) and the configurations for higher sectors are projected out [see Fig. 3(c)]. When ϕ is still close to π (), the ground state is approximately expressed aswhere is a basis in real space, namely each site state is represented by either spin-up (↑) or spin-down (↓). The basis is here restricted to a subspace, so that m is summed over all possible combinations of spin configuration with up and down spins in a L lattice sites. is the total number of the spin configurations, i.e., . The wave function (4) becomes exact in the isotropic spin limit ϕ = π+. Accordingly, the spin-spin correlations have long-range FM ordering for all three spin components: and for any . Taking the spin isotropic Hamiltonian at ϕ = π+ as an unperturbed one, the unperturbed ground state is given by Eq. (4). When , the perturbed Hamiltonian can be written as and the lowest-order energy correction is . Therefore, with increasing ϕ from π, the antiferromagnetic (AFM) fluctuations increase and the long-range FM ordering is weakened. Nonetheless, the correlations and retain the same asymptotic behaviors indicating the long-range FM ordering until , characterized by the saturation to a finite negative value in the limit ; whereas, decays in a power law with . This is because the AFM fluctuations are mainly introduced along the z-direction. It is confirmed by a slow decrease of the total spin as a function of ϕ, as seen in Fig. 2(c). We thus call this state FM-xy state. A collapse of the long-range FM ordering is detected by a drop-off of the total spin at , suggesting a first-order transition.

Figure 3

Schematic pictures of five states realized in the 1D KH model, except the TLL state at . In the phases (a,c,d), the spins lie mostly on the xy-plane and these states are rotational invariant around the z-axis. In phases (b,e), the spins almost align along the z-direction and a state having opposite spin directions is degenerate.

Ferromagnetic-z phase

A long-range FM ordered state stabilizes also at , where all the exchange tensors are FM since J + K is negative despite positive K. However, in contrast to the FM-xy state, spin alignment along the z-direction is favored because of the easy-axis-XXZ-like interactions . Therefore, the most dominant spin configurations are given by the highest sectors, i.e., the ground state is expressed aswith and . This wave function is exact in the isotropic limit of ϕ = π−. We call this type of long-range FM ordering FM-z state. Let us take Eq. (5) as the unperturbed ground state. For , the perturbed Hamiltonian can be written as . This is an Ising-like AFM correlation and it clearly disturbs the FM-z state. However, the lowest-order correction to the ground-state energy is . It means that the perturbation acts only gradually with being away from the isotropic point ϕ = π. As a result, the total spin is not that much reduced around ϕ = π in the FM-z phase. To estimate the lower bound of the FM-z phase, we shall define the FM-z order parameter. A state with long-range FM ordering is a state with broken spin symmetry along the z-direction; macroscopically, there are two degenerate ground states and . Applying open boundary conditions, one of the two ground states is picked imposing initial conditions on the calculation. The long-range ordered state is thus directly observable as a symmetry-broken state in our DMRG calculations. The order parameter is defined as

The finite-size scaling analysis is performed. The extrapolated values in the thermodynamic limit are shown in Fig. 4. Notably, the FM-z state survives even at , where a part of the interactions is AFM, i.e., J + K > 0. Nevertheless, the long-ranged FM order is drastically suppressed by the AFM fluctuations at and completely destroyed at . The order parameter has no jump at the transition point, suggesting a second-order phase transition.

Figure 4

(a) Finite-size scaling analyses of the FM-z order parameter with polynomial fitting functions. (b) Extrapolated values of the FM-z order parameter to the thermodynamic limit L → ∞ as a function of ϕ.

Spiral-xy phase

As shown above, the FM-z phase remains down to . Then, we ask which kind of phase comes next at . In this range, K is AFM, J is FM, and . On the x-links, the x-components of spins tend to be antiparallel due to the strong AFM interaction along the x direction; whereas, their y-components tend to be parallel due to FM J term. One can think about the y-links in the same way. Eventually, the spins lie on the xy plane and rotate by 90° from one site to the next. The magnetic unit cell is twice as large as the structural unit cell, i.e., including four lattice sites [see Fig. 3(a)]. We call this state spiral-xy state. To confirm this magnetic structure, we calculate the static spin structure factor, defined aswhere the length of the structural unit cell, i.e., two lattice spacings, is taken to be unity. The DMRG results with L = 60 cluster are shown in Fig. 5. Near the AFM Kitaev point , has a large peak indicating a periodicity of four lattice sites on the xy-plain. On the other hand, is almost zero for all q since the FM z interaction J is tiny compared to the dominant AFM x or y interactions . With increasing ϕ, the peak becomes lower and a peak in develops. Basically, the spins are tilted in one direction along the z-axis with keeping the periodicity on the xy-plain. Those peak heights are reversed around the transition point ϕ = 0.65π from the spiral-xy to the FM-z phases. We note that the height of the peak in coincides with the FM-z order parameter. As shown in Fig. 5(e), the spin-spin correlations decay as a power law for all the spin components in the spiral-xy phase. It means that the peak in disappears in the thermodynamic limit: the spiral-xy structure is not long-range ordered.

Figure 5

Static spin structure factors for (a) ϕ = 0.51π, (b) ϕ = 0.6π, (c) ϕ = 0.65π, and (d) ϕ = 0.75π. (e) Log-log plot of the spin-spin correlation functions as a function of distance at ϕ = 0.6π.

Néel-z ordered phase

Next, we turn to the parameter region (, ). We start from the SU(2) symmetric limit ϕ = 2π, where the system is the original 1D AFM Heisenberg model. The ground-state wave function can be exactly obtained and it is know to have . So, the first perturbative correction is given by , which provides an easy-axis anisotropy to the system. Therefore, the possibility of a continuous Néel-z order transition is conceived by analogy to the easy-axis anisotropic XXZ chain. For small , the magnetization is expect to grow gradually with

Let us confirm it numerically. In the long-range Néel-z ordered state, the translational symmetry is broken in a finite system due to the Friedel oscillation under the open boundary conditions, so that the Néel-z state can be directly observed by extracting one of the degenerate states, like in the FM-z state. Generally, the Friedel oscillations in the center of the system decay as a function of the system length. If the amplitude at the center of the system persists for arbitrary system lengths, it corresponds to a long-range ordering. Here, we are interested in the formation of alternating spin flip along the z-direction. Thus, the Néel-z order parameter is defined as

This quantity is equivalent to the magnetization M in the thermodynamic limit. The finite-size scaling analyses of was performed using the results for systems with up to and the extrapolated values to the thermodynamic limit were obtained. They are shown in Fig. 6. We indeed see a slow increase of near the SU(2) symmetric limit . Further with decreasing ϕ, the order parameter develops up to and drops down to 0 at .

Figure 6

(a) Finite-size scaling analyses for the Néel-z order parameter with polynomial fitting functions. (b) The extrapolated values of to the themodynamic limit as a function of ϕ. (c) Central charge, calculated with periodic chains, as a function of ϕ for several chain lengths.

In fact, this lower ϕ-boundary of the Néel-z ordered phase can be estimated analytically. At ϕ = tan-1 , the exchange term disappears in the Hamiltonian (1) and the system is just written as a sum of double-spin-flip and Ising parts

Each of the partition Hamiltonians and is exactly solvable. For , the system is regarded as noninteracting fermions with “pair hopping” and the ground-state wave function iswhere is summed over all possible spin configurations created by the double spin flips starting from or (), i and γ( =+ or −) are taken to create the configuration . While, for , the system is a simple Ising one and the ground-state wave function iswhere is the vacuum state. The ground-state wave functions (11) and (12) are orthogonal since they do not share the same spin configurations. It means the ground state is two-fold degenerate at the critical point . The ground-state energy is . This degeneracy is lifted away from . Equation (12) is the ground state for larger-ϕ side, namely, in the Néel-z phase; while, Eq. (11) for lower-ϕ side (see Sec. 6). This is consistent with the fact that reaches exactly 0.5 at , as shown in Fig. 6(b). Therefore, the ground-state wave function is completely changed at and it clearly suggests a first-order transition. This is also confirmed by the jump of as well as by the singularity in at .

Staggered-xy ordered phase

Let us then consider a region at , where the signs of Heisenberg and Kitaev terms counterchange from the spiral-xy state. This may mean the local spin structure is similar to the spiral-xy state. As discussed above, the system has a first-order transition at . In the lower-ϕ vicinity () the ground state is exactly given by Eq. (11), which derives asymptotic behaviors of the spin-spin correlations: , , , and . Since the system is rotation invariant around the z-axis, the coefficients α and β can take arbitrary real numbers under a condition . As illustrated in Fig. 3(d), all the spins lie on the xy-plane and the magnetic unit cell contains four lattice sites as in the spiral-xy state. However, the crucial difference from the spiral-xy state is that this xy state exhibits a long-range ordering, which is clearly indicated by the correlation functions. We call this state staggered-xy state, since it resembles a Néel-like state with a period of four sites.

We then investigate the ϕ-dependence of this xy state. Applying a staggered field along the x-direction on both open system edge sites, the presence or absence of the long-range staggered-xy order can be determined by studying the decay of the x-component of local spin as the Friedel oscillation. Thus, the staggered-xy order parameter is defined at the center of the system as

In the upper and lower vicinity of , and are obtained from Eqs (11) and (12), respectively. For the other ϕ values the profile of the x-component of local spin and the finite-size scaling of the order parameter are shown in Fig. 7(a,b). The extrapolated values of to the thermodynamic limit are plotted in Fig. 7(c). We can see that the order parameter jumps at the both phase boundaries suggesting first-order transitions. The collapse of the staggered-xy ordering at the lower boundary is clearly confirmed by the profile of , as shown in Fig. 7(a).

Figure 7

(a) The x-component of local spin, where opposite magnetic fields ±10 are applied at either end of the system. (b) Finite-size scaling analyses for the staggered-xy order parameter with polynomial fitting functions. (c) The extrapolated values of to the themodynamic limit as a function of ϕ.

Tomonaga-Luttinger-liquid phase

The remaining region is . At ϕ = 0 the system is equivalent to the spin-isotropic AFM Heisenberg chain, which is a gapless spin liquid. It is known that the low-energy physics is described by a Tomonaga-Luttinger (TL) model with a boson field, which is equivalent to the unity central charge (c = 1) conformal field theory (CFT)[49]. Let us now consider what happens when we move away from ϕ = 0. At the deviating interactions from the spin-isotropic Hamiltonian are written as

Obviously, these interactions cannot produce any explicit magnetic order, since they only enhance the quantum fluctuations. However, it is a nontrivial question whether the c = 1 CFT is conserved. To examine it, we directly calculate the central charge, which can be obtained from the von Neumann entanglement entropy in the DMRG procedure as follows: Let us consider a quantum 1D periodic system with length L. The von Neumann entanglement entropy of its subsystem with length l is given as , where is the reduced density matrix of the subsystem and ρ is the full density matrix of the whole system. Using the CFT, the entropy of the subsystem with length l for a fixed system length L has been derived[50-52]:where s1 is a non-universal constant. In the DMRG calculations this quantity can be accurately estimated by the second derivative of Eq. (15) with respect to l, namely,

The results for periodic chains with L = 10, 20, and 30 are plotted in Fig. 8. We see a quick convergence to c = 1 with increasing system size in the whole range of . Hence, we confirm the system remains in the TL liquid state as far as both J and K are AFM. Similarly, we studied the central charge in the Néel-z phase. The results for periodic chains with L = 16, 24, and 32 are plotted in Fig. 6(c). Although c = 0 should be retained in the Néel-z phase, we see a very slow converge to c = 0 with increasing L near ϕ = 2π, reflecting the tiny magnetization.

Figure 8

Central charge, calculated with periodic L = 10, 20, and 30 chains, as a function of ϕ.

Kitaev points

As discussed above, the Kitaev points are singular, not adiabatically connected to the neighboring phases. At these points , the Heisenberg interaction J vanishes in our Hamiltonian Eq. (1) and it is reduced to

In order to make the notation easier, we here focus on the case of ϕ = −π/2. Note that the case of can be similarly considered. It is convenient to fermionize the Hamiltonian (17). By applying the Jordan-Wigner transformation, we obtainwhere c and c are the fermion annihilation operators at the black and white labeling sites in the i-th unit cell, illustrated in Fig. 1(d), respectively. Generally, it is possible to write fermion operators in terms of Majorana fermions defined as:

With these operators Eq. (19) is transformed to

We immediately notice that the Majoranas b1, w2 are not included in the Hamiltonian (22). We then define new, non local, fermion operators aswhere the indices i run over the unit cells instead of over the lattice sites [see Fig. (1d)]. Using these operators, (22) can be expressed as

This describes a p-wave paired superconductor[53]. The exact solution for the ground state is easily obtained by a Fourier transformation:where and . Using a Bogoliubov transformation, this can be diagonalized and the quasiparticle excitation is given by

Thus, we can confirm the system is still in the gapless region. The ground-state energy is calculated as

To consider the meaning of the Majoranas b1, w2 which do not appear explicitly in the Hamiltonian (22), we take a possible recombination of them into new fermion operators:

We can now trace these back to the initial spin operators. Using Eqs (20) and (21), we get

Moreover, the inverse of Jordan-Wigner transformation used to obtain Eq. (19) derives

Thus, the Majorana operators can be expressed by spin operators, like

We now confirm that the system has a free spin per unit cell at the Kitaev points. Therefore, the dimensionality of the ground-state manifold, or number of zero Majorana modes, is for periodic chain and for open chain. Some more details are given in the Appendix.

Phase diagram

In Fig. 9(a) we summarize the ground-state phase diagram as a function of ϕ in a pie chart form. Notably, each phase has another exactly-symmetrically-placed phase in the pie chart, i.e, TLL and FM-xy, Néel-z and FM-z, spiral-xy and staggered-xy. The paired phases have the common features. The TLL and FM-xy states, where the exchange components are either all AFM or all FM, are basically understood within the framework of the isotropic Heisenberg chain; the FM-xy state is long-range ordered and the TLL state is critical. The Néel-z and FM-z states are described by the easy-axis XXZ Heisenberg chain. The dominant features are determined by the Ising term. Depending on the sign of the Ising term, the spins are aligned ferromagnetically or antiferromagnetically along the z axis. The spiral-xy and staggered-xy states can be basically interpreted as the easy-plane XXZ Heisenberg chain affected by the double-spin-flip term. The system has a four-site periodicity, and exhibits a critical behavior and a long-range ordering for the former and latter states, respectively, in analogy to the relation between the TLL and FM-xy states.

Figure 9

(a) Ground-state phase diagram of the 1D KH model. (b) Ground-state phase diagram of the honeycomb-lattice KH model, where the same notations as in Eq. (1) are used for K and J. (c) Brick-wall lattice created by coupling the 1D KH chains by the z-bond. This is topologically equivalent to the honeycomb lattice.

It is now possible to establish a connection between our ground-state phase diagram and that of the two-dimensional KH model on a honeycomb lattice [Fig. 9(b), normalized to the current notation from[10]]. To consider the phase-to-phase correspondence, it is helpful to regard the honeycomb lattice as coupled KH chains, as shown in Fig. 9(c). Here, the interchain coupling is equivalent to the z-bond in the honeycomb-lattice KH model, written aswhere k, l are summed over all connected sites between the KH chains.

Hereafter, we report how the introduction of the z-bonds coupling affects the (quasi-)ordered states present in the one-dimensional case, mapping them to those of the honeycomb system:

At the interchain coupling is AFM. Our TLL state has no long-range order but strong AFM fluctuations. Therefore, a Néel order is intuitevely expected once the chains are antiferromagnetically coupled. Hence, our TLL phase corresponds to the Néel phase in the honeycomb case.

At our system is in the FM-z state. While considering this interval, we need to examine two different cases as the sign of the interaction on the z-bonds changes from AFM to FM. We obtain a zigzag state in the honeycomb case when the FM-z chain state is coupled by AFM interchain couplings; whereas, a FM state in the case of FM interchain couplings. The change in sign of the interchain coupling is assumed to happen at , where J + K = 0. This ϕ value is reasonably close to the transition point between the zigzag and FM states in the honeycomb case. Thus, our FM-z phase is distributed to either zigzag or FM phases of the honeycomb KH model depending on the AFM or FM nature on the z-bond.

At our system is in the Néel-z state. In the same way as above, we need to consider the cases of AFM and FM z-bonds. We easily find that our Néel-z state is base for the Néel and stripy states of the honeycomb KH model[46]. Namely, our Néel-z phase is distributed to either Néel phase in the presence of AFM interchains coupling or the stripy phase in the presence of FM interchains coupling of the honeycomb KH model. The value of giving is again near the transition point between the Néel and stripy states in the honeycomb case.

Therefore, it is possible to understand all the long-range ordered phases of the honeycomb-lattice KH model in the framework of the coupled KH chains. In brief, an extracted zigzag lattice line from any ordered state of the honeycomb-lattice KH model corresponds to one of the phases in the 1D KH chain; and the remaining degrees of freedom derive from the fact that the zigzag lines are coupled either ferromagnetically or antiferromagnetically by .

Two xy phases in our phase diagram are left. Since these states are stabilized by the dominant xy-Kitaev term, they are connected to the Kitaev spin liquid states in the honeycomb case. In hindsight, this means that our staggered-xy ordered state can collapse to the spin liquid state due to strong FM fluctuations along the z-direction.

Low-lying excitations

To examine the low-energy excitations for each phase, we calculated the dynamical spin structure factors , defined in Eq. (2). Under periodic boundary conditions the momentum is taken as since the unit cell contains two lattice sites and the number of unit cells is in a system with L sites. We study chains with L = 24 using the Lanczos ED method and the results are shown only for as .

As confirmed above, the low-energy physics at is described by a gapless TL model. The low-lying excitations are expected to be basically equivalent to those of the easy-plane XXZ chain because [see Eq. (2)]. If we ignore the double-spin-flip term , the main dispersion (lower bound of the spectrum) is given bywith . This was obtained by using the transfer matrix of the six-vertex model[54,55]. Nevertheless, the effect of is unknown. Therefore, to investigate the effect of we employ a standard SWT for the bipartite system. Using the Holstein-Primakovff representation the spin operators of for the black labeling sites in Fig. 1(d) are replaced as and ; similarly, for the white labeling sites and , where a and b are canonical boson annihilation operators at site i. Then, applying the Fourier transform we obtain:

This off-diagonal term gives a splitting in the main dispersion (36). Then, we can speculate the total dispersions as

Figure 10(a) shows the dynamical structure factors S(q, ω) and S−(q, ω) for . The largest peak appears in reflecting the AFM fluctuations. The intensities in S(q, ω) are weaker than those in S−(q, ω) due to the easy-plane xy anisotropy. The split dispersions are well described by Eq. (38). The splitting is largest in the vicinity of the Kitaev point (J → 0+, K → 1):

Figure 10

Dynamical structure factors calculated with a periodic 24-site chain for (a) TLL , (b) spiral-xy [], and (c) FM-z [] states. The left and right panels show S(q, ω) and S−(q, ω), respectively. The red lines are analytical dispersions: Eq. (38) for and Eq. (41) for . The arrow indicates an appearance of the spiral-xy fluctuation in the FM-z state.

The main dispersion (lower bound of the spectrum) of the isotropic Heisenberg model at ϕ = 0 (K = 0) is also reproduced by Eq. (38).

Across the Kitaev point from the TL to spiral-xy phases the momentum of the Fermi point is shifted from q = π to q = 0. Thus, in the vicinity of the Kitaev point (J →0−, K → 1), the momentum of the excitation dispersion is shifted from Eq. (39) by :for S−(q, ω). In contrast, S(q, ω) exhibits a q-independent continuum between ω = 0 and ω = K since the z-component of exchange interaction is much smaller than the other interactions. With increasing ϕ from , the FM Ising interaction and double-spin-flip fluctuations increase and the exchange interaction decreases becoming zero at the critical boundary to the FM-z phase. Figure 10(b) shows the dynamical structure factors S(q, ω) and S−(q, ω) for . The main dispersion in S−(q, ω) is scaled as . Although the largest peak at indicates a four-site periodicity in the spiral-xy state, it diminishes with increasing system size because of no long-range ordering. While in S(q, ω), a dispersion appears and a peak toward the FM-z state develops at . Both the dispersion widths are roughly scaled by .

Figure 10(c) shows the dynamical structure factors and for . Near ϕ = π () in the FM-z phase, the system is effectively described by a FM XXZ Heisenberg chain with easy-axis anisotropy and the ground state is approximately expressed by Eq. (5). An applied operator does not change the ground state in Eq. (2) and the final state has only a zero energy excitation from . Therefore, S(q, ω) has a very sharp peak at and almost no spectral weight appears at the other momenta. This peak keeps its weight constant with increasing system length, indicating the long-range FM-z ordering. On the other hand, in Eq. (2) the operator dopes one magnon into the FM-z alignment so that the SWT is expected to give a good approximation for the excitation dispersion of S−(q, ω):

We can confirm that a gap opens at q = 0 reflecting the easy-axis anisotropy. With approaching the neighboring spiral-xy phase, the double-spin-flip fluctuations grow gradually in influence; accordingly, the q = 0 peak in S(q, ω) shrinks and a peak develops at q = π, in S−(q, ω).

Figure 11(a) shows the dynamical structure factors S(q, ω) and S−(q, ω) for . Near ϕ = π () in the FM-xy phase, the system is effectively described by a FM XXZ Heisenberg chain with easy-plane anisotropy and the ground state is approximately expressed by Eq. (4). Thus, the excitation spectrum is expected to be gapless. Since the total S of the ground state is zero, unlike the case of FM-z state both S(q, ω) and S−(q, ω) have the same excitation dispersion as

Figure 11

Dynamical structure factors calculated with a periodic 24-site chain for (a) FM-xy [], (b) staggered-xy [], and (c) Néel-z [] states. The left and right panels show S(q, ω) and S−(q, ω), respectively. The red lines are analytical dispersions: Eq. (42) for , Eq. (43) in S(q, ω) and Eq. (44) in S−(q, ω) for , and Eq. (45) for .

Note that a sharp peak at in S−(q, ω) keeps its weight constant with increasing system length, indicating the FM-xy long-range ordering. On the other hand, in the vicinity of the Kitaev point , the excitation dispersion of S−(q, ω) is well described by Eq. (39).

In the vicinity of the Kitaev point , the excitation spectra are the same as in the other Kitaev point , namely, Eq. (40) for S−(q, ω) and a q-independent continuum for S(q, ω). Whereas in the vicinity of the Néel-z ordered phase , the ground state is expressed by Eq. (11). Thus, the dispersions are described by a single magnon excitation:for S(q, ω), andfor S−(q, ω). Figure 11(b) shows the dynamical structure factors S(q, ω) and S−(q, ω) for . The excitation dispersions are well reproduced by Eqs (43) and (44). We also find a sharp peak at (q, ω) = (π, 0) in S−(q, ω), which indicates the staggered-xy long-range ordering. This peak keeps its weight constant with increasing the system length, in contrast to the similar peak for the spiral-xy state.

In the vicinity of the staggered-xy ordered phase , the ground state is expressed by Eq. (12). Accordingly, S(q, ω) has only a delta peak at (q, ω) = (0, 0). Whereas, S−(q, ω) is exactly explained by a single magnon dispersion. It is obtained by the SWT as

Although this is equivalent to Eq. (43), the spectral weight is uniform for all q values. Since the transition at is first ordered, there is no peak indicating a connection to the staggered-xy ordered phase. On the other hand, near ϕ = 2π, the system can be basically regarded as an easy-axis AFM XXZ Heisenberg chain so that the excitation dispersion is described by Eq. (38) for both S(q, ω) and S−(q, ω). Figure 11(c) shows the dynamical structure factors S(q, ω) and S−(q, ω) for the intermediate region . The main dispersion of S−(q, ω) is basically described by Eq. (45) but some features from Eq. (38) seem to be somewhat mixed.

Discussion

We have established the presence of a variety of phases in the 1D KH system. Especially, it is surprising that most of the ϕ ranges are covered by long-range ordered phases despite considering a pure 1D system. In this context, we now consider the K-intercalated α-RuCl3, namely, K0.5RuCl3. One should be aware of the fact that several different parameter sets have been suggested for undoped α-RuCl3: (i) K = −5.6, J = 1.2 ()[25], (ii) K = 7.0, J = 4.6 ()[22], (iii) K = 8.1, J = 2.9 ()[22], (iv) K = −6.8, J = 0 ()[56] in unit of meV. If we assume that the charge ordering pattern in K0.5RuCl3 is that illustrated in Fig. 1(c), the parameter sets (i)–(iv) correspond to the staggered-xy, FM-z, spiral-xy, and FM Kitaev point, respectively. In practice, the charge ordering could cause a significant change of the parameter since they are very sensitive to the Ru-Cl-Ru bond angle[25]. In other words, once the magnetic properties of K0.5RuCl3 are observed, we may easily speculate the possible parameter set of K0.5RuCl3 and even the charge ordering pattern by comparing then to our rich phase diagram. To gain deeper insights, theoretical and experimental studies under magnetic fields are also required.

It is relevant to seek other possible realizations of 1D KH system. Even if the Kitaev interaction in 1D systems is present, it is considerede to be very small. However, as shown above, even a tiny Kitaev interaction can stabilize the ordered state. It might also be intriguing to reconsider quasi-1D materials having two sublattices like Ni2(EDTA)(H2O)4, [Ni(f-rac-L)(CN)2], LiCuSbO4, and Rb2Cu2Mo3O12 from the point of view of the 1D Kitaev system.

The ϕ-dependent phase diagram of our model is similar to that for the honeycomb-lattice KH model. Remarkably, all the magnetically ordered states of the honeycomb-lattice KH model can be interpreted in terms of the coupled 1D KH chains. In other words, the key elements to derive the magnetic ordering in the 2D honeycomb-lattice KH model are already contained in the 1D KH chain. This will give us a deeper insight in the understanding of transition from an ordered state to a disordered state such as the Kitaev QSL.

Methods

In this section, we present a detailed derivation of (33), (34), and the related additional information. We used a Jordan-Wigner transformationto rewrite the spin operators in Eq. (17) by fermion operators, and we obtained Eq. (19).

By applying the inverse of the Jordan-Wigner transformationthe fermion operator was traced back to the initial spin operators, and we obtained Eq. (32) from Eq. (31).

We can also derive the same results as Eqs (33) and (34) by taking a different combination of the Majoranas from Eqs (29) and (30):

From Eqs (20) and (21), we have

Using (50)–(52), we obtain

Then, we obtain

In a similar way, we find

Actually, Eqs (59) and (60) are identical to Eqs (33) and (34).