Nematic chiral spin liquid in a Kitaev magnet under external magnetic field

Quantum spin liquid (QSL) is an exotic quantum state of matter where the spins of electrons remain disordered even at absolute zero temperature [1, 2]. Unlike conventional phases of matter, QSL states do not fit the Landau-Ginzburg-Wilson framework, and exhibit long-range quantum entanglement and fractionalized excitations [3, 4, 5]. The exactly soluble Kitaev honeycomb model plays the milestone role for understanding the physics in QSL. In particular, a non-Abelian chiral spin liquid (CSL) can be induced by breaking time-reversal symmetry in the Kitaev model, which possesses the non-Abelian anyon that has potential application in topological quantum computing [6, 7]. Experimentalists have been searching for the Kitaev QSL candidate materials in various compounds with $4d$ and $5d$ transition metal-based honeycomb lattice. These compounds have strong spin-orbit coupling that may induce the bond-dependent Kitaev interaction [8]. Some examples of these compounds are A₂IrO₃ with A = Li, Na [9], H₃LiIr₂O₆, which is a hydrogen intercalated modification of A₂IrO₃ [10], and $\alpha$-RuCl₃ [11, 12, 13, 14]. Recently, the cobalt-based $3d^{7}$ honeycomb magnets such as A₂Co₂TeO₆ [15, 16], A₃Co₂SbO₆ [17, 18], and BaCo₂(AsO₄)₂ [19, 20] are also thought to be candidates of Kitaev materials.

Among these compounds, $\alpha$-RuCl₃ has been widely recognized as a prominent Kitaev material. Due to the non-Kitaev interactions in the material, the compound exhibits a zigzag magnetic order at low temperature [11, 12]. By suppressing the zigzag order with external magnetic fields [21, 22, 23], in the intermediate fields experimental evidences for a QSL have been reported, e.g., the excitation continuum in the inelastic neutron scattering spectrum [13, 14, 24], the half-quantized thermal Hall conductivity [25, 26, 27] suggesting the presence of gapless Majorana edge modes [28, 29], and the quantum oscillations in the low-temperature longitudinal heat conductivity possibly due to a spinon Fermi surface [30]. Despite these intriguing progresses, nature of this state are still under debate [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46].

Motivated by the exciting experimental progress, there have been a plenty of theoretical studies in an effort to find QSL in related spin models. Amongst others, the $K$-$\Gamma$-$\Gamma^{\prime}$ model with the Kitaev interaction $K$ and off-diagonal couplings $\Gamma,\Gamma^{\prime}$ has been considered as a minimal model to describe a broad class of Kitaev materials including $\alpha$-RuCl₃ [47, 48, 49, 50, 51, 52, 53]. Besides the well known Kitaev spin liquids, other promising QSL candidates have also been proposed in extended Kitaev models, including the QSL in the antiferromagnetic (AFM) Kitaev model under an intermediate magnetic field [54, 55, 56, 57, 58, 59, 60], different QSL phases in the $K$-$\Gamma$ model with or without a magnetic field [61, 62, 63], the multinode QSL [64] and nematic QSL in the $K$-$\Gamma$-$\Gamma^{\prime}$ model [65, 66, 67]. However, it is hard to accurately estimate the microscopic parameters that can match the experimental data and thus it is unclear how to connect the found QSLs in model studies to the compounds in experiment.

Recently, a $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model has been proposed for $\alpha$-RuCl₃. Multiple accurate many-body simulations have been combined to examine this model, which reproduces the major experimental features in equilibrium and dynamical measurements, giving the dominant ferromagnetic Kitaev interaction, additional Heisenberg coupling $J$, and other parameters [68]. Interestingly, this model under the $[111]$-direction magnetic field shows an emergent QSL phase between the low-field zigzag order and high-field polarized phase [68]. While the variational Monte Calro (VMC) study proposes a gapped CSL, the exact diagonalization results seem to suggest a gapless state [68]. Remarkably, the theoretically predicted two-transition scenario in Ref. [68] has been witnessed by recent high magnetic field experiment on $\alpha$-RuCl₃ at $\sim 35$ T and $\sim 83$ T, respectively, with the field along the out-of-plane $c^{*}$ axis [69]. A field-induced intermediate phase is found experimentally in the field-angle phase diagram [69]. The experimental findings and the good agreements with model study open a new route to search for QSL in $\alpha$-RuCl₃ and call for further understanding on the QSL phase in the $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model.

