Chinese Physics Letters, 2023, Vol. 40, No. 12, Article code 126403 A Possible Quantum Spin Liquid Phase in the Kitaev–Hubbard Model Shaojun Dong1, Hao Zhang2,3, Chao Wang1, Meng Zhang2,3, Yong-Jian Han1,2,3*, and Lixin He1,2,3,4* Affiliations 1Institute of Artificial Intelligence, Hefei Comprehensive National Science Center, Hefei 230088, China 2CAS Key Laboratory of Quantum Information, University of Science and Technology of China, Hefei 230026, China 3Synergetic Innovation Center of Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China 4Hefei National Laboratory, University of Science and Technology of China, Hefei 230088, China Received 12 September 2023; accepted manuscript online 9 November 2023; published online 6 December 2023 *Corresponding authors. Email: smhan@ustc.edu.cn; helx@ustc.edu.cn Citation Text: Dong S, Zhang H, Wang C et al. 2023 Chin. Phys. Lett. 40 126403    Abstract The quantum spin liquid (QSL) state of Kitaev-like materials, such as iridium oxides $A_2$IrO$_3$ and $\alpha$-RuCl$_3$, has been explored in depth. The half-filled Kitaev–Hubbard model with bond-dependent hopping terms is used to describe the Kitaev-like materials, which is calculated using the state-of-the-art fermionic projected entangled pair states method. We find a QSL phase near the Mott insulator transition, which has a strong first-order transition to the semi-metal phase with the decrease of Hubbard $U$. We suggest that a promising approach to finding QSL states is by finding iridium oxides that are near the Mott insulator transition.
cpl-40-12-126403-fig1.png
cpl-40-12-126403-fig2.png
cpl-40-12-126403-fig3.png
cpl-40-12-126403-fig4.png
cpl-40-12-126403-fig5.png
cpl-40-12-126403-fig6.png
DOI:10.1088/0256-307X/40/12/126403 © 2023 Chinese Physics Society Article Text A quantum spin liquid (QSL)[1-3] state is a quantum state that lacks any long-range magnetic order even down to zero temperature. QSLs have nontrivial topological properties that may host exotic excitations with fractional statistics, such as spinons and visions, which may have important applications in quantum computing[4,5] and may play a crucial role in high-temperature superconductivity. The Kitaev model[5] is an exactly solvable model on a 2D honeycomb lattice, which hosts a QSL ground state. Several iridium oxides $A_2$IrO$_3$, as well as $\alpha$-RuCl$_3$, have been proposed to exhibit the Kitaev QSL state.[6-14] These materials have a honeycomb structure and the strong spin–orbit coupling leads to an effective $j_{\rm eff}=\frac{1}{2}$ spin model with bond-dependent anisotropic exchange interactions,[8,15] which are the essential ingredients of the Kitaev model. In addition to the Kitaev exchange interactions, there are also Heisenberg interactions in these materials.[16] The Kitaev–Heisenberg model has been intensively studied and it has been shown that the QSL state can only survive within a rather small parameter space.[16-20] Indeed, Na$_2$IrO$_3$ and $\alpha$-RuCl$_3$ were found to have a zigzag antiferromagnetic (AFM) order through resonant x-ray magnetic scattering and inelastic neutron scattering experiments. Tremendous efforts have been made to find the QSL state in these materials.[6,9,21-25] Recently, evidence of the QSL phase in the presence of large external magnetic fields was described.[26,27] An important question is that, given the extremely small parameter space for the QSL state in the Kitaev–Heisenberg model, is it possible to find Kitaev QSL state in real materials when no external field is applied? The Kitaev–Hubbard model provides a more realistic framework for describing iridium oxides. When the Hubbard $U$ is small, higher-order interactions become important, which may introduce exotic states. The Kitaev–Hubbard model has been studied by mean-field theories.[28-30] It has been shown that there is a QSL phase in the region of $t' < t$ when $U$ is small, where $t$ and $t'$ are the isotropic and spin-dependent hopping terms, respectively. A further decrease in $U$ results in a semi-metal (SM) phase. However, these calculations were based on mean-field approximations,[28-30] which need to be examined with more rigorous methods. Furthermore, these studies focus on the $t' < t$ region and the phase diagram for $t'>t$ is missing. Experimentally, iridium oxide materials, e.g., Na$_2$IrO$_3$ and $\alpha$-RuCl$_3$, are believed to have strong spin-dependent hopping terms[9,31-33] in the $t'>t$ region. The projected entangled pair states method (PEPS),[34-39] and its generalization to fermionic systems (fPEPS)[40-44] provide systematically improvable variational wave functions for the many-body problems, which allow more rigorous treatment of the Kitaev–Hubbard model. In this Letter, we apply this recently developed and highly accurate fPEPS method to explore the phase diagram of the half-filled Kitaev–Hubbard model. The results show that the QSL state is absent in the $t' < t$ region, in contrast to the previous mean-field results.[28-30] Instead, we find a QSL phase in the $t'>t$ region when $U$ is small. We show that the phase transition from the SM phase to the QSL phase is a first-order transition, whereas the QSL-zigzag transition is a continuous transition. Given that iridium oxide materials and $\alpha$-RuCl$_3$ are in the $t'>t$ region, it is possible to find suitable materials that may host the QSL state. Although the multi-orbital Hubbard model offers a more realistic description of iridium oxides,[45] its calculation using the PEPS method remains exceedingly challenging at present. Consequently, we turn to the single-orbital Kitaev–Hubbard model, i.e., a simplified version of the multi-orbital Hubbard model, to provide valuable insights into the complex physics of these materials. The Kitaev–Hubbard model reads \begin{align} H=\sum_{\langle i ,j \rangle_\alpha,s} \Big\{\hat{c}_{i,s}^† \Big(\frac{t+t'\sigma^\alpha}{2}\Big)\hat{c}_{j,s} + {\rm H.c.}\Big\} + U\sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow}, \tag {1} \end{align} where $\hat{c}_{i,s}$ represents the annihilation operator that destroys an electron with spin $s$ at site $i$, and $\hat{n}_{i,s}=\hat{c}_{i,s}^† \hat{c}_{i,s}$ is the number operator; $\sigma^\alpha$, with $\alpha=x$, $y$, $z$ are the Pauli matrices; $\langle i,j \rangle_\alpha$ denotes the nearest-neighbor pairs in the three hopping directions of the lattice, as sketched in Fig. 1. The $t$ and $U$ terms are the hopping and the on-site Coulomb interaction terms in the normal Hubbard model, whereas the $t'$ is the spin-dependent hopping due to spin-orbit coupling.
cpl-40-12-126403-fig1.png
Fig. 1. (a) The honeycomb lattice with bond-dependent interactions. (b) The honeycomb lattice is mapped to a PEPS defined on a square lattice, where two physical indices on the honeycomb lattice are combined in a single site.
The fPEPS method[34-40] is one of the most promising approaches for investigating strongly correlated electron systems. We simulate the honeycomb lattice with an open boundary condition using a square tensor network. We map the honeycomb lattice from Fig. 1(a) to a square lattice of Fig. 1(b), where a unit cell of the honeycomb lattice containing two sites is treated as a single site in the square lattice (except those sites at the corner, which remain a single site in the square lattice). Each tensor network site contains two physical indices. The physical properties can still be easily calculated on the physical (honeycomb) lattice. We optimize the wave functions using the so-called fPEPS++ method developed by our group,[39,40,46] i.e., the fPEPS wave functions are optimized via a stochastic gradient method, whereas the energy and energy gradients are calculated using a Monte Carlo sampling method. This method significantly reduces the computational complexity with respect to the bond dimension $D$, thereby allowing a much larger bond dimension to be used and resulting in highly accurate and converged results for large finite systems. All the parameters in the fPEPS wave functions are independent and subject to optimization. We first obtain the quantum wave functions with the simple update method. The fPEPS wave functions are further optimized using the stochastic gradient method until the results fully converge.[39,40,46] In our calculations, the bond dimension $D=14$ and the truncation bond dimension $D_c=42$ are used, which show good convergence. The magnetic ordering of the system is determined by the spin structure factor, which is defined as \begin{align} S(\{{\boldsymbol k}_\alpha\})=\sum_{i,j,\alpha} e^{-i{{\boldsymbol k}_\alpha} \cdot ({{\boldsymbol r}_i}-{{\boldsymbol r}_j})}S^\alpha_i \cdot S^\alpha_j , \tag {2} \end{align} where ${{\boldsymbol r}_i}$ is the coordinate of site $i$ on the lattice, and $S^\alpha_i$ is the spin operator at site ${{\boldsymbol r}_i}$, along the $\alpha$ direction ($\alpha=x,y,z$). The AFM order is characterized by the non-zero value of the spin structure factor at ${\boldsymbol k}_x={\boldsymbol k}_y={\boldsymbol k}_z=(\frac{4\pi}{3},\,0)$, whereas the zigzag order is detected by a non-zero value of the spin structure at ${{\boldsymbol k}_x}=(\frac{2\pi}{3},\,0)$, ${{\boldsymbol k}_y}=(-\frac{\pi}{3},\,\frac{\sqrt{3}\pi}{3})$ and ${{\boldsymbol k}_z}=(\frac{\pi}{3},\,\frac{\sqrt{3}\pi}{3})$. The SM phase is characterized by a vanishing charge gap \begin{align} \varDelta=E_{N+1}+E_{N-1}-2E_N, \tag {3} \end{align} where $E_N$ is the total energy of the system with $N$ electrons. When the system is in the insulating phase and none of the aforementioned magnetic orders are present, the system could represent a highly promising candidate for a QSL phase. We first discuss the Kitaev–Hubbard model Eq. (1) in the large-$U$ limit. Without loss of generality, we take $t=1$ throughout this study. At half filling, the model can be reduced to the Kitaev–Heisenberg spin model to the leading order of $1/U$,[28,30] \begin{align} H_{\rm eff}=\sum_{\langle i ,j \rangle_\alpha} \Big(\frac{(1-t'^2)}{U} {\boldsymbol S}_i \cdot {\boldsymbol S}_j + \frac{2t'^2}{U}S_i^\alpha S_j^\alpha\Big), \tag {4} \end{align} where $S_i^\alpha$, $\alpha=x,y,z$, are the spin operators at site $i$, ${\boldsymbol S}_i=(S_i^x,\,S_i^y,\,S_i^z)$, and $\langle i,j \rangle_\alpha$ denotes the nearest-neighbor pairs in the three hopping directions of the lattice (see Fig. 1). The Kitaev–Heisenberg model has been studied intensively and the phase diagram of the model is well known.[16,18,19,47-55] In Ref. [17], the authors extended the original model to its full parameter space, i.e., \begin{align} H_{\rm KH}=\sum_{\langle i ,j \rangle_\alpha} [J {\boldsymbol S}_i \cdot {\boldsymbol S}_j + K S_i^\alpha S_j^\alpha ], \tag {5} \end{align} where $J=\cos\phi$ is the Heisenberg coupling strength and $K=2\sin\phi$ is the Kitaev coupling strength. The phase angle $\phi$ may vary from 0 to $2\pi$. The phase boundaries have been obtained by exact diagonalization of the Hamiltonian on a 24-site hexagonal lattice with periodic boundary conditions. Compared with the large-$U $ effective Hamiltonian of Eq. (4), the angle $\phi$ can be related to $t'$ as $\cot \phi=\frac{1-t'^2}{t'^2}>-1 $ and $\sin\phi=\frac{t'^2}{U}>0$, and therefore we have a constrain of $0\le \phi < 3\pi/4$ for Eq. (5). In this parameter region, there is a QSL phase on $\phi\in (88^{\circ},\,92^{\circ})$, corresponding to $t' \sim (0.99970,\,1)$ in Eq. (4), which is almost a single point in the parameter space. It has been shown that when $t' < 1$, the ground state is an AFM phase and it is the zigzag phase for $t'>1$. Therefore, in the large $U$ limit, the AFM phase and the zigzag phase are separated by a QSL state that survives (almost) only at the $t'=1$ line. The decrease in $U$ may introduce higher-order spin interactions,[28] which may stabilize the QSL phase in a larger region. Several studies have shown some insight into this problem, where the authors claim that an algebraic QSL state lies between $t'\sim 0.7$ and $t'=1$ when $U $ is small.[28-30] However, these calculations were based on mean-field approximations,[28-30] which need to be examined by more rigorous method. Furthermore, the studies focused on the $t' < 1$ region and the phase diagram for $t'>1$ was not studied.
cpl-40-12-126403-fig2.png
Fig. 2. The ground state phase diagram of the Kitaev–Hubbard model on the $t'$–$U$ plane, where we have set $t=1$. Four phases have been identified, including an SM phase, an AFM phase, a zigzag phase and a QSL phase. The scatters represent the parameters calculated using the fPEPS method. The phase boundary between SM and QSL state is around $U=4.5$.
We calculate the phase diagram of the Kitaev–Hubbard model in the $t'$–$U$ plane using the fPEPS method and the results are shown in Fig. 2. Four phases have been identified in the phase diagram. On the left side of the diagram, where $U$ is small, there is a large SM phase that adiabatically connects to the phase at $t'=0$, and $U=0$, i.e., the electronic structure of graphene. In the large $U$ limit, the system is in an AFM phase for $t' < 1$ and a zigzag phase when $t'>1$. Remarkably, within the parameter range of $t'\in (1,\,1.2)$ and $U\in (4.5,\,7)$, both the AFM and the zigzag spin orders vanish. As $U>7$, the region of this phase reduces to the $t'=1$ line. Since this phase adiabatically connects to the QSL phase in the large $U$ limit of the Kitaev–Hubbard model[17] and the possible ordered phases from the spin models[16,18,19,47-55] and the electron models[28-30] have been excluded by the order parameters, it is plausible to consider this phase as a candidate QSL phase. The SM phase is accompanied by the vanishing of the charge gap. Figures 3(a) and 3(b) depict the charge gaps at $t'=0.95$ and 1.05 on the honeycomb lattices of 126, 198, and 286 sites. The charge gaps decrease with decreasing $U$. For both $t'=0.95$ and $t'=1.05$, the charge gaps become zero between $U=4.0$–4.5, which are the phase boundaries between the SM and Mott insulator phases. The SM phase is further ensured with the vanishing of the local magnetic moment.
cpl-40-12-126403-fig3.png
Fig. 3. The charge gaps $\varDelta$ on different lattice sizes for (a) $t'=0.95$ and (b) $t'=1.05$.
cpl-40-12-126403-fig4.png
Fig. 4. The finite size scaling of the AFM order parameters for (a) $U=5$ and (b) $U=8$, and the finite size scaling of the zigzag order parameters for (c) $U=5$ and (d) $U=8$.
We now focus on the insulating region. To determine the magnetic order, we calculate the AFM and zigzag order parameters in the thermodynamic limit via finite size scaling. Figures 4(a) and 4(c) depict AFM and zigzag order parameters for $U=5$, with $t'=0.95$, 1, 1.05, and 1.2 versus the square root of the number of lattice sites used in the calculations, whereas Figs. 4(b) and 4(d) show the results for $U=8$. To reduce the boundary effects, the order parameters are calculated using only the central region of the lattice.[39,46,56] For $U=5$ and $t'=0.95$, the system shows a finite AFM order in the thermodynamic limit. When $t'=1.2$, the system shows a zigzag order. However, for $t'=1$ and $t'=1.05$, both AFM and zigzag orders vanish in the thermodynamic limit. In contrast, for $U=8$ and $t'=0.95$, the system also has an AFM order. In addition, for $t'=1.05$ and $t'=1.2$, the system shows a zigzag order. Only when $t'=1$, do both AFM and zigzag orders vanish, as expected from the Kitaev–Heisenberg model.[17]
cpl-40-12-126403-fig5.png
Fig. 5. (a) AFM order in the thermodynamics limit as a function of $U$ for $t'=0.95$. (b) AFM order and zigzag order versus $t'$ at $U=5$.
Figure 5(a) depicts the AFM order parameter (in the thermodynamic limit) as a function of $U$ for $t'=0.95$. The AFM order disappears at approximately $U \sim4.2$, which is coincident with the disappearance of the charge gap $\varDelta$. This result suggests that there is no other phase between the SM and AFM phases. The calculated phase boundary between the SM and AFM phases is in agreement with the mean field results[28-30] for $t' < 0.8$. However, it shows a distinguishable difference for $t'$ between 0.8 and 1. Mean field calculations suggest that there is a QSL state in this region,[28-30] which is absent in more rigorous fPEPS calculations. Figure 5(b) depicts the AFM and the zigzag order parameters along the line of $U=5$. The system has an AFM order when $t' < 1$ and a zigzag order when $t'>1.15$. Both orders disappear in between, which suggests that it is possibly a QSL phase. To determine the order of the phase transitions, we calculate the first-order energy derivative with respect to $U$, $\frac{d E}{d U}$, using the Hellmann–Feynman theorem, i.e., \begin{align} \frac{d E}{d U}=\langle \varPsi |\frac{1}{N} \sum_i {\partial H \over \partial U} |\varPsi \rangle= \langle \varPsi | \frac{1}{N}\sum_i \hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}|\varPsi \rangle, \tag {6} \end{align} where $|\varPsi\rangle$ is the ground state fPEPS wave function, and $N$ is the total number of lattice sites that are used to calculate the total energy. To reduce the boundary effects, only the central region of lattice is used to calculate the total energies.[39,46,56] The results are shown in Fig. 6(a) for $t'=1.05$. The second derivative of the energy with respect to $U$ (by the finite difference method) is shown in Fig. 6(b). Clearly, there is a sharp discontinuity of $\frac{d E}{d U}$ at the SM–QSL transition around $U \sim4.5$, which suggests that this is a strong first-order transition. In contrast, the transition between the QSL and zigzag phases is rather continuous. There might be a small discontinuity in the second derivative of the energy, which is not captured due to the relatively large energy fluctuation near the phase transition, thereby hindering the accurate determination of the second derivative. Note that $\langle \hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}\rangle$ also characterizes the electron double occupancy on site $i$. The QSL–SM transition is driven by the sudden increase of the double occupancy. The QSL phase also benefits from the increase in electron double occupancy.
cpl-40-12-126403-fig6.png
Fig. 6. (a) The first-order energy derivative with respect to $U$, i.e., $\frac{d E}{d U}$, which is the same value of the double occupy number, and (b) the second-order energy derivative $\frac{d^2 E}{d U^2}$ on the lattices of different sizes with $t'=1.05$.
The values of the Heisenberg coupling $J$ and the Kitaev coupling $K$ in real materials, such as Na$_2$IrO$_3$ and $\alpha$-RuCl$_3$, have been estimated in Refs. [9,31-33]. It has been suggested that, in these materials, $|J| < |K|$ and $J/K < 0$, i.e., $t'>t$. We find a large region of zigzag phase in $t'>t$, $U>4.5$, which is consistent with the zigzag phase in Na$_2$IrO$_3$ and $\alpha$-RuCl$_3$, determined by the resonant x-ray magnetic scattering and inelastic neutron scattering experiments.[6,9,21,22] However, since iridium oxides are in the $t'>t$ region, it is a promising routine to find the QSL state in iridium oxides that are close to the Mott insulator transitions. Discussion and Outlook. In this study, we explore the ground state phase diagram of the Kitaev–Hubbard model under half-filling conditions. Our analysis reveals the presence of distinct phases, including the SM phase, AFM phase, and zigzag phase. Of particular significance, we identify the potential existence of a QSL phase, which emerges near the Mott insulator transition. This intriguing observation prompts the speculation that materials such as iridium oxides, which exhibit robust spin-dependent hopping and are close to Mott insulator transitions, may harbor QSL behavior. However, we note that direct evidence of the QSL phase is still lacking. The correlation functions undergo rapid decay and thus require extremely high accuracy at a large distance, which is presently beyond our capability. On the other hand, the multi-orbital Hubbard model is a more realistic model for the iridium oxides.[45] However, treating the multi-orbital Hubbard model using the PEPS method still presents considerable challenges. Very recently, there was a strong indication of the QSL phase in $\alpha$-RuCl$_3$ under external fields.[27] In this work, we only investigate the ground state under zero magnetic field. Exploring the model under a magnetic field presents a promising avenue for further research. Acknowledgements. This work was supported by the National Natural Science Foundation of China (Grant Nos. 12104433 and 11874343), and the Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0301200).
References Resonating valence bonds: A new kind of insulator?The Resonating Valence Bond State in La2CuO4 and SuperconductivitySpin liquids in frustrated magnetsFault-tolerant quantum computation by anyonsAnyons in an exactly solved model and beyondSpin Waves and Revised Crystal Structure of Honeycomb Iridate Na 2 IrO 3 Relevance of the Heisenberg-Kitaev Model for the Honeycomb Lattice Iridates A 2 IrO 3 α RuCl 3 : A spin-orbit assisted Mott insulator on a honeycomb latticeMagnetic order in α RuCl 3 : A honeycomb-lattice quantum magnet with strong spin-orbit couplingEvidence for a Field-Induced Quantum Spin Liquid in α - RuCl 3 Field-induced quantum criticality in the Kitaev system α RuCl 3 Phase diagram of α RuCl 3 in an in-plane magnetic fieldField-induced transitions in the Kitaev material α RuCl 3 probed by thermal expansion and magnetostrictionThermodynamic Perspective on Field-Induced Behavior of α RuCl 3 Electronically highly cubic conditions for Ru in α RuCl 3 Kitaev-Heisenberg Model on a Honeycomb Lattice: Possible Exotic Phases in Iridium Oxides A 2 IrO 3 Zigzag Magnetic Order in the Iridium Oxide Na 2 IrO 3 Quantum phase transition in Heisenberg-Kitaev modelFinite-temperature phase diagram of the Heisenberg-Kitaev modelOrder-by-disorder and spin-orbital liquids in a distorted Heisenberg-Kitaev modelLong-range magnetic ordering in Na 2 IrO 3 Direct evidence of a zigzag spin-chain structure in the honeycomb lattice: A neutron and x-ray diffraction investigation of single-crystal Na 2 IrO 3 Evidence for Magnetic Fractional Excitations in a Kitaev Quantum-Spin-Liquid Candidate α-RuCl3Rare-Earth Chalcohalides: A Family of van der Waals Layered Kitaev Spin Liquid CandidatesQuantum Spin Liquid Phase in the Shastry-Sutherland Model Detected by an Improved Level Spectroscopic MethodNeutron Spectroscopy Evidence for a Possible Magnetic-Field-Induced Gapless Quantum-Spin-Liquid Phase in a Kitaev Material $\alpha$-RuCl$_3$Gapless Spin Excitations in the Field-Induced Quantum Spin Liquid Phase of α RuCl 3 Stable Algebraic Spin Liquid in a Hubbard ModelDistinct-symmetry spin-liquid states and phase diagram of the Kitaev-Hubbard modelTopological phases of the Kitaev-Hubbard model at half fillingProximate Kitaev quantum spin liquid behaviour in a honeycomb magnetSpin waves and stability of zigzag order in the Hubbard model with spin-dependent hopping terms: Application to the honeycomb lattice compounds Na 2 IrO 3 and α - RuCl 3 Generic Spin Model for the Honeycomb Iridates beyond the Kitaev LimitUnifying projected entangled pair state contractionsA practical introduction to tensor networks: Matrix product states and projected entangled pair statesMatrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systemsAccurate Determination of Tensor Network State of Quantum Lattice Models in Two DimensionsRenormalization algorithms for Quantum-Many Body Systems in two and higher dimensionsGradient optimization of finite projected entangled pair statesGradient optimization of fermionic projected entangled pair states on directed latticesGrassmann tensor network states and its renormalization for strongly correlated fermionic and bosonic statesSimulation of strongly correlated fermions in two spatial dimensions with fermionic projected entangled-pair statesFermionic projected entangled pair statesContraction of fermionic operator circuits and the simulation of strongly correlated fermionsThree-band Hubbard model for Na 2 IrO 3 : Topological insulator, zigzag antiferromagnet, and Kitaev-Heisenberg materialGapless spin liquid ground state of the spin- 1 2 J 1 J 2 Heisenberg model on square latticesGlobal phase diagram of a doped Kitaev-Heisenberg modelEnergy dynamics in the Heisenberg-Kitaev spin chainHoneycomb-Lattice Heisenberg-Kitaev Model in a Magnetic Field: Spin Canting, Metamagnetism, and Vortex CrystalsDynamics of the Kitaev-Heisenberg ModelTopological excitations in the ferromagnetic Kitaev-Heisenberg modelSpin liquid fingerprints in the thermal transport of a Kitaev-Heisenberg ladderHeisenberg-Kitaev model in a magnetic field: 1 / S expansionFinite-temperature properties of the Kitaev-Heisenberg models on kagome and triangular lattices studied by improved finite-temperature Lanczos methodsVariational study of the Kitaev-Heisenberg-Gamma modelStudying Two-Dimensional Systems with the Density Matrix Renormalization Group
[1] Anderson P 1973 Mater. Res. Bull. 8 153
[2] Anderson P W 1987 Science 235 1196
[3] Balents L 2010 Nature 464 199
[4] Kitaev A 2003 Ann. Phys. 303 2
[5] Kitaev A 2006 Ann. Phys. 321 2
[6] Choi S K et al. 2012 Phys. Rev. Lett. 108 127204
[7] Singh Y et al. 2012 Phys. Rev. Lett. 108 127203
[8] Plumb K W et al. 2014 Phys. Rev. B 90 041112
[9] Sears J A et al. 2015 Phys. Rev. B 91 144420
[10] Baek S H et al. 2017 Phys. Rev. Lett. 119 037201
[11] Wolter A U B et al. 2017 Phys. Rev. B 96 041405
[12] Sears J A et al. 2017 Phys. Rev. B 95 180411
[13] Gass S et al. 2020 Phys. Rev. B 101 245158
[14] Bachus S et al. 2020 Phys. Rev. Lett. 125 097203
[15] Agrestini S et al. 2017 Phys. Rev. B 96 161107
[16] Chaloupka J et al. 2010 Phys. Rev. Lett. 105 027204
[17] Chaloupka J et al. 2013 Phys. Rev. Lett. 110 097204
[18] Schaffer R et al. 2012 Bhattacharjee S & Kim Y B Phys. Rev. B 86 224417
[19] Reuther J et al. 2011 Phys. Rev. B 84 100406
[20] Sela E et al. 2014 Phys. Rev. B 90 035113
[21] Liu X et al. 2011 Phys. Rev. B 83 220403
[22] Ye F et al. 2012 Phys. Rev. B 85 180403
[23] Ran K J, Wang J H, Bao S et al. 2022 Chin. Phys. Lett. 39 027501
[24] Ji J T, Sun M J, Cai Y Z et al. 2021 Chin. Phys. Lett. 38 047502
[25] Wang L, Zhang Y L, and Sandvik A W 2022 Chin. Phys. Lett. 39 077502
[26] Zhao X X, Ran K J, Wang J H et al. 2022 Chin. Phys. Lett. 39 057501
[27] Zheng J C et al. 2017 Phys. Rev. Lett. 119 227208
[28] Hassan S R et al. 2013 Phys. Rev. Lett. 110 037201
[29] Liang L, Wang Z, and Yu Y 2014 Phys. Rev. B 90 075119
[30] Faye J P L, Sénéchal D, and Hassan S R 2014 Phys. Rev. B 89 115130
[31] Banerjee A et al. 2016 Nat. Mater. 15 733
[32] Mohapatra S and Singh A 2019 J. Magn. Magn. Mater. 479 229
[33] Rau J G, Lee E K H, and Kee H Y 2014 Phys. Rev. Lett. 112 077204
[34] Lubasch M, Cirac J I, and Bañuls M C 2014 New J. Phys. 16 033014
[35] Orús R 2014 Ann. Phys. 349 117
[36] Verstraete F, Murg V, and Cirac J 2008 Adv. Phys. 57 143
[37] Jiang H C, Weng Z Y, and Xiang T 2008 Phys. Rev. Lett. 101 090603
[38] Verstraete F and Cirac J I 2004 arXiv:cond-mat/0407066 [cond-mat.str-el]
[39] Liu W Y, Dong S J, Han Y J, Guo G C, and He L 2017 Phys. Rev. B 95 195154
[40] Dong S J, Wang C, Han Y, Guo G C, and He L 2019 Phys. Rev. B 99 195153
[41] Gu Z C, Verstraete F, and Wen X G 2010 arXiv:1004.2563 [cond-mat.str-el]
[42] Corboz P and Orús R, Bauer B, and Vidal G 2010 Phys. Rev. B 81 165104
[43] Kraus C V and Schuch N, Verstraete F, and Cirac J I 2010 Phys. Rev. A 81 052338
[44] Barthel T, Pineda C, and Eisert J 2009 Phys. Rev. A 80 042333
[45] Laubach M et al. 2017 Phys. Rev. B 96 121110
[46] Liu W Y et al. 2018 Phys. Rev. B 98 241109
[47] Okamoto S 2013 Phys. Rev. B 87 064508
[48] Steinigeweg R and Brenig W 2016 Phys. Rev. B 93 214425
[49] Janssen L, Andrade E C, and Vojta M 2016 Phys. Rev. Lett. 117 277202
[50] Gohlke M and Verresen R, Moessner R, and Pollmann F 2017 Phys. Rev. Lett. 119 157203
[51] Joshi D G 2018 Phys. Rev. B 98 060405
[52] Metavitsiadis A, Psaroudaki C, and Brenig W 2019 Phys. Rev. B 99 205129
[53] Cônsoli P M and Janssen L, Vojta M, and Andrade E C 2020 Phys. Rev. B 102 155134
[54] Morita K and Tohyama T 2020 Phys. Rev. Res. 2 013205
[55] Zhang S S and Halász G B, Zhu W, and Batista C D 2021 Phys. Rev. B 104 014411
[56] Stoudenmire E and White S R 2012 Annu. Rev. Condens. Matter Phys. 3 111