Introduction

The search for quantum spin liquids (QSLs) has been an active endeavor in condensed matter physics for decades1,2,3,4,5,6,7,8. In contrast to magnetically ordered states, which are characterized by spontaneous symmetry breaking and local order parameters, this novel state of matter manifests itself in other exotic ways, such as long-range quantum entanglement and fractional spin excitations, which are represented by emergent particles (spinons) and gauge fields. In particular, a gapped QSL must be associated with a non-trivial topological order9,10. A chiral spin liquid (CSL) is defined as a topological QSL that breaks time-reversal symmetry (TRS) and parity symmetry. It was first proposed by Kalmeyer and Laughlin11,12 that the bosonic fractional quantum Hall state with filling factor ν = 1/2 could be the ground state of some frustrated Heisenberg antiferromagnets in two dimensions. Soon after that, Wen, Wilczek, and Zee13 introduced the spin chirality \({E}_{123}\equiv {{{{\bf{S}}}}}_{1}\cdot \left({{{{\bf{S}}}}}_{2}\times {{{{\bf{S}}}}}_{3}\right)\) to characterize such a TRS-breaking QSL state, which they named CSL.

With strong geometric frustration and concrete quantum material realizations, the S = 1/2 Kagome Heisenberg antiferromagnet (KHAF) has been regarded as one of the promising candidates to host QSLs14,15,16,17,18. Theoretically, it has been generally accepted that the ground state of the KHAF with nearest-neighbor (NN) interactions is a QSL state6,18, although a valence-bond crystal (VBC) with a 36-site unit cell has also been proposed19,20,21,22. However, despite extensive research23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47, the nature of the QSL state is still a mystery. In particular, the issue of “to gap or not to gap” has not yet been settled down unambiguously. While exact diagonalization24,25,26,27,28,29,30, density matrix renormalization group (DMRG)32,33,34,35,36, and a tensor network variational calculation37 indicate a finite energy gap, variational Monte Carlo38,39,40,41, infinite DMRG42, and further tensor network studies43,44 favor a gapless ground state, i.e., the U(1) Dirac QSL38.

Meanwhile, some early mean-field theory calculations proposed the possibility of a CSL19,48,49,50. Further studies have shown that a CSL can be stabilized in the KHAF with chiral three-spin interactions51, or with 2nd and 3rd NN Ising52 or Heisenberg couplings53,54,55,56, where the nature of this CSL was found to be of the ν = 1/2 Kalmeyer-Laughlin type.

In this work, we reexamine the evolution of the ground state of the KHAF with equal 2-nd and 3-rd NN Heisenberg couplings. We use the DMRG method57,58,59 (up to bond dimension D = 18,000) and analytical analyses to map out its ground-state phase diagram (Fig. 1a), which consists of a \(\sqrt{3}\times \sqrt{3}\) ordered phase, two distinct VBC phases, and a CSL phase. The most striking observation is that the ground state of the KHAF with only NN interactions lies in the aforementioned CSL phase, which is supported by the calculations of energy susceptibility, wave function fidelity, spin chirality, ground-state degeneracy, and the existence of topologically degenerate ground states in the semion sector. It is further shown that the CSL phase extends into the region where the 2nd and 3rd NN couplings are ferromagnetic.

Fig. 1: Basic description of the model.
figure 1

a The phase diagram of the \(J\,{{\mbox{-}}}\,{J}^{{\prime} }\) model consisting of a \(\sqrt{3}\times \sqrt{3}\) ordered phase, two VBC phases (VBC I and VBC II with different patterns, see Fig. 4e and f), and a CSL phase with finite spin chirality χ. Note that the KHAF at \({J}^{{\prime} }=0\) (red star) also lies in the CSL phase with χ ≈ 2 × 10−4. b YC8 cylinder for the Kagome lattice, where the dashed red line represents periodic boundary condition, and the flat right edge is highlighted by the solid red line. The colored area denotes a unit cell of the \([\frac{\pi }{2},0]\) mean-field ansatz, with π/2 (zero) flux through the elementary triangles (hexagons). c Entanglement spectra of the left half of the cylinder, labeled by \({K}_{y}^{L}\) (momentum along the y-direction) and \({S}_{z}^{L}\) (total z-component spin), for ground states in the (upper panel) identity and (lower panel) semion sectors on the 152-site YC8 cylinder with \({J}^{{\prime} }=0.4\).