In this work, we study the nature of the intermediate QSL phase in the $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model using three methods: density matrix renormalization group (DMRG), VMC, and exponential tensor renormalization group (XTRG) for finite-temperature calculations. The DMRG calculations on different cylinder geometries up to the circumference of $6$ unit cells identify the existence of an intermediate phase, with a robust nonzero scalar spin chiral order $\langle\mathbf{S}_{i}\cdot(\mathbf{S}_{j}\times\mathbf{S}_{k})\rangle$ and a lattice nematic order characterizing a $C_{3}$ symmetry breaking. In addition, the algebraic-like spin correlation decay observed in the ground state and the power-law behavior of the specific heat $C_{m}\sim T^{2}$ obtained at low temperature strongly suggest a gapless QSL state.

We further investigate the intermediate phase using the VMC with extended variational parameters, which unveils that allowing the $C_{3}$ symmetry breaking can further lower the variational energy of the previously identified gapped CSL [68] and enlarge the spinon gap, giving a variational energy very close to that obtained in DMRG. Using quantitative analysis, we examine how the mean field dispersion varies and find that the spinon gap is very small compared to the total band width ($\sim 0.02$), which may explain why the spin correlation decay looks like algebraic on our studied cylinders. The spinon density of states which is approximately proportional to energy above the spinon gap can agree with the $C_{m}\sim T^{2}$ behavior at low temperature. The Van Hove singularity at the higher energy $E/|K|\sim 0.4$ also supports the double peak structure found in the specific heat. The gapless or gapped nature of the QSL may require further simulations on larger sizes, but the predicted low-temperature properties above $T\approx 2$K, which are less affected by finite-size effects, could be tested in high-field experiments on $\alpha$-RuCl₃.

This paper is organized as follows. In Sec. II, we introduce the $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model and numerical methods employed in this study. The DMRG simulations are demonstrated in Sec. III, and the XTRG thermodynamics results are presented in Sec. IV. In Sec. V, we show the nematic order obtained in the VMC calculations and compare the results with the DMRG and XTRG results. We present our summary and discussion in Sec. VI.

We study the honeycomb-lattice $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model in an out-of-plane magnetic field along the $[111]$ direction, which has the Hamiltonian defined as

where $\langle i,j\rangle_{\gamma}$ denotes the nearest-neighbor (NN) $\gamma$ bond with $\{\alpha,\beta,\gamma\}$ being $\{x,y,z\}$ under a cyclic permutation, $\mathbf{S}_{i}=\{S_{i}^{x},S_{i}^{y},S_{i}^{z}\}$ are the pseudo spin-$1/2$ operators at the site $i$. $K$ denotes the FM Kitaev interaction, $J$ is the FM Heisenberg interaction, $\Gamma$ and $\Gamma^{\prime}$ are the off-diagonal spin couplings. Following the microscopic model determined in Ref. [68], we choose the Kitaev coupling as the energy unit $K=-1.0$, $J=-0.1$, $\Gamma=0.3$ and $\Gamma^{\prime}=-0.02$. Without external magnetic field, this system has the finite magnetic point group $D_{3d}\times Z_{2}^{T}$, where $Z_{2}^{T}=\{E,T\}$ is the time-reversal symmetry group, and each element in $D_{3d}$ stands for a combination of spin rotation and lattice rotation owing to the spin-orbit coupling. In this paper, we focus on the system under an external magnetic field along the $[111]$ direction, which gives the Zeeman coupling $h\sum_{i,\alpha}S^{\alpha}_{i}$ in Eq. (1).

