Possible chiral spin liquid state in the S = 1 / 2 kagome Heisenberg model

,

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 order [9,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 Laughlin [11,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 Zee [13] introduced the spin chirality E 123 ≡ S 1 • (S 2 × S 3 ) to characterize such a TRS-breaking QSL state, which they named CSL.
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 method [57][58][59] (up to bond dimension D = 18000) and analytical analyses to map out its ground-state phase diagram [see Fig. 1 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.

RESULT Model
We consider the KHAF with SU(2)-symmetric longerrange couplings (referred to as the J-J ′ model hereafter) where ⟨i j⟩ γ (γ = 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 equal [52][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 ′ = 0. Henceforth we set J = 1 as the energy unit and focus on the regime of −0.15 ≤ J ′ ≤ 0.4 in the YC cylinder geometry [see Fig. 1(b)].

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 ′ ≳ 0.08 [53,54], for which much insight can be gained from a Gutzwiller projected wave function [55,56].This wave function is constructed in a fermionic parton representation for S = 1/2 spins, S a i = 1 2 c † i σ a c i (a = x, y, z), where σ a are Pauli matrices and c i = (c i,↑ , c i,↓ ) T are fermion annihilation operators.The physical Hilbert space of S = 1/2 is restored by imposing a local single-occupancy constraint σ=↑,↓ c † iσ c iσ = 1.In the parton representation, a mean-field Hamiltonian, parameterized by link variables χ i j = e ϕ i j on the 1st NN bonds, takes the form of [38] Due to the SU(2) gauge redundancy [60][61][62][63], Gutzwiller projected states obtained from Eq. ( 2) are distinguished by the gauge fluxes threading through elementary triangles and hexagons on the Kagome lattice [see Fig. Rσ , where L (R) denotes the left (right) boundary of the cylinder.For the CSL state, the minimally entangled states (anyon eigenbasis) [64] are constructed as follows [65,66]: where |Φ⟩ is a Fermi sea state with all negative energy modes of Eq. ( 2) being occupied, and PG is the Gutzwiller projection operator imposing the single-occupancy constraint.With the MPO-MPS method [66,67], the CSL anyon eigenbasis in Eq. ( 3) can be converted into matrix product states (MPSs) with high fidelity.
For the J-J ′ model in the CSL phase, two topologically degenerate ground states, denoted by |Ψ I ⟩ (|Ψ S ⟩) in the identity (semion) sector, can be obtained by initializing DMRG with their parton counterpart |Ψ 1 ⟩ (|Ψ 2 ⟩) [68,69].For a YC8 cylinder with 152 sites, the total ground-state energy difference (including boundary contributions) between two topological sectors is ∆E 0 ≈ 0.33 (0.27) for J ′ = 0.4 (0. Spontaneous TRS breaking at J ′ = 0 After examining the CSL phase at relatively large J ′ , we now show that this CSL phase is adiabatically connected to the ground state at J ′ = 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-DMRG [68]).For Boosted-DMRG, we adopt either the parton ansatz (3) or postoptimized MPSs of the neighborhood Hamiltonian in parameter space (see also a related adiabatic DMRG approach [53,72]).Nevertheless, the ground-state energies obtained by Random-and Boosted-DMRG are almost identical (see Supplementary Note 1), both of which agree with the results in Ref. [54].For the KHAF with J ′ = 0, the ground-state energy per site is found to be −0.4383(6), which also agrees with the best available DMRG results [33,34].
For −0.15 ≤ J ′ ≤ 0.4, the first-order derivative of the ground-state energy E 0 does not show any singularities, as shown in Fig. 2(a).However, we observe that the second-order derivative of E 0 exhibits three sharp peaks [see Fig. 2(a)].The location of these peaks indicates that there is no (at least neither first nor second-order) phase transition within the region of 0 ≤ J ′ ≤ 0.1.We emphasize that these results are stable against the system-size and bond-dimension scaling (see Supplementary Note 6) and indeed, are consistent with a previous DMRG result [54].In contrast, previous DMRG studies [53,54] suggest a critical point at J ′ c ≈ 0.08 that is identified by vanishing spin chirality.This contradiction can be resolved as follows (see Supplementary Note 2): Consider the spin chirality χ i jk (≥ 0) defined on an elementary triangle △ i jk [see Fig. 1(b)], we find that (i) χ i jk is sizable for a state deep inside the CSL phase, e.g., χ i jk ∼ 0.06 at J ′ = 0.2; (ii) χ i jk will decrease rapidly as J ′ decreases and exceed J ′ = 0.1, e.g., χ i jk ∼ 6 × 10 −4 (2 × 10 −4 ) at J ′ = 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 Supplementary Note 2); (iii) χ i jk continues to decrease to a numerical zero, χ i jk ∼ 6 × 10 −7 at J ′ = −0.04,around which CSL to VBC II phase transition occurs [see Figs.1(a) and 2].Moreover, in Fig. 3, we show the real-space chiral-chiral correlation function at J ′ = 0. Comparing J ′ = 0 and J ′ = 0.2 (see the results of J ′ = 0.2 in Supplementary Note 2 and also in Ref. [54]) data with various parameters, one can conclude that (i) the longranged 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 ′ = 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 ′ = 0, the chiralchiral correlation tends to be long-ranged on cylinders with sufficiently large circumferences, provided that the bond dimensions are adequately sized to catch this feature.
In addition to the energetics, we have also calculated the wave function fidelity F(J ′ ) = |⟨Ψ(J ′ )|Ψ(J ′ − ∆J ′ )⟩| and the corresponding susceptibility χ F ≡ dF/dJ ′ , where |Ψ(J ′ )⟩ is the ground state obtained by Randomor Boosted-DMRG initialized with |Ψ(J ′ + ∆J ′ )⟩.Here we adopt ∆J ′ = 0.01.For J ′ ≤ − 0.04, the kinks in F and singularities in χ F , as shown in Fig. 2(b), are in agreement with the location of the sharp peaks in the second-order derivative of E 0 .
For −0.03 < J ′ ≤ 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 Therefore, DMRG has no preference over a particular θ which is not under control and varies with J ′ .This would lead to a fluctuating fidelity F in Fig. 2(b).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 ′ is carried out by initializing DMRG with the post-optimized MPS for coupling J ′ + ∆J ′ .This adiabatic procedure can be done recursively from J ′ = 0.1 to J ′ = −0.03and, indeed, the fidelity F now becomes a plateau close to 1, as shown in Fig. 2(b).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 twodimensional ground-state subspace spanned by |Φ I ⟩ and | ΦI ⟩.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 ≤ J ′ ≤ 0.08, the relative difference between these two energies is smaller than 3 × 10 −5 (see Supplementary Note 1).Also, we compute the energy variance σ 2 ≡ ⟨H 2 ⟩ − ⟨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 Supplementary Note 1), 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 studies [30,33].

Topologically degenerate semion sector
Another essential signal of CSL around J ′ = 0 is that |Ψ S ⟩, the ground state in the semion sector on a 296-site YC8 cylinder, can be continuously evolved from J ′ = 0.2 to J ′ = 0 by using Boosted-DMRG.As shown in Fig. 2(c), the wave function fidelity F between two neighboring-J ′ |Ψ S ⟩'s is always close to unity for 0.02 ≤ J ′ ≤ 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 |Ψ S ⟩'s also show qualitative consistency in the region 0 ≤ J ′ ≤ 0.2, as demonstrated in Fig. 2(d).Moreover, the tiny per-site energy variances of |Ψ S ⟩'s, σ 2 /N ∼ 10 −4 , suggests that |Ψ S ⟩ 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 6).
We would like to delve into the sudden drop of the wave function fidelity starting at J ′ ≈ 0.01.It turns out that when J ′ is closer and closer to J ′ = 0 (or be precise, J ′ = J c ≈ −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 3), which increases the likelihood of these semions "colliding" with each other on finite cylinders.Thereby, it is more difficult to stabilize |Ψ S ⟩ as a "metastable" higher-energy state in DMRG, and the collapse of |Ψ S ⟩ into the identity-sector ground state |Ψ I ⟩ is more likely to occur due to the DMRG energy optimization.This agrees with our observation that the evolution of |Ψ S ⟩ can be unstable when the length of the cylinder is not long enough.In fact, the evolution breaks down at around J ′ < 0.02 when examining a shorter 152-site YC8 cylinder (see Supplementary Note 3).
For the most contentious KHAF with J ′ = 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 4).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 ′ ≲ −0.04, we calculate the NN spin-spin correlation ⟨S i • S j ⟩ and the spin structure factor which are illustrated in Fig. 4. For J ′ < −0.11, S (q) exhibits clear peaks at the K points in the extended Brillouin zone [see Fig.  1(a).Finally, in the CSL phase, S (q) is featureless [see Fig. 4(d)] and the NN correlation ⟨S i • S j ⟩ is also spatially uniform (not shown), which are consistent with a QSL state.

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 √ 3 × √ 3 magnetically ordered phase emerge in the phase diagram (in order of decreasing J ′ ).
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 insertion [53,74].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 ′ = 0. Since the semion sector is a fingerprint of a Kalmeyer-Laughlintype CSL, we argue that a gapped CSL appears to be the most promising candidate for the ground state at J ′ = 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 ′ = 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.

Supplemental Materials
This Supplemental Material includes more numerical details of the DMRG calculations.In Supplementary Note 1, we show an energetic comparison between Random-DMRG and Boosted-DMRG.In Supplementary Note 2, we perform a systematical study on the spin chirality near the point of J ′ = 0.In Supplementary Note 3, we provide more details on the characterizations of the ground state in the semion sector.In Supplementary Note 4, we compare the performance of DMRG simulations at J ′ = 0 with initial MPSs prepared from random states and several parton ansatz.In Supplementary Note 5, we show the bond-bond correlation functions to further characterize two VBC phases.In Supplementary Note 6, we discuss the system-size and bond-dimension scalings of our DMRG calculations.

Supplementary Note 1: Energetic comparisons and re-diagonalization
We introduce the relative energy difference between two close energies, E 1 and E 2 , as to quantify their closeness.As mentioned in the main text, we obtained two almost degenerate ground states, |Ψ RD ⟩ and |Ψ BD ⟩, 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 |Ψ RD ⟩ (|Ψ BD ⟩) as E RD (E BD ), the relative energy difference δ(E RD , E BD ) is extremely small, as listed in Table S1.Meanwhile, with two non-orthogonal MPSs |Ψ RD ⟩ and |Ψ BD ⟩ 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: and By solving the generalized eigenvalue problem, namely, H 2 x = EN 2 x, we can obtain two re-diagonalized eigenenergies, Ẽ1 and Ẽ2 .The relative energy difference δ( Ẽ1 , Ẽ2 ), with respect to newly obtained Ẽ1 and Ẽ2 , is also listed in Table S1.We can see that these energies are still very close to each other, implying that |Ψ RD ⟩ and |Ψ BD ⟩ are indeed the two states in the ground state manifold.Supplementary Fig. S2.Chiral-chiral correlation function as a function of distance r for (a) J ′ = 0.2 on YC8 and YC12 cylinders, (b) J = 0 on YC8 cylinders, and (c) J ′ = 0 on YC12 cylinders.Note that, for YC8 cases, we choose the reference point to x = N x /4 and measure the correlator within the middle half of the cylinder to minimize the boundary effort, and for YC12 cases, we fix the reference point to x = 2. Real number wavefunction is adopted.
We further show details of the chiral-chiral correlation function obtained in our DMRG calculations for both well-defined CSL state (at J ′ = 0.2) and J ′ = 0 state with varying system geometries and MPS bond dimensions in Supplementary Fig. S2.
At J ′ = 0.2, on wider YC12 cylinders, the correlation function is flat when the bond dimension D is larger than 12000, establishing the well convergency of our calculation and the existence of a long-range chiral order.On YC8 cylinders, while the ground state at J ′ = 0.2 has been accepted as a CSL [53,54], the chiral-chiral correlation function still present a decaying fashion, indicating the difficulty of faithful evaluation of a CSL phase in the KHAF by simply considering chiral-chiral correlation function.We emphasize that the data in Supplementary Fig. S2 are quantitatively consistent with previous works [53,54].
At J ′ = 0, for all the systems considered in Supplementary Figs.S2(b-c), the correlation continuously enhances with an increasing bond dimension, implying the existence of a long-range chiral order.Notice that these data are also quantitatively the same as the data provided in Ref. [34].Moreover, because of a weaker and weaker chirality when approaching J ′ = 0 from a positive J ′ (which implies a smaller and smaller energy gap), the required system size to observe a flat chiral-chiral correlation function is larger and larger, which makes it difficult to be accessible in DMRG calculations.until J ′ < 0.02 and there is a large jump of F around J ′ = 0.01.However, we emphasize that this jump of F does not correspond to a phase transition at J ′ = 0.01.By definition, F is highly sensitive to the local information of two wave functions and below we will show that the large jump is caused by the slight deviations of semions from the boundaries of cylinders.
We find that F at J ′ = 0 increases significantly as the bond dimension increases from D = 5000 to D = 8000.As mentioned above, the edge semions are more localized with larger D, implying that the increase of F results from the shift of semions to the boundary sides.
The absence of phase transition around J ′ = 0.02 can also be revealed by studying the density-matrix fidelity where ρ A and ρ B , being two reduced density matrices for a column of the cylinder, correspond to two different states |Ψ A ⟩ and |Ψ B ⟩.This fidelity measures how close a specific part of two wave functions is, while the boundary effects are precluded as much as possible.F d 's for the semion sector on each column are illustrated in Supplementary Fig. S7.We find that F d 's between the semion states at J ′ = 0.02 and J ′ = 0.01 as well as J ′ = 0 (orange and green closed circles, respectively) are close to 1 in the middle of the cylinders, indicating the bulk parts of those states are quite close.However, these F d 's around two boundaries are only about 0.6, due to the aforementioned slight shifts of edge semions.For comparison, we also show the F d 's between two truly orthogonal states, e.g., the semion and identity sectors.The wavefunction fidelities (i.e., F's) between those two states at the same J ′ are round 10 −10 , and corresponding F d 's are only round 0.65 in the middle of cylinders.
FIG. 1.(a) The phase diagram of the J-J ′ model consisting of a √ 3 × √ 3 ordered phase, two VBC phases [VBC I and VBC II with different patterns, see Figs. 4(e) and (f)], and a CSL phase with finite spin chirality χ.Note that the KHAF at J ′ = 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 [ π 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 L y (momentum along the y-direction) and S L z (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 ′ = 0.4.

FIG. 3 .
FIG. 2.(a) The first-(blue dots) and second-order (red triangles) derivatives of the ground-state energy E 0 versus J ′ on the 224-site YC8 cylinders.The inset shows a zoom-in of the second-order energy derivatives around the J ′ = 0 point.(b) The neighboring wave function overlap ⟨Ψ(J ′ )|Ψ(J ′ − ∆J ′ )⟩ (∆J ′ = 0.01) and the corresponding susceptibility versus J ′ .The empty circles are obtained by Random-DMRG, and the full squares are obtained by Boosted-DMRG.Two states at the same J ′ have quasi-degenerated energies (see Supplementary Note 1).(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 section-sector ground states remain qualitatively similar as J ′ decreases.
4(a)], and the NN correlation ⟨S i • S j ⟩ is spatially uniform (not shown), which can be identified as the phase[23,31, 73].In the two intermediate phases (−0.11 ≤ J ′ ≤ −0.04), S (q) are featureless [see Figs.4(b) and (c)], indicating the absence of longrange magnetic order, and the NN correlation ⟨S i • S j ⟩ exhibits clear translational symmetry breaking patterns [see Figs.4(e) and (f)], which are two different VBC phases (see more details in Supplementary Note 5) and dubbed as VBC I and VBC II in Fig.

1 FIG. 4 .
FIG. 4. Top: The spin structure factor S (q) on 224-site YC8 cylinders for (a) the √ 3 × √ 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 ⟨S i • S j ⟩ 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 [54].The cylinders are shown with only central columns.

χ
Rediagonalization Real 224 sites Real 296 sites D = 5000 Real 296 sites D = 8000 Real 296 sites D = 12000 Boosted-I 296 sites D = 5000 Boosted-S 296 sites D = 5000 Supplementary Fig. S1.The spin chirality χ calculated with different bond dimension D and system size N. (i) The blue down-triangle is obtained by the re-diagonalization method.(ii) The orange circle, green square, and red triangle are evaluated from the chiral-chiral correlation function.(iii) The purple cross and brown plus are directly calculated by Boosted-DMRG (complex-number MPS).Here Boosted-I (Boosted-S) stands for the ground state in the identity (semion) sector.