Full size image

Results

Model

We consider the KHAF with SU(2)-symmetric longer-range couplings (referred to as the \(J\,{{\mbox{-}}}\,{J}^{{\prime} }\) model hereafter)

$$H=J\mathop{\sum}\limits_{{\langle ij\rangle }_{1}}{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+{J}^{{\prime} }\left(\mathop{\sum}\limits_{{\langle ij\rangle }_{2}}{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}+\mathop{\sum}\limits_{{\langle ij\rangle }_{3}}{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}\right)\,,$$
(1)

where 〈ij〉γ (γ = 1, 2, 3) denote the γ-th NN bonds on the Kagome lattice (the 3rd NN bonds are defined within hexagons only). While it is also interesting to explore the general phase diagram when the 2-nd and 3-rd NN bond couplings are not equal52,53,54, we restrict ourselves to the model in Eq. (1), which, as we shall see below, proves sufficient for understanding the interpolation to the KHAF at \({J}^{{\prime} }=0\). Henceforth we set J = 1 as the energy unit and focus on the regime of \(-0.15\le {J}^{{\prime} }\le 0.4\) in the YC cylinder geometry (Fig. 1b).

Chiral spin liquid phase

Let us start the discussion of the Hamiltonian in Eq. (1) from its relatively well-understood CSL phase with moderately large \({J}^{{\prime} }\gtrsim 0.08\)53,54, for which much insight can be gained from a Gutzwiller projected wave function55,56. This wave function is constructed in a fermionic parton representation for S = 1/2 spins, \({S}_{i}^{a}=\frac{1}{2}{{{{\boldsymbol{c}}}}}_{i}^{{\dagger} }{\sigma }^{a}{{{{\boldsymbol{c}}}}}_{i}\) (a = x, y, z), where σa are Pauli matrices and \({{{{\boldsymbol{c}}}}}_{i}={({c}_{i,\uparrow },{c}_{i,\downarrow })}^{T}\) are fermion annihilation operators. The physical Hilbert space of S = 1/2 is restored by imposing a local single-occupancy constraint \({\sum }_{\sigma = \uparrow ,\downarrow }{c}_{i\sigma }^{{\dagger} }{c}_{i\sigma }=1\).

In the parton representation, a mean-field Hamiltonian, parameterized by link variables \({\chi }_{ij}={e}^{{\phi }_{ij}}\) on the 1st NN bonds, takes the form of38

$${H}_{{{{\rm{MF}}}}}=\mathop{\sum}\limits_{{\langle ij\rangle }_{1}}\mathop{\sum}\limits_{\sigma =\uparrow ,\downarrow }({\chi }_{ij}{c}_{i\sigma }^{{\dagger} }{c}_{j\sigma }+\,{{\mbox{H.c.}}}\,)\,.$$
(2)

Due to the SU(2) gauge redundancy9,60,61,62, Gutzwiller projected states obtained from Eq. (2) are distinguished by the gauge fluxes threading through elementary triangles and hexagons on the Kagome lattice (Fig. 1b). A CSL state is of particular interest that corresponds to the ansatz with \(\frac{\pi }{2}\)-flux through triangles and zero-flux through hexagons, dubbed as \([\frac{\pi }{2},0]\) state. Note that other parton Ansätze in this work will be labeled in the same scheme.

With a cylindrical boundary condition, the \([\frac{\pi }{2},0]\) state exhibits four exact boundary zero modes, \({d}_{L\sigma }^{{\dagger} }\) and \({d}_{R\sigma }^{{\dagger} }\), where L (R) denotes the left (right) boundary of the cylinder. For the CSL state, the minimally entangled states (anyon eigenbasis)63 are constructed as follows64,65:

$$\left\vert {\Psi }_{1}\right\rangle ={\hat{P}}_{G}{d}_{L\uparrow }^{{\dagger} }{d}_{L\downarrow }^{{\dagger} }\left\vert \Phi \right\rangle ,\quad \left\vert {\Psi }_{2}\right\rangle ={\hat{P}}_{G}{d}_{L\uparrow }^{{\dagger} }{d}_{R\downarrow }^{{\dagger} }\left\vert \Phi \right\rangle \,,$$
(3)

where \(\left\vert \Phi \right\rangle\) is a Fermi sea state with all negative energy modes of Eq. (2) being occupied, and \({\hat{P}}_{G}\) is the Gutzwiller projection operator imposing the single-occupancy constraint. With the MPO-MPS method65,66, the CSL anyon eigenbasis in Eq. (3) can be converted into matrix product states (MPSs) with high fidelity.

For the \(J\,{{\mbox{-}}}\,{J}^{{\prime} }\) model in the CSL phase, two topologically degenerate ground states, denoted by \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) (\(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)) in the identity (semion) sector, can be obtained by initializing DMRG with their parton counterpart \(\left\vert {\Psi }_{1}\right\rangle\) (\(\left\vert {\Psi }_{2}\right\rangle\))67,68. For a YC8 cylinder with 152 sites, the total ground-state energy difference (including boundary contributions) between two topological sectors is ΔE0 ≈ 0.33(0.27) for \({J}^{{\prime} }=0.4\,(0.2)\), and the two ground states are orthogonal to each other, as indicated by their overlap ∣〈ΨI∣ΨS〉∣~10−11(10−9). We note that the time-reversal partner of \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) (\(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)), denoted by \(\left\vert {{\Psi }_{{{{\rm{I}}}}}}^{* }\right\rangle\) (\(\left\vert {{\Psi }_{{{{\rm{S}}}}}}^{* }\right\rangle\)), has the same energy as \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) (\(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)) due to the spontaneous TRS breaking. As shown in Fig. 1c, the entanglement spectrum69,70 of \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) (\(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)) displays the characteristic counting {1, 3, 4, 7, . . . } ({2, 2, 6, 8, . . . }), which suggests the chiral SU(2)1 Wess-Zumino-Witten model being the edge conformal field theory. These results firmly validate the existence of a CSL phase of Kalmeyer-Laughlin type in this model.