In DMRG calculation of the ground state [70], we study the system on a cylinder geometry with the periodic boundary condition along the circumference direction ($Y$ direction) and the open boundary condition along the axis direction ($X$ direction). To avoid possible bias from the lattice geometry, we consider two geometries, namely the XC and YC cylinders shown in Figs. 1(a) and 1(b). For the XC cylinder, there is a $\pi/3$ angle between the $Y$ and $X$ direction. For the YC cylinder, the $Y$ direction is perpendicular to the $X$ direction. We denote the system size as XC$L_{y}$-$L_{x}$ and YC$L_{y}$-$L_{x}$, where $L_{y}$ and $L_{x}$ denote the numbers of unit cells along the two directions, respectively. When we show the results on both cylinders together, we use the length along the circumference direction to demonstrate different system sizes, which is denoted as $W_{y}$ with the lattice constant as the unit. In our DMRG simulation, we keep up to $D=2000$ bond states, which allow us to obtain accurate results on the $L_{y}=6$ cylinder with the truncation error at about $10^{-6}$. To reduce the finite-size effect, we also implement the infinite-DMRG (iDMRG) [71] by using the TeNPy library [72] to check our results, where the system size along the $X$ direction can be taken as long enough to eliminate the size effect in this direction. In the iDMRG simulations, we denote the lattice as XC$L_{y}$ and YC$L_{y}$.

To compute the thermodynamic properties, we employ the exponential tensor renormalization group (XTRG) method [73, 74] to simulate the model on the YC cylinder with $L_{y}=4$ and $L_{x}$ up to $12$. The retained bond dimension $D$ is up to $1000$, which ensures the truncation error smaller than $10^{-4}$ down to the lowest temperature $T/|K|=0.008$.

In this section, we demonstrate our DMRG simulation results, focusing on the QSL phase. Since the phase boundaries have been investigated in a recent study [68], we leave our results on the determination of phase transition points in Appendix A, which identify the transitions at $h\simeq 0.1$ and $h\simeq 0.4$ as shown in Fig. 1(c), in consistent with the previous finding [68]. With increased magnetic field, an intermediate QSL emerges between a zigzag magnetic order at low field and a polarized phase at high field.

The three phases have different features of spin structure factor. In Figs. 1(d-f), we show the intra-sublattice structure factor defined as

where the summation is for the sites in the same sublattice. The obtained results of $\tilde{S}({\bf k})$ are the same for the two sublattices of the honeycomb lattice. In the zigzag phase [Fig. 1(d)], the prominent peaks at the ${\bf M}$ points characterize the zigzag order. In the QSL phase [Fig. 1(e)], a dominant peak appears at the ${\bf\Gamma}$ point with a broad peak at a ${\bf M}$ point. In the polarized phase [Fig. 1(f)], the featureless structure factor reflects the very short-range spin correlation. In the following, we focus on the characterization of the QSL phase.

To further explore the nature of the intermediate QSL phase, we turn to study the specific heat $C_{\rm m}$ and thermal entropy $S$. In a previous work by some of the authors [68], a prominent double-peak/hump feature has been revealed in the curve of the specific heat $C_{\rm m}$ versus temperature, for the intermediate field regime, which characterizes two temperature scales $T_{\rm H}$ (the higher-temperature scale) and $T^{\prime\prime}_{\rm L}$ (the lower-temperature scale). Below $T^{\prime\prime}_{\rm L}$, the calculated specific heat data seem to fall into a power-law decay [68].

To clarify this possible power-law behavior, here we improve the XTRG calculation to larger system sizes, namely YC4-$L_{x}$ cylinders ($L_{x}=6,8,10,12$) and down to lower temperature. As shown in Fig. 7, we find that both the position and height of the $C_{\rm m}$ hump at $T^{\prime\prime}_{\rm L}\simeq 0.07$ (under the magnetic field $h=0.34$) are quite robust with growing system size. Below $T^{\prime\prime}_{L}$, the $C_{\rm m}$ curves exhibit a power-law decay $C_{\rm m}\sim T^{\alpha}$ with the power exponent approaching $\alpha\simeq 2$, reflecting the gapless excitations (or with a very small excitation gap) in this low-energy scale.

Besides, we also show the thermal entropy $S$ in the inset of Fig. 7, where the low-temperature data has well converged versus bond dimension $D$ and also follows a power-law decay close to $S\sim T^{2}$. It is noteworthy that there are considerable amount of entropies below $T^{\prime\prime}_{L}$, indicating that the low-temperature states are strongly fluctuating and have large density of states in low-energy excitations.

Finally, we investigate the intermediate phase of the system using VMC, and try to provide a physical picture to help understand part of the above numerical results.

We introduce the fermionic slave-particle representation for spin operator as $S_{i}^{\alpha}=\frac{1}{2}C_{i}^{\dagger}\sigma^{\alpha}C_{i}$ with $C_{i}^{\dagger}=(c_{i\uparrow}^{\dagger},c_{i\downarrow}^{\dagger})$, $\alpha\equiv x,y,z,$ and $\sigma^{\alpha}$ the Pauli matrices. To ensure that the Hilbert space of the fermions faithfully represents that of the physical spin, a constraint for particle number $\hat{N}_{i}\equiv c_{i\uparrow}^{\dagger}c_{i\uparrow}+c_{i\downarrow}^{\dagger}c_{i\downarrow}=1$ is imposed for every site. This complex fermion representation is equivalent to the Majorana representation using $b_{x},b_{y},b_{z},c$ fermions [7]. For the model described by Eq. (1), the most general expression of the mean-field Hamiltonian reads [83]

where $\Psi_{i}=(C_{i},\bar{C}_{i})$, $\bar{C}_{i}=(c_{i\downarrow}^{\dagger},-c_{i\uparrow}^{\dagger})^{T}$, and $R^{\gamma}_{\alpha\beta}=-i(\sigma^{\alpha}+\sigma^{\beta})/\sqrt{2}$ is a rotation matrix. Considering the SU(2) gauge symmetry in this fermionic representation, the particle-number constraint $\hat{N}_{i}=1$ is extended to the SU(2) gauge-invariant three-component form $\operatorname{Tr}(\Psi_{i}\tau\Psi^{\dagger}_{i})=0$, which is implemented by the Lagrangian multipliers $\lambda^{x,y,z}$ with $\tau^{x,y,z}$ being the generators of the SU(2) gauge group. The mean-field matrices $U^{0,1,2,3}_{ji}$ can be expanded with the identity matrix and $\tau^{x,y,z}$. For the studied model, the SU(2) gauge symmetry breaks down to its $Z_{2}$ subgroup named the invariant gauge group (IGG) [84].

While a QSL state preserves the space group symmetry, the symmetry group of its mean-field Hamiltonian is an enlarged one called the projective symmetry group (PSG), whose group elements are the space group operations followed by certain SU(2) gauge transformation. We adopt the Kitaev PSG [85, 86], namely the PSG of the mean-field description of the pure Kitaev model, then the general form of the coefficients $U_{ji}^{(m)}$ can be expressed as