Spontaneous TRS breaking at \({J}^{{\prime} }=0\)

After examining the CSL phase at relatively large \({J}^{{\prime} }\), we now show that this CSL phase is adiabatically connected to the ground state at \({J}^{{\prime} }=0\) point. For this purpose, we carry out DMRG calculations that are initialized with randomly generated MPSs (dubbed as Random-DMRG) as well as specific wave functions (dubbed as Boosted-DMRG67). For Boosted-DMRG, we adopt either the parton ansatz (3) or post-optimized MPSs of the neighborhood Hamiltonian in parameter space (see also a related adiabatic DMRG approach53,71). Nevertheless, the ground-state energies obtained by Random- and Boosted-DMRG are almost identical (see Method), both of which agree with the results in ref. 54. For the KHAF with \({J}^{{\prime} }=0\), the ground-state energy per site is found to be −0.4383(6), which also agrees with the best available DMRG results33,34.

For \(-0.15\le {J}^{{\prime} }\le 0.4\), the first-order derivative of the ground-state energy E0 does not show any singularities, as shown in Fig. 2a. However, we observe that the second-order derivative of E0 exhibits three sharp peaks (Fig. 2a). The location of these peaks indicates that there is no (at least neither first nor second-order) phase transition within the region of \(0\le {J}^{{\prime} }\le 0.1\). We emphasize that these results are stable against the system-size and bond-dimension scaling (see Supplementary Note 5) and indeed, are consistent with a previous DMRG result54. In contrast, previous DMRG studies53,54 suggest a critical point at \({J}_{c}^{{\prime} }\approx 0.08\) that is identified by vanishing spin chirality. This contradiction can be resolved as follows (see method and Supplementary Note 1): Consider the spin chirality χijk(≥0) defined on an elementary triangle Δijk (Fig. 1b), we find that (i) χijk is sizable for a state deep inside the CSL phase, e.g., χijk ~ 0.06 at \({J}^{{\prime} }=0.2\); (ii) χijk will decrease rapidly as \({J}^{{\prime} }\) decreases and exceed \({J}^{{\prime} }=0.1\), e.g., χijk ~ 6 × 10−4 (2 × 10−4) at \({J}^{{\prime} }=0.08\,(0)\), which is small but does not vanish numerically (note that the truncation error in our DMRG calculations is ≲10−6), and is stable against finite-size scaling (see Method and Supplementary Note 1); (iii) χijk continues to decrease to a numerical zero, χijk ~ 6 × 10−7 at \({J}^{{\prime} }=-0.04\), around which CSL to VBC II phase transition occurs (Figs. 1a and 2). Moreover, in Fig. 3, we show the real-space chiral-chiral correlation function at \({J}^{{\prime} }=0\). Comparing \({J}^{{\prime} }=0\) and \({J}^{{\prime} }=0.2\) (see the results of \({J}^{{\prime} }=0.2\) in Supplementary Note 1 and also in ref. 54) data with various parameters, one can conclude that (i) the long-ranged correlation of spin chirality becomes evident only when the circumference and bond dimensions are sufficiently large. This phenomenon holds true even when the system is deep inside the CSL phase, namely, \({J}^{{\prime} }=0.2\); (ii) Increasing the circumference of cylinders and/or bond dimensions can enhance the chiral-chiral correlation obtained by DMRG, thereby amplifying the strength of spin chirality. These observations strongly suggest that even at \({J}^{{\prime} }=0\), the chiral-chiral correlation tends to be long-ranged on cylinders with sufficiently large circumferences, provided that the bond dimensions are adequately sized to catch this feature.

Fig. 2: Evidence for Kramers’ degeneracy and topological degeneracy.
figure 2

a The first- (blue dots) and second-order (red triangles) derivatives of the ground-state energy E0 versus \({J}^{{\prime} }\) on the 224-site YC8 cylinders. The inset shows a zoom-in of the second-order energy derivatives around the \({J}^{{\prime} }=0\) point. b The neighboring wave function overlap \(\langle \Psi ({J}^{{\prime} })| \Psi ({J}^{{\prime} }-\Delta {J}^{{\prime} })\rangle\) (\(\Delta {J}^{{\prime} }=0.01\)) and the corresponding susceptibility versus \({J}^{{\prime} }\). The empty circles are obtained by Random-DMRG, and the full squares are obtained by Boosted-DMRG. Two states at the same \({J}^{{\prime} }\) have quasi-degenerated energies (see Method). c Neighboring wave function overlaps (green squares) and energy variances (blue triangles) of semion-sector ground states in 296-site YC8 cylinders obtained by Boosted-DMRG. d The low-lying entanglement spectra of semion-sector ground states remain qualitatively similar as \({J}^{{\prime} }\) decreases.

Full size image
Fig. 3: Chiral-chiral correlation functions.
figure 3

Semi-Log plot of chiral-chiral correlation function as a function of distance r for the \({J}^{{\prime} }=0\) model on YC8 and YC12 cylinders. Inset: Zoom-in of the area bounded by gray box with a linear plot.

Full size image

In addition to the energetics, we have also calculated the wave function fidelity \(F({J}^{{\prime} })=| \langle \Psi ({J}^{{\prime} })| \Psi ({J}^{{\prime} }-\Delta {J}^{{\prime} })\rangle |\) and the corresponding susceptibility \({\chi }_{F}\equiv \,{{\mbox{d}}}F/{{\mbox{d}}}{J}^{{\prime} }\), where \(\left\vert \Psi ({J}^{{\prime} })\right\rangle\) is the ground state obtained by Random- or Boosted-DMRG initialized with \(\left\vert \Psi ({J}^{{\prime} }+\Delta {J}^{{\prime} })\right\rangle\). Here we adopt \(\Delta {J}^{{\prime} }=0.01\). For \({J}^{{\prime} }\le -0.04\), the kinks in F and singularities in χF, as shown in Fig. 2b, are in agreement with the location of the sharp peaks in the second-order derivative of E0.