with $\tau^{c}={1\over\sqrt{3}}(\tau^{x}+\tau^{y}+\tau^{z})$. Hence we have the complete set of variational parameters $\{\rho_{a},\rho_{c},\rho_{d},\rho_{f},\eta_{0},\eta_{3},\eta_{5},\eta^{x}_{7},\eta^{z}_{7},\theta,\lambda\}$, in which the subset $\{\rho_{a},\rho_{c},\rho_{d},\rho_{f},\eta_{0},\eta_{3},\eta_{5}\}$ preserves all the symmetries of the original Hamiltonian Eq. (1). To make the VMC compatible with the spacial nematic order, we introduce three additional parameters $\{\eta^{x}_{7},\eta^{z}_{7},\theta\}$ in the mean-field Hamiltonian to break the $C_{3}$ symmetry but preserve a residual mirror symmetry. The mirror plane is parallel to the $z$ bond and perpendicular to the lattice plane. Then Gutzwiller projection is performed to the mean-field ground state to enforce the particle number constraint. Thus we obtain the trial wave function $|\Psi(\{x\})\rangle=P_{G}|\Psi_{\rm mf}(\{x\})\rangle=\sum_{\alpha}f(\alpha)|\alpha\rangle$, where $\{x\}$ denote the variational parameters and $\alpha$ stands for the Ising bases in the many-body Hilbert space. The energy of the trial state $E(\{x\})=\langle\Psi(\{x\})|H|\Psi(\{x\})\rangle/\langle\Psi(\{x\})|\Psi(\{x\})\rangle=\sum_{\alpha}\frac{|f(\alpha)|^{2}}{\sum_{\gamma}|f(\gamma)|^{2}}(\sum_{\beta}\langle\beta|H|\alpha\rangle\frac{f(\beta)*}{f(\alpha)*})$ is computed using Monte Carlo sampling. Then the optimal parameters $\{x\}$ are determined by minimizing the energy $E(\{x\})$. Our VMC calculations are performed on the torus geometry up to the size $9\times 9\times 2$ with $162$ sites.

In the absence of the variational parameters $\{\eta^{x}_{7},\eta^{z}_{7},\theta\}$, an intermediate CSL phase between the zigzag phase and the polarized trivial phase is identified with Chern number $\nu=2$ [68]. This CSL is an Abelian QSL with the quasiparticle excitations $a$, $\bar{a}$, $\epsilon$, $1$, where $1$ denotes the vacuum, $\epsilon$ is the fermion, $a$ and $\bar{a}$ are two types of vortices with topological spin $e^{i\pi/4}$ [7]. These excitations follow the fusion rules $\epsilon\times\epsilon=1$, $a\times\bar{a}=1$, $a\times\epsilon=\bar{a}$, $\bar{a}\times\epsilon=a$, and $\bar{a}\times\bar{a}=\epsilon$ [7]. The gapless edge of this CSL contains $2$ branches of chiral Majorana excitations with total chiral central charge $c_{-}=1$.

By considering the additional parameters $\{\eta^{x}_{7},\eta^{z}_{7},\theta\}$, our VMC calculations still support a CSL state, but with a nematic order which further lowers the variational energy. For example, on the $8\times 8\times 2$ size at the magnetic field $h=0.3$, the variational energy of the CSL state with the $C_{3}$ symmetry is $-0.3691$. By allowing the $C_{3}$ symmetry breaking, the energy is lowered to $-0.3746$ with an emergent nematic order. This optimized variational energy is very close to the energies found in the DMRG simulation, e.g., $-0.3790$ in the bulk of the XC6 cylinder. The details of the VMC identification of the nematic order can be found in Appendix E. Besides the spontaneous nematic order in the new ground state, the obtained low-lying energies also have a drastic change. In Fig. 8, we show the dispersion of the fermion bands in the optimized mean-field Hamiltonian Eq. (6) with $h=0.3$, which indicates that with the $C_{3}$ symmetry breaking a small excitation gap is opened.

The mean field dispersion is plotted for a system with $120\times 120\times 2$ sites. There are totally eight bands due to the two sublattices, two spin components and the particle-hole redundancy. A remarkable feature of the spinon dispersion is that six of the eight bands are quit flat and are close to zero energy, and the rest two bands are very dispersive. The existence of the nearly flat bands is very similar to the mean field spectrum of the pure Kitaev model, where the dispersive band mainly come from the $c$ fermions and the flat bands result from the $b_{x,y,z}$ fermions. Our variational calculation indicates that four parameters are dominating in the intermediate QSL phase, namely $\rho_{a},\rho_{c},\eta_{0}$ and $\eta_{3}$, where $\rho_{a},\rho_{c}$ come from the Kitaev interactions, $\eta_{0}$ comes from the Heisenberg interaction and $\eta_{3}$ comes from the $\Gamma$ interaction.