For \(-0.03\, <\, {J}^{{\prime} }\le 0.09\), the fidelity F obtained from Random-DMRG (using real MPSs) fluctuates strongly without any regular pattern. This seemingly odd result indicates that the system either has degenerate ground states (from which DMRG chooses their superposition, which varies randomly at different \({J}^{{\prime} }\)) or hosts a pile of phase transitions within this narrow region. As we mentioned earlier, \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) and its time-reversal partner \(\left\vert {{\Psi }_{{{{\rm{I}}}}}}^{* }\right\rangle\) are Kramers degenerate ground states in the CSL phase, whose real combinations are \(\left\vert {\Phi }_{{{{\rm{I}}}}}\right\rangle =(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle +\left\vert {{\Psi }_{{{{\rm{I}}}}}}^{* }\right\rangle )/\sqrt{2}\) and \(\left\vert {\tilde{\Phi }}_{{{{\rm{I}}}}}\right\rangle =i(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle -\left\vert {{\Psi }_{{{{\rm{I}}}}}}^{* }\right\rangle )/\sqrt{2}\). As usual, Random-DMRG using real MPSs would find their superposition \(\cos (\theta )\left\vert {\Phi }_{{{{\rm{I}}}}}\right\rangle +\sin (\theta )\left\vert {\tilde{\Phi }}_{{{{\rm{I}}}}}\right\rangle\). Though DMRG is typically biased towards a low-entanglement combination, here the state is an equal weight superposition of two minimally entangled states, i.e., \(({e}^{i\theta }\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle +{e}^{-i\theta }\left\vert {{\Psi }_{{{{\rm{I}}}}}}^{* }\right\rangle )/\sqrt{2}\). Therefore, DMRG has no preference over a particular θ which is not under control and varies with \({J}^{{\prime} }\). This would lead to a fluctuating fidelity F in Fig. 2b. Furthermore, we note that θ can be fixed (or slowly varying) by performing a series of Boosted-DMRG calculations, i.e., the ground-state search for coupling \({J}^{{\prime} }\) is carried out by initializing DMRG with the post-optimized MPS for coupling \({J}^{{\prime} }+\Delta {J}^{{\prime} }\). This adiabatic procedure can be done recursively from \({J}^{{\prime} }=0.1\) to \({J}^{{\prime} }=-0.03\) and, indeed, the fidelity F now becomes a plateau close to 1, as shown in Fig. 2b. This is sharp evidence for the existence of degenerate ground states and also rules out the multiple-transition scenario.

Generally, the two MPSs obtained from Random-DMRG and Boosted-DMRG are not orthogonal to each other and can be used to determine the two-dimensional ground-state subspace spanned by \(\left\vert {\Phi }_{{{{\rm{I}}}}}\right\rangle\) and \(\left\vert {\tilde{\Phi }}_{{{{\rm{I}}}}}\right\rangle\). This can be formulated as a variational problem by diagonalizing the Hamiltonian in this subspace, which yields two orthogonal states and their energies. For \(-0.03\le {J}^{{\prime} }\le 0.08\), the relative difference between these two energies is smaller than 3 × 10−5 (see method). Also, we compute the energy variance σ2 ≡ 〈H2〉 − 〈H〉2 with respect to DMRG-optimized ground states and find that σ2/N is smaller than 2 × 10−4. Combining these observations with quasi-degenerate energy (see method), we conclude that the states obtained by Random- and Boosted-DMRG indeed form a two-dimensional ground-state space. It is also worth mentioning that the energy difference between these two states is smaller than the gap estimates in previous studies30,33.

Topologically degenerate semion sector

Another essential signal of CSL around \({J}^{{\prime} }=0\) is that \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\), the ground state in the semion sector on a 296-site YC8 cylinder, can be continuously evolved from \({J}^{{\prime} }=0.2\) to \({J}^{{\prime} }=0\) by using Boosted-DMRG. As shown in Fig. 2c, the wave function fidelity F between two neighboring-\({J}^{{\prime} }\) \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)’s is always close to unity for \(0.02\le {J}^{{\prime} }\le 0.2\) and then drops to a smaller but finite value (F ≳ 0.4). Note that this sudden drop is caused by a slight shift of edge semions in real space and does not correspond to a phase transition, as will be discussed in next paragraph. Furthermore, the entanglement spectra of \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)’s also show qualitative consistency in the region \(0\le {J}^{{\prime} }\le 0.2\), as demonstrated in Fig. 2d. Moreover, the tiny per-site energy variances of \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\)’s, σ2/N ~ 10−4, suggests that \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\) is an eigenstate of Hamiltonian (1). A study of the scaling behavior of σ2 further confirms the robustness of our results (for details, see Supplementary Note 5).

We would like to delve into the sudden drop of the wave function fidelity starting at \({J}^{{\prime} }\approx 0.01\). It turns out that when \({J}^{{\prime} }\) is closer and closer to \({J}^{{\prime} }=0\) (or be precise, \({J}^{{\prime} }={J}_{c}\approx -0.04\)), the finite-size energy gap between identity and semion sectors becomes smaller and smaller. Meanwhile, the edge semions at each boundary are more and more delocalized (see Supplementary Note 2), which increases the likelihood of these semions “colliding” with each other on finite cylinders. Thereby, it is more difficult to stabilize \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\) as a “metastable” higher-energy state in DMRG, and the collapse of \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\) into the identity-sector ground state \(\left\vert {\Psi }_{{{{\rm{I}}}}}\right\rangle\) is more likely to occur due to the DMRG energy optimization. This agrees with our observation that the evolution of \(\left\vert {\Psi }_{{{{\rm{S}}}}}\right\rangle\) can be unstable when the length of the cylinder is not long enough. In fact, the evolution breaks down at around \({J}^{{\prime} } < 0.02\) when examining a shorter 152-site YC8 cylinder (see Supplementary Note 2).

For the most contentious KHAF with \({J}^{{\prime} }=0\), we have also considered two more Gutzwiller Ansätze associated with Eq. (2): (i) [0, π] state [a gapless U(1) Dirac QSL] and (ii) [0, 0] state (a QSL with spinon Fermi surface)38. After converting them into MPSs, we find that their overlaps with the DMRG-optimized ground states are almost zero. When using these two Ansätze for initializing DMRG, neither of them can improve the performance of DMRG (see Supplementary Note 3). Although this does not rule out the U(1) Dirac and spinon Fermi surface QSL scenarios, it gives a hint that neither [0, π] nor [0, 0] ansatz is a good description for the ground state of the KHAF.

Phase diagram

To investigate the nature of the other three phases for \({J}^{{\prime} }\lesssim -0.04\), we calculate the NN spin-spin correlation 〈Si ⋅ Sj〉 and the spin structure factor

$$S({{{\bf{q}}}})=\frac{1}{N}\mathop{\sum}\limits_{ij}\langle {{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}\rangle {e}^{i{{{\bf{q}}}}\cdot ({{{{\bf{r}}}}}_{i}-{{{{\bf{r}}}}}_{j})}\,,$$
(4)

which are illustrated in Fig. 4. For \({J}^{{\prime} } < -0.11\), S(q) exhibits clear peaks at the K points in the extended Brillouin zone (Fig. 4a), and the NN correlation 〈Si ⋅ Sj〉 is spatially uniform (not shown), which can be identified as the \(\sqrt{3}\times \sqrt{3}\) spin-ordered phase23,31,72. In the two intermediate phases (\(-0.11\le {J}^{{\prime} }\le -0.04\)), S(q) are featureless (Fig. 4b, c), indicating the absence of long-range magnetic order, and the NN correlation 〈Si ⋅ Sj〉 exhibits clear translational symmetry breaking patterns (Fig. 4e, f), which are two different VBC phases (see more details in Supplementary Note 4) and dubbed as VBC I and VBC II in Fig. 1a. Finally, in the CSL phase, S(q) is featureless (Fig. 4d) and the NN correlation 〈Si ⋅ Sj〉 is also spatially uniform (not shown), which are consistent with a QSL state.

Fig. 4: Spin correlation functions.
figure 4

Top: The spin structure factor S(q) on 224-site YC8 cylinders for a the \(\sqrt{3}\times \sqrt{3}\) phase, b the VBC I phase, c the VBC II phase, and d the CSL phase. The solid (dashed) hexagon is the first (extended) Brillouin zone of the Kagome lattice. Bottom: Normalized absolute value of the NN spin-spin correlation 〈Si ⋅ Sj〉 on a 224-site YC8 cylinder for e the VBC I phase and f the VBC II phase. Note that the VBC II phase is similar with the VBC phase in ref. 54. The cylinders are shown with only central columns.

Full size image

Discussion

To summarize, we revisit the KHAF by combining the DMRG method and analytical analyses and provide clear evidence that the CSL phase in the KHAF with longer-range interactions extends to a broader region in the phase diagram and, most strikingly, includes the KHAF with only NN interactions. In the vicinity of the CSL phase, two distinct VBC phases and a \(\sqrt{3}\times \sqrt{3}\) magnetically ordered phase emerge in the phase diagram (in order of decreasing \({J}^{{\prime} }\)).

Moreover, Gutzwiller projected wave function provides us a good initial ansatz to systematically prepare and study ground states in distinct topological sectors, without involving other technique tricks, such as pinning fields and flux insertion53,73. Thanks to this great advantage, the ground state in the semion sector can be properly prepared, and then adiabatically evolved by DMRG all the way up to \({J}^{{\prime} }=0\). Since the semion sector is a fingerprint of a Kalmeyer-Laughlin-type CSL, we argue that a gapped CSL appears to be the most promising candidate for the ground state at \({J}^{{\prime} }=0\).

Our work may also attract new future attention to this pendent issue. It would be desirable to observe the degenerate ground states at \({J}^{{\prime} }=0\) by, e.g., larger-scale DMRG calculations on wider cylinders or tensor network calculations in the thermodynamic limit. From the experimental side, the signature of semion excitations could be explored in real materials, even when the Kagome compounds, having interactions beyond the KHAF with NN interactions, may no longer have a CSL as their ground states.

Methods

Energetic comparisons and re-diagonalization

We introduce the relative energy difference between two close energies, E1 and E2, as

$$\delta ({E}_{1},{E}_{2})\equiv \left\vert \frac{{E}_{1}-{E}_{2}}{({E}_{1}+{E}_{2})/2}\right\vert ,$$

to quantify their closeness. As mentioned above, we obtained two almost degenerate ground states, \(\left\vert {\Psi }_{{{{\rm{RD}}}}}\right\rangle\) and \(\left\vert {\Psi }_{{{{\rm{BD}}}}}\right\rangle\), by using Random-DMRG and Boosted-DMRG, respectively. In the CSL phase, these two states are generally not orthogonal to each other. Denoting the variational energy of \(\left\vert {\Psi }_{{{{\rm{RD}}}}}\right\rangle\) (\(\left\vert {\Psi }_{{{{\rm{BD}}}}}\right\rangle\)) as ERD (EBD), the relative energy difference δ(ERD, EBD) is extremely small. Namely, for \(-0.03\le {J}^{{\prime} }\le 0.08\), the relative difference between these two energies is smaller than 2 × 10−6 (see Supplementary Note 6).

Meanwhile, with two non-orthogonal MPSs \(\left\vert {\Psi }_{{{{\rm{RD}}}}}\right\rangle\) and \(\left\vert {\Psi }_{{{{\rm{BD}}}}}\right\rangle\) in hand, we carry out a further variational optimization to obtain re-orthogonalized states and corresponding energies. This optimization is formulated as a generalized eigenvalue problem by constructing the following 2 × 2 matrices:

$${H}_{2}=\left(\begin{array}{cc}\langle {\Psi }_{{{{\rm{RD}}}}}| H| {\Psi }_{{{{\rm{RD}}}}}\rangle &\langle {\Psi }_{{{{\rm{RD}}}}}| H| {\Psi }_{{{{\rm{BD}}}}}\rangle \\ \langle {\Psi }_{{{{\rm{BD}}}}}| H| {\Psi }_{{{{\rm{RD}}}}}\rangle &\langle {\Psi }_{{{{\rm{BD}}}}}| H| {\Psi }_{{{{\rm{BD}}}}}\rangle \end{array}\right),$$
(5)

and

$${N}_{2}=\left(\begin{array}{cc}\langle {\Psi }_{{{{\rm{RD}}}}}| {\Psi }_{{{{\rm{RD}}}}}\rangle &\langle {\Psi }_{{{{\rm{RD}}}}}| {\Psi }_{{{{\rm{BD}}}}}\rangle \\ \langle {\Psi }_{{{{\rm{BD}}}}}| {\Psi }_{{{{\rm{RD}}}}}\rangle &\langle {\Psi }_{{{{\rm{BD}}}}}| {\Psi }_{{{{\rm{BD}}}}}\rangle \end{array}\right).$$
(6)

By solving the generalized eigenvalue problem, namely, H2x = EN2x, we can obtain two re-diagonalized eigen-energies, \({\tilde{E}}_{{{{\rm{1}}}}}\) and \({\tilde{E}}_{{{{\rm{2}}}}}\). The relative energy difference \(\delta ({\tilde{E}}_{{{{\rm{1}}}}},{\tilde{E}}_{{{{\rm{2}}}}})\), with respect to newly obtained \({\tilde{E}}_{{{{\rm{1}}}}}\) and \({\tilde{E}}_{{{{\rm{2}}}}}\), can be also found in Supplementary Note 6. We can see that these energies are still very close to each other (for \(-0.03\le {J}^{{\prime} }\le 0.08\), the relative difference between these two energies is smaller than 3 × 10−5), implying that \(\left\vert {\Psi }_{{{{\rm{RD}}}}}\right\rangle\) and \(\left\vert {\Psi }_{{{{\rm{BD}}}}}\right\rangle\) are indeed the two states in the ground state manifold.

Spin chirality order and correlation functions

In order to characterize the spin chirality, we introduce the three-spin operator on an elementary triangle Δijk (see the elementary triangles in Fig. 1b) as

$${E}_{ijk}={{{{\bf{S}}}}}_{i}\cdot \left({{{{\bf{S}}}}}_{j}\times {{{{\bf{S}}}}}_{k}\right).$$
(7)

A nonzero spin chirality in the ground state provides direct evidence of the spontaneous breaking of time-reversal symmetry, which is useful for characterizing the chiral spin liquid phase in the \(J-{J}^{{\prime} }\) model. Here, we calculate the spin chirality by the following three independent methods.

SC-Method I: Re-diagonalize two degenerated states

Notice that Eijk is purely imaginary and Hermitian, indicating that 〈ψ∣Eijk∣ψ〉 = 0 for any real state \(\left\vert \psi \right\rangle\). In the CSL phase, we can obtain two orthogonal degenerate ground states, \(\left\vert {\tilde{\Phi }}_{1}\right\rangle\) and \(\left\vert {\tilde{\Phi }}_{2}\right\rangle\), by the re-diagonalization method introduced above. Because both \(\left\vert {\tilde{\Phi }}_{1}\right\rangle\) and \(\left\vert {\tilde{\Phi }}_{2}\right\rangle\) are real wave functions, the expectation values of Eijk with respect to those two states must be zero. Nevertheless, we can evaluate the spin chirality by introducing the following 2 × 2 matrix:

$${\chi }_{2}(ijk)=\left(\begin{array}{cc}0&\langle {\tilde{\Phi }}_{1}| {E}_{ijk}| {\tilde{\Phi }}_{2}\rangle \\ \langle {\tilde{\Phi }}_{2}| {E}_{ijk}| {\tilde{\Phi }}_{1}\rangle &0\end{array}\right).$$
(8)

It is obvious that χ2(ijk) is Hermitian and has two eigenvalues, ± χijk (χijk > 0), corresponding to chiral and anti-chiral ground states, respectively. Hereafter we choose the elementary triangle Δijk in the center of a YC8 cylinder to compute χijk. Notice that there are two types of elementary triangles corresponding to up- and down-triangles and they give rise to almost the same results. To simplify the notation, we drop the subscripts and use χ to represent the spin chirality averaged by measuring typical elementary triangles.

SC-Method II: Extract spin chirality from chiral-chiral correlation functions

Although the scalar chirality order parameter is zero for real states, the spin chirality can be approximately detected by measuring the chiral-chiral correlation function as

$$\chi =\sqrt{| \langle {E}_{{i}_{0}{j}_{0}{k}_{0}}{E}_{{i}_{d}{j}_{d}{k}_{d}}\rangle | }\,,$$
(9)

where i0j0k0 (idjdkd) indicates three sites in the starting (ending) elementary triangle. Therefore, d indicates the distance between these two elementary triangles. In the thermodynamic limit, χ is finite for a time-reversal symmetry-breaking state even if d is large. On the finite-length cylinders, we choose two elementary triangles with a sufficiently large distance. To suppress the boundary effect, we chose the case with the starting elementary triangle being at Nx/4 and the ending elementary triangle 3Nx/4, where Nx is the length of the cylinder. Note that the same method is also employed in ref. 54 as the major approach to determine spin chirality.

SC-Method III: Directly measure spin chirality in complex-number MPS

In principle, the spin chirality χ does not necessarily vanished in a complex-number MPS. Therefore, the spin chirality can be determined by directly calculating χijk = 〈ΨI(S)∣Eijk∣ΨI(S)〉, where \(\vert {\Psi }_{{{{\rm{I(S)}}}}}\rangle\) is the ground state obtained by Boosted-DMRG initialized with a parton ansatz of CSL with [π/2, 0] (see \(\vert {\Psi }_{1(2)}\rangle\) in the main text). Because the Hamiltonian preserves the time-reversal symmetry, in finite DMRG calculations, the time-reversal partner, \(\vert {\Psi }_{{{{\rm{I(S)}}}}}^{* }\rangle\), will be gradually mixed into the state by the numerical fluctuations introduced by the optimization process in DMRG. Therefore, the spin chirality measured using this SC-Method III is smaller than the actual value of spin chirality for a CSL state, hence providing a lower bound.