It will be illustrative to compare the VMC results with those from DMRG and XTRG. The size of the spinon gap with respect to the total band width is $0.02|K|$ (we have set the band width of the spinons to be $2|K|$) which is very small. Practically, it is not easy to detect such a small gap using correlation functions. This qualitatively interprets the algebraic-like spin correlation decay observed in DMRG results. Furthermore, above the energy $0.02|K|$, the fermionic spinons start to show density of states which is proximately proportional to energy in a small energy window. This agrees with the $C_{m}\propto T^{2}$ law of the specific heat at the temperature region $0.02|K|\sim 0.04|K|$. With further increasing energy, a Van Hove singularity appears in the dispersive band at $E\sim 0.5|K|$ which yields a large density of states. Hence a peak (or shoulder) is expected to appear in the specific heat $C_{m}$ at around $0.5|K|$. We notice that the existence of two-temperature scale structure in the specific heat is indeed observed in the XTRG results shown in Fig. 7.

By combining the DMRG, XTRG, and VMC calculations, we investigate the nature of the intermediate QSL phase of a honeycomb $K$-$J$-$\Gamma$-$\Gamma^{\prime}$ model in an out-of-plane magnetic field along the $[111]$ direction, which is sandwiched between a low-field zigzag order phase and a high-field polarized phase.

By using DMRG simulations on different cylinder geometries, we confirm the absent magnetic order and unveil the algebraic-like spin correlation function decay $C(r)\sim 1/r^{K_{s}}$ with small power exponent $K_{s}\lesssim 1.5$. Besides, we identify a robust spin scalar chiral order $\langle\mathbf{S}_{i}\cdot(\mathbf{S}_{j}\times\mathbf{S}_{k})\rangle$ and a lattice nematic order characterizing the $C_{3}$ symmetry breaking. We also extend the XTRG simulation to longer 4-leg cylinders and down to lower temperature. Below the lower-temperature scale $T^{\prime\prime}_{\rm L}\simeq 0.07$, we observe a robust power-law behavior of the specific heat $C_{\rm m}\sim T^{2}$, as well as considerable amount of thermal entropy. By allowing $C_{3}$ symmetry breaking in the VMC calculation, we obtain the same gapped CSL reported in Ref. [68] but with an additional nematic order that can further lower the variational energy.

While our DMRG results also find the chiral order and nematic order that can agree with the nematic CSL found by VMC, the ground-state and thermal measurements suggest a nearly gapless QSL. Yet, the VMC results provide another picture to understand the finite-size simulation results. The small ratio of the gap with respect to the spinon band width ($\sim 0.02$) indicates a large correlation length that may go beyond the resolution of present DMRG calculation. The finite density of states for spinons, which are nearly proportional to energy in a small energy window above the gap, also agrees with the behavior of the specific heat $C_{m}\sim T^{2}$. Identifying the nature of gapped or gapless may need further studies on larger system size.

It is noteworthy that the QSL predicted here for $\alpha$-RuCl₃ with relatively large low-energy excitation could be further studied experimentally in high-field NMR [87], ESR spectroscopy [88], and specific heat measurements. Besides, this intermediate QSL may also relate to other recently proposed field-induced states in the Kitaev candidates, e.g., the broad magnetic continuum observed by time-domain terahertz spectroscopy measurements on BaCo₂(AsO₄)₂ [20].

We also notice that this QSL phase might be related to the nematic QSL-like phase in the $K$-$\Gamma$-$\Gamma^{\prime}$ model with dominant FM Kitaev interaction and off-diagonal exchanges $\Gamma>0$, $\Gamma^{\prime}<0$ [65, 66, 67]. While this phase was found not connected to the Kitaev spin liquid, many-body simulations found the $C_{3}$ symmetry breaking [66, 67] and gapless-like excitations [67]. It would be interesting to explore the potential connection of the two phases and whether the phase in the $K$-$\Gamma$-$\Gamma^{\prime}$ model may also be a nematic CSL.