Signatures of Dirac cones in a DMRG study of the Kagome Heisenberg model

The antiferromagnetic spin-$1/2$ Heisenberg model on a kagome lattice is one of the most paradigmatic models in the context of spin liquids, yet the precise nature of its ground state is not understood. We use large scale density matrix normalization group simulations (DMRG) on infinitely long cylinders and find indications for the formation of a gapless Dirac spin liquid. First, we use adiabatic flux insertion to demonstrate that the spin gap is much smaller than estimated from previous DMRG simulation. Second, we find that the momentum dependent excitation spectrum, as extracted from the DMRG transfer matrix, exhibits Dirac cones that match those of a $\pi$-flux free fermion model (the parton mean-field ansatz of a $U(1)$ Dirac spin liquid)

Understanding the ground state of the antiferromagnetic spin-1/2 Heisenberg model on a kagome lattice (KAH) has proved to be one of the most vexed issues in quantum magnetism . The KAH is one of the simplest models with strong frustration, and is a reasonable starting point for understanding various layered magnets such as Herbertsmithite [24,25]. The possibility of a quantum spin liquid (QSL) [26] on the KAH was proposed more than two decades ago [1], and was more recently confirmed numerically through density matrix renormalization group (DMRG) simulations [8]. However, there is a multitude of possible QSLs [27,28], and despite tremendous efforts, the precise nature of the QSL for KAH remains unknown. Thus we will refer to it simply as the kagome spin liquid.
Currently the most promising candidates are the gapped Z 2 spin liquid [1,[29][30][31] ( Z 2 SL) and the gapless U (1) Dirac spin liquid (DSL) [3]. Recent DMRG studies appear to support a gapped Z 2 SL scenario [8,11,12]. However, characteristic properties of the gapped Z 2 SL, e.g. four-fold topological degeneracy on a torus and fractional statistics of spinons, have not been observed. On the other hand, there are several indications favoring a DSL over a Z 2 SL. First, extensive variational Monte Carlo studies suggest a DSL [5,[14][15][16][17]. Second, the kagome spin liquid was found to be proximate to a chiral spin liquid stabilized by the addition of longer ranged spin exchange interactions [32,33], with indications the transition between them is continuous [19]. While there is no known theory that could describe a continuous transition between a Z 2 SL and the chiral spin liquid [34,35], such a transition occurs naturally if the kagome spin liquid is a DSL. Third, a recent theoretical work [36] sug-arXiv:1611.06238v1 [cond-mat.str-el] 18 Nov 2016 gests a DSL is a natural possibility by investigating a lattice gauge theory formulation [37,38] of the easy-axis (XXZ) kagome model [19,20].
Most of the early experimental studies on candidate materials, such as Herbertsmithite, suggested a gapless spin liquid scenario based on spin susceptibility and specific heat measurements [24]. However, the effective model for those materials is thought to be more complicated than the KAH, and it was difficult to determine if these gapless signatures were intrinsic properties of the kagome spin liquid or were due to magnetic impurities. A recent NMR study of Herbertsmithite found evidence for a finite spin gap (0.03J ∼ 0.07J) by using the Knight shift to filter out impurity contributions [39].
In this paper, we revisit the kagome spin liquid problem using the DMRG [40][41][42] method, which remains one of the most unbiased and powerful numerical methods to deal with this problem.We systematically investigate the energy gap and excitation spectrum of the kagome spin liquid phase using extensions to the previously used algorithms: (i) We provide new insight into the heavily debated spin gap issue of the KAH by computing its dependence on boundary conditions, and show that the spin gap from our DMRG simulations is consistent with a gapless QSL (e.g. DSL). (ii) We obtain the momentumresolved spectrum of correlation lengths of the KAH, which is closely related to the excitation spectrum [43].
In particular, this spectrum shows signatures of Dirac cones at the locations expected for a U(1) gapless DSL [3]. The method we use here can also be directly applied to explore QSL phases in other lattice models.
The paper is organized as follows. We begin by reviewing some promising spin liquid candidates and previous DMRG studies in Sec. II. We then discuss the expected behavior of various QSLs when placed on the cylinder geometry in Sec. III. We present our numerical DMRG data in Sec. IV. First we show that the spin gap drops significantly compared to previously reported values as we twist the boundary conditions. Then we extract the momentum dependent excitation spectrum from the DMRG transfer matrix. The triplet excitations reveal a Dirac cone structure, with the Dirac point located at the M point of the Brillouin zone, as expected for a π-flux DSL. We finally conclude with a summary and discussion in Sec. V.

A. Previous DMRG studies
Yan et al. [8] performed an extensive DMRG study of the kagome Heisenberg model on cylinders with circumference sizes from L y = 2 (e.g. YC4) to L y = 6 (e.g. YC12) unit cells. Most importantly, it was found that the ground state is a symmetric spin liquid state that has a much lower energy than competing valence bond crystal states [4,7]. By performing a careful finite size scaling analysis on different geometries, an energy per site of E 0 = −0.4379(3) was estimated. Several observations were made regarding the nature of the spin liquid state: (i) The ground state is a spin singlet protected by a small spin gap to the lowest lying spin-1 state. A spin gap ∆ S=1 = 0.125(9)J, with J being the strength of the exchange interaction, was found for the L y = 6 (unit cells) cylinder (XC12-2). This gap is smaller than the gap of ∆ S=1 = 0.164J extracted from earlier exact diagonalization studies of a 36 site cluster [44]. (ii) The singlet gap, separating the ground-state from the lowest spin-0 excited state, is estimated to be much smaller, ∆ S=0 = 0.054 (9) for the L y = 6 (unit cells) cylinder (XC12-2). This differs strongly from the exact diagonalization results that estimate the singlet gap to be less than 0.01J [44]. The difference is attributed to strong finite size effects for the small clusters used for the exact diagonalization studies. (iii) The ground state is found to have a very short correlation length. The findings of Yan et al. were further corroborated by the SU(2)-invariant DMRG study of Depenbrock, et al. [11]. It was further argued based on the entanglement properties that the spin liquid is likely a Z 2 spin liquid [11,12].

B. Parton constructions for kagome spin liquids
For the discussion that follows, it will prove very useful to have a picture of the competing phases within the language of the fermionic parton construction. Here, we briefly review the parton construction for a 2D plane [27,28,45]; in Sec. III, we discuss the novel phenomena which arise when a QSL is wrapped onto a cylinder.
In the parton construction, the physical spin operator is expressed as a bilinear of a fictitious S = 1/2 fermion,Ŝ i = σ,σ f † i,σ S σ,σ f i,σ . In order to reproduce the correct properties forŜ i , one enforces the constraint 1 = σ f † i,σ f i,σ . This introduces an extensive redundancy: for example, the spin operator is left invariant under U(1) 'gauge transformations' f i,σ → e iφi f i,σ (in fact, there is actually a larger SU(2) redundancy). In the resulting effective field theory, the fermionic partons are coupled to an emergent gauge field which both enforces the constraint and implements the microscopic redundancy as a gauge invariance.
Rather than elaborating on the field theory, we use the parton construction as a variational ansatz for the ground state and its low-lying excitations. Letting |MF be an ansatz free-fermion wavefunction for the f , we can obtain a spin-1/2 wavefunction by projecting onto single occupancy, Various non-trivial spin-liquid phases result by choosing |MF to be the ground state of certain free-fermion H f , as will be discussed below.
The parton picture simultaneously suggests an ansatz for the excitations. If γ σ is a low-lying excitation of H f , with σ =↑ / ↓, we can obtain a corresponding ansatz for a low-lying excitation of the spin-system, This excitation is called a fermionic 'spinon' excitation, since it carries a S=1/2 moment. Consequently the spinon is a topological excitation, e.g., it cannot be made by acting with the local S i operators in some patch. Triplet excitations are obtained from the twospinon states. If the mean-field description is well behaved, the energy of the ansatz excitation |γ σ should have some qualitative relation to the energy of γ in H f ; for example, the location of band minima or gapless points. We now turn to the free-fermion ansatz H f thought to be relevant for the kagome spin liquid.

Gapped Z2 spin-liquid
We obtain a gapped Z 2 spin-liquid [1,29,30] by supposing that the partons form a BCS superconductor: Various choices of sign-structures for the t, ∆ (corresponding, for example, to patterns of π-flux through plaquettes) in fact lead to eight known gapped Z 2 -SLs, which differ in how the crystal symmetries act on the spinon excitations [1,[46][47][48][49]. In all eight cases the BdG particles γ σ are gapped, leading to a spin gap. There is also a second type of excitation in a Z 2 spin-liquid: a π-flux of the emergent gauge field which is coupled to the partons (the other gauge excitations are at very high energy due to the Higgs mechanism). An ansatz for this π-flux (alias 'vison') excitation is given by appropriately modifying the pairing t ij → −t ij , ∆ ij → −∆ ij along a semi-infinite line to mimic π-flux piercing the plane at the endpoint, finding the resulting free-fermion ground state, and Gutzwiller projecting. The π-flux is also a topological excitation, but it does not carry spin. The existence of the π-flux leads to important phenomena on the cylinder we will return to.

U (1)-Dirac spin-liquid
In the U(1)-Dirac spin-liquid, we choose a hoppingonly ansatz [3,5,[14][15][16][17]: The uniform nearest-neighbor ansatz t ij = 1 is a poor choice for the kagome model, since it leads to a flat valence band: with an odd number of sites in the unit cell at a density of n = 1 per site, this band would be half full. Instead, one considers the hopping pattern shown in Fig. 1(b), which has π-flux piercing the hexagons. This π-flux ansatz results in a halved magnetic Brillouin zone with two Dirac cones at momenta Q = (π/2, π/2), Q = (−π/2, −π/2). Combined with spin, there are N f = 4 gapless Dirac cones, and hence a vanishing spin gap. Since the emergent gauge field is not 'Higgsed', there are also low energy gauge fluctuations, and the resulting effective theory is essentially N f = 4 QED 3 .

Chiral spin-liquid
The chiral spin-liquid is also a hopping-only ansatz [3], but we choose the t ij to have phases such that the occupied bands have a total Chern number of C = 1 per spin species. From the point of view of the DSL, it is obtained by turning on the same-sign mass term, ψ Q,σ ψ Q,σ , on all N f = 4 Dirac cones. By choosing C = 1, the ansatz breaks both time-reversal T and reflection P , though their combination P T is preserved. Generically this leads to a non-vanishing expectation value for the chiral order parameter S i · S j × S k , where i, j, k are three nearby sites. Since the Chern band is gapped, there is a finite energy cost for spinon excitations, and the emergent gauge field is gapped by the effective Chern-Simons term. The Chern-Simons term leads to a spin-Hall coefficient of σ (s) H = 1 2 . Note that the Dirac-SL can be considered the 'parent state' for many of these spin-liquids: the Z 2 gapped states are obtained by turning on some superconducting terms, while the CSL arises if a staggered B-field spontaneously forms in addition to the π-flux. This picture gives a natural scenario for continuous phase transitions into these neighboring phases.

III. SPIN LIQUIDS ON A CYLINDER
Before providing numerical DMRG results for the kagome spin liquid, it is worthwhile to pause and consider what one should expect for a spin liquid in the DMRG simulation, which compactifies the kagome lattice into a cylinder. Since spin liquids are topological phases, they display a number subtle features on a cylinder.

A. Kagome cylinder geometries
A cylinder is defined by identifying sites that differ by some r, which defines the different compactifications. We work with the YC2n-2m cylinders, where x is identified with x + n a 1 − m a 2 and a 1 , a 2 are the kagome Bravais vectors defined in Fig. 1(a). We will have n > m, so we call L y = n the 'circumference' of the cylinder, while m amounts to a shift of the most naive identification m = 0. We note that YC2n-0, YC2n-n are actually the YC2n and XC2n geometries introduced in Ref. [8]. The M points are (labeled by (k1, k2)), M1 = (π, 0), M2 = (0, π) and M3 = (π, π). The two Dirac points of the DSL are at ±Q. The DSL has a halved magnetic BZ (dashed), since it is a π-flux ansatz.
Since one direction of the system is compactified, the momentum along that direction is discretized, specifically Note that the other direction is infinite, so k 1 and k 2 can take continuous values so long as they satisfy the above relation. The different geometries (different n and m) thus provide different cuts through the Brillouin zone, as shown in Fig. 1(c)-(d).
A second parameter one can tune is the boundary condition θ of the spins around the circumference, the same way one would measure spin-stiffness. We obtain a twist boundary condition by modifying Heisenberg coupling according to S + x,y S − x ,y+1 → e iθ/Ly S + x,y S − x ,y+1 . By analogy to flux threading, we call θ the spin-flux in the cylinder , as shown in Fig. 2 (a).

B. Gapped SL on a cylinder
Gapped spin liquids have fractionalized quasiparticles that have non-trivial statistics (e.g., the fermionic spinon). A direct consequence of such fractionalization is the topological degeneracy on a torus; there is one ground state per quasiparticle type [29,50]. Precisely the same degeneracy arises for a long (infinite) cylinder -one can picture the ends of the cylinder as identified into a torus at infinity. The energy splitting between these 'topological sectors' is exponentially small in the circumference.
In the Z 2 -SL, the topological degeneracy can be understood in terms of the boundary conditions of the fermionic partons. While the boundary conditions of the spins are set by θ, the fermionic spinons are coupled to an emergent gauge field, and the flux 'φ' of this gauge field through the cylinder effectively changes their boundary conditions. In a superconductor, π-flux is invisible to the condensate, so periodic (PBC) and anti-periodic (APBC) will give nearly degenerate energies -the splitting should decay exponentially in the cylinder circumference vs. coherence length.
There is an interesting interplay between the topological degeneracy and spin-flux θ. Spin rotations S ± i → e ±iΩ S ± i are implemented on the partons as f i,σ → e iΩσ/2 f i,σ , where σ = ±1 corresponds to spin-up and spin-down components (the eigenvalue of 2S z ). As a consequence, the spin flux θ will change the boundary conditions of the up/down partons by ±θ/2. A superconductor has no spin-stiffness, so the change in the energy will again be exponentially small in the circumference of the cylinder. However, at θ = 2π, H f (in Eq.(3)) does not return to itself: because they carry half-integer spin, the boundary conditions for the f σ have been changed by ±2π/2 ≡ π. Equivalently, threading 2π-spin flux adiabatically exchanges between the PBC and APBC topological sectors, as shown in Fig. 2(b). Note that the location of the crossing is pinned to θ = π by time-reversal, where it acts by exchanging the two topological sectors.
In summary, there is a striking signature of a gapped Z 2 spin liquid: as spin-flux is inserted, there should be an exponentially small change (L x O[e −L cyl /ξ ]) in the ground state energy, and at θ = 2π the ground state should not return to itself, but instead should thread the emergent π-flux that corresponds to the other topological degenerate groundstate [51,52].

C. Dirac spin liquid on a cylinder
Like the Z 2 -SL, the DSL has an emergent gauge field which can dynamically change the boundary conditions of spinons (including either PBC or APBC). However, because of the semi-metallic parton band structure the response is quite different: any twist of the parton boundary conditions will shift the energy by v F L x /L 2 y , where v F is the spinon velocity, due to the energy of the filled parton bands below the Dirac point. In particular, a boundary condition in which the allowed k-modes avoid the Dirac point will be the lowest in energy, since this effectively opens up a gap of order v F /L y . Thus when a Dirac spin liquid is placed on a narrow cylinder, the gauge field will generically adjust to open up a gap. Therefore, a numerical observation of a non-vanishing gap on a single long cylinder generally does not rule out the possibility of a gapless DSL.
Here, flux threading can be used to find fingerprints of a DSL on a cylinder. As before, adiabatically threading spin-flux θ through the cylinder twists the up/down parton boundary conditions by ±θ/2 as shown in Fig. 2(a). This twist is in addition to the (shared) emergent flux φ. For symmetry reasons, the internal gauge flux will either be φ = 0, π. As we thread spin flux, the energy spectrum will again generically look like Fig. 2(b), but in contrast to a gapped SL the splitting is only algebraically small. Beyond the crossing, it will become difficult to adiabatically track the state due to the large splitting and small gap (since the allowed k-modes becomes close to the Dirac point); at some point π-emergent flux will enter to partially cancel the spin-flux, at which point there will be a discontinuous jump in the ground state energy of the DMRG simulation [53].
The expected behavior of the gap depends on the geometry, which we consider in more detail at the meanfield level in order to interpret our numerical results. At spin-flux θ, the momenta of the partons f ↑/↓ on a YC2n-2m cylinder are restricted to where φ = 0, π is the emergent gauge flux. The ± correspond to spinons f ↑ , f ↓ respectively. There are two classes of cylinders which behave very differently with θ: • Type I cylinder: YC2n-4k. For YC2n-4k cylinders, when θ = 0 the flux φ = π avoids both Dirac points, while both are present for φ = 0. In this case, we would (naively) expect the gap to decrease like (2π − θ)/L x until adiabaticity is lost after θ > π and the emergent π-flux tunnels in. Hence we can't force the system to go gapless.
Here, for φ = 0 the ↑ Q and ↓ Q components are gapless, while for φ = π the ↓ Q and ↑ Q are. All other values of spin-flux have a gap which should decrease as |π − θ|/L x . Thus at θ = π two of the four Dirac fermions are present regardless of the emergent flux, and hence we can force the system to go gapless.

IV. NUMERICAL RESULTS
In the following, we will discuss the numerical results for the KAH obtained using DMRG simulations on infinite cylinders (iDMRG) [42]. Besides the nearestneighbor Heisenberg interactions, we also include a smallsecond neighbor interaction for some simulations, We will study the behavior of the spin gap and transfer matrix spectrum (an analog of excitation spectrum) as we adiabatically twist the boundary conditions (spinflux θ). We implement the adiabatic twist by using the previous DMRG wavefunction as the initial step for the next θ value [53]. The ground state energy during the insertion is provided in the Appendix. Note that while DMRG generally finds the absolute ground state, when passing adiabatically through a level crossing there may be a small regime in which the DMRG follows the higher, metastable level until 'tunneling' into the lower state, a phenomena we encounter below. We find that the behavior of all the geometries can be grouped into the two types discussed in Sec. III C. For the Type I cylinders (YC8-0 and YC8-4), the adiabaticity fails at θ ≈ 4π/3, and the simulation collapses to the other topological sector with lower energy; [53] note for π < θ 4π/3 we are tracking the higher, metastable state, not the absolute ground state. For Type II (YC6-2, YC8-2, YC10-2, YC12-2, YC8-6), the adiabaticity fails around θ ≈ π, at which point we find an instability of the kagome spin liquid towards an ordered state. For a DSL, this instability corresponds to a spontaneous generation of mass iΨσ 3 µ 3 Ψ (see the Appendix Sec. C for more details).

A. Spin gap
The spin (triplet) gap ∆ S=1 is obtained by creating a S z = 1 excitation in the bulk of the cylinder and then calculating the energy difference to the S z = 0 sector [32,54] (see Appendix Sec. A for details). We stress here that the spin gap is different from the (extensive) energy splitting between the different topological sectors. Fig. 3(a) provides raw data for the spin gap as function of the twist angle θ for DMRG bond dimension (a) Upper bound on the spin gap ∆S=1 under the twist boundary condition θ for J2 = 0; the estimate is obtained using DMRG bond dimension m = 4000 and using the lowest-energy topological sector. Data ends at the failure of adiabaticity. The qualitative behavior depends on the type of cylinder (I or II), which in a DSL would be gapless at θ = 2π, π respectively. (b) Dependence of the spin gap ∆S=1 (solid line) and S z = 1 correlation length ξ (dashed line) on the spin flux θ and DMRG bond dimension m. Data is taken with J2 = 0.05 on the Type II YC8-2 cylinder. Generally the estimated gap decreases with the bond dimension m; the larger the system size is, the more likely we will overestimate on the spin gap. m = 4000. Generally the spin gap decreases with increasing bond dimension, as seen in Fig. 3b, and thus the data shown provides only an upper bound for the spin gap. When θ = 0, the spin gap is rather large, for example YC8-0 and YC8-4 have a spin gap of ∆ S=1 ≈ 0.15J, completely consistent with previous DMRG simulations [8,11]. Strikingly, the spin gap shows a significant decrease when θ is increased, and right before the failure of adiabaticity, the spin gap drops to a very small value ∆ S=1 ≈ 0.02J − 0.04J. The approximately linear decrease of the gap for larger twist angles θ is suggestive of a DSL gapless spin liquid. As discussed in the previ-ous section, the spin gap for a gapped spin liquid should have an exponentially small (O(e −Ly/ξ ) dependence on the twist angle θ, though of course ξ can be large.
We emphasize that, based on our data alone, we cannot conclude whether the gap will vanish in the thermodynamic limit. For large enough L y , the spin gap of a DSL at θ = 0 should decrease with the circumference size L y as ∆ S=1 ∼ v F /L y . This behavior has not been observed in previous studies [8,11], nor in our current simulations. On the one hand, this could be an artifact of the finite bond dimension m of DMRG simulations; finite m tends to overestimate the spin gap ( Fig. 3(b)), necessitating a careful extrapolation in 1/m α , where α is an unknown power. Since m must increase exponentially with circumference to achieve the same accuracy, it becomes very challenging to accurately extract the spin gap for large circumferences. Alternatively, for small circumferences it is unclear whether the gap will follow the naive mean field expectation ∆ S=1 ∼ v F /L y ; a possible reason is discussed in Sec. IV B 1.

B. Transfer matrix spectrum
We examine the transfer matrix spectrum on the infinite cylinder (i.e., the correlation-length spectrum) and compare the KAH with a free fermion π-flux model. In the ansatz wavefunction of the DMRG (matrix product states), the correlation functions of charge-q operators (e.g., the spin-spin correlations C S=1 (r) = S + 0 S − r con ) can be expanded as a sum of exponentials, Here r is the distance along the long direction of the cylinder, q is a quantum number, m q is a bond dimension, and λ q,j are eigenvalues of the DMRG 'transfer matrix' [55]. Using the quantum numbers q, we can distinguish, for example, a triplet excitation from a singlet excitation (or more precisely, a S z = 1 excitation from a S z = 0 excitation). The eigenvalues λ q,j = e ikq,j −ξ −1 q,j have a real part, corresponding to a correlation length ξ q,j , and an imaginary part, corresponding to a momentum k q,j . The largest, ξ q = max(ξ q,j ), bounds the correlation length of all charge-q operators.
A recent work by Zauner et al. [43] pointed out a relation between the energy spectrum of the physical excitations and the spectrum of the transfer matrix. A more familiar statement is that the largest correlation lengths ξ q set an upper bound for the lowest excitation gaps ∆ q (up to a factor), and for a Lorentz invariant systems it holds that ∆ q ∝ 1/ξ q . The corresponding k q gives the momentum of the excitation along the length of the cylinder. These relations actually hold nicely for the KAH at different twist angles, as demonstrated in Fig. 3b. Figures 4a-b show the S z = 1 transfer matrix spectrum of the KAH and the π-flux free fermion model Eq. (3) as (d) function of the twist angle θ. We consider the Type II YC8-2 geometry, and included a small J 2 to stabilize the adiabatic twist up to θ = π (the Appendix shows results for other geometries). The three different colors label three 'bands:' we observe that the momenta k S=1,j cluster into three distinct groups, and we plot the largest several ξ S=1,j from each momenta group. The momentum can be resolved into its lattice components k 1 , k 2 , providing an alternative way to plot the data shown in Figures 4c-f.

Interpretation as DSL
The KAH and free fermion spectra are remarkably similar. The excitation spectrum can be understood based on the free fermion π-flux model. S z = 1 excitations arise from particle-hole excitations near the Dirac points; a momentum p − q spin flip takes the form S + (p − q) = f † ↑ (p)f ↓ (q). The π-flux state has two Dirac points at Q = (π/2, π/2) and −Q = (−π/2, −π/2). We group the particle-hole excitations into intra-valley forward (blue); intra-valley backward (yellow); and inter-valley forward (red), as illustrated in the cartoon 4(g) (inter-valley backward scattering is higher in energy). The Dirac points are avoided on the θ = 0 YC8-2 cylinder, as shown in Fig.  1(d), but as θ increases, the allowed momenta shift and eventually pass through the Dirac point; the f ↑ and f ↓ feel opposite flux, hence move oppositely. As can be seen in the cartoon, this shift affects the three modes in a qualitatively different fashion. The dispersion of the red mode follows a Dirac behavior and becomes gapless at the twist angle θ = ±π. In terms of momenta, the gapless point occurs at (2k 1 , k 2 ) = (0, π) as expected from the displacement between Q and Q . The yellow mode has a constant energy under the twist angle θ. The blue mode has similar response as the red mode, but remains gapped when the system hits the Dirac points (θ = π).
The spectrum of the KAH and the π-flux free fermion model show surprisingly good agreement: (i) the red mode has a linear sharp Dirac cone structure; (ii) the yellow mode is almost flat; (iii) the modes occur with the predicted momenta. The qualitative difference between two models is that the yellow and blue modes in the KAH are lower compared with the free fermion model. It may that even though the DSL theory should have an emergent SU(4) symmetry in 2D, in the quasi-1D geometry intra-valley interactions are stronger.
The existence of the renormalized flat yellow band also explains the 'kink' in the θ-dependence of the triplet gap ∆ S=1 : for small θ it drops below the linear red band, which then cross. This implies that gaps obtained in previous DMRG studies, which all worked at θ = 0, were probing the yellow intra-valley excitation. Since the yellow band is subject to strong interaction effects, this may relate to the non-observation of v F /L y gap scaling on accessible cylinders.
We want to remark that within our DMRG simulations the correlation spectrum of the KAH still has a finite "gap" even at the Dirac point. This is a necessary consequence of DMRG, since the finite bond dimension m induces a finite correlation length. We find the ξ increase with m, as expected for a DMRG simulation of a critical system. In fact, a similar behavior is found also in the free fermion model. The correlation length estimated from DMRG (with m = 250 in Fig. 4) is finite even for the gapless free fermion model with an infinite correlation length. For a very large bond dimension (m = 3000, see Appendix Fig. 8), the correlation length becomes very large (ξ ∼ 1000 sites), supporting that the finite correlation length at the Dirac point (in Fig. 4) is purely an artifact of small bond dimension.

C. S z = 0 spectrum and scalar chirality
We now turn to the S z = 0 transfer matrix spectrum shown in Fig. 5(a), which includes singlet excitations for which the interactions have a more drastic effect. Note that because our numerics do not explicitly preserve SO(3) symmetry, the S z = 0 spectrum contains both SO(3) singlets and elements of SO(3) multiplets. The correlation length in the S z = 0 sector shows a critical feature: at the the twist angle θ = π, the S z = 0 correlation length ξ ∼ 100 unit cells, which (holding fixed the bond dimension m = 6000) is much larger than the correlation length ξ ∼ 5 in the S z = 1 sector. This effect is beyond a mean-field analysis, since for a free fermion, there will be no difference between S z = 0 and S z = 1 sectors.
The large correlation length in the S z = 0 sector arises from long-range correlations of the scalar chirality χ = S i · ( S j × S k ). To confirm this, in Fig. 5(b), we show a fit for the correlation function of the scalar chirality, χ(0)χ(r) ∼ e −r/ξ . Indeed, the correlation length ξ is almost identical to the largest correlation length obtained from the transfer matrix in the S z = 0 sector. The large correlation length of the scalar chirality is somewhat surprising given our knowledge of N f = 4 QED 3 (U (1) DSL). The scalar chirality corresponds to the SU(4) singletψψ, which was suggested to have higher scaling dimension than the SU(4) adjoint fermion bilinears [56,57]. Naively, an operator with a lower scaling dimension in a critical theory would give a larger correlation length in DMRG calculation with a finite bond dimension, in an apparent contradiction with our result. However, our result does not necessarily indicate the scalar chirality has the lowest scaling dimension in the critical theory, since the magnitude of the correlation length is not directly equivalent to the scaling dimension of an operator. Another possible explanation is that the cylinder geometry drastically changes the scaling analysis, since (if it was stable) QED 3 dimensionally reduces to N f − 1 coupled Tomonaga-Luttinger liquids [58], for which the scaling analysis may differ substantially from the 2D limit. In addition, we note that the SU(4)-invariant massψψ generates the CSL ( it gives the bands a net Chern-number), and we know that the KAH is proximate to the CSL [19,32,33,59]. It was also suggested that the scalar chirality contains a monopole operator, which might have a lower scaling dimension than the fermion bilinears [6].

V. CONCLUSION AND DISCUSSION
We used large scale DMRG simulations to study the quantum spin liquid phase in the S = 1/2 kagome antiferromagnetic Heisenberg model. DMRG studies of the KAH are most natural on a cylinder. We point out that, even if a gapless QSL such as the Dirac spin liquid is realized in the two-dimensional bulk limit, the system generally acquires a non-vanishing gap on the cylinder. The predicted size of this gap depends sensitively on the geometry and boundary conditions of the cylinder, complicating the finite size scaling analysis. Thus the observation of a gap in the previous DMRG studies of the KAH, which is also reproduced in our study, might not rule out a gapless QSL until the gap can be accurately measured for a sequence of 'equivalent' geometries (e.g., YC8, YC12, YC16 · · · ), which is extremely challenging due to the exponential blow-up in DMRG bond dimension.
To better identify the nature of the QSL in the KAH using DMRG, we insert spin-flux through the cylinder, changing the boundary condition around the circumference. Using adiabatic flux insertion we find that the spin gap on the cylinder geometry is much smaller than estimated from previous DMRG simulations. Second, we found that the momentum dependent excitation spectrum as estimated from the DMRG transfer matrix exhibits Dirac cones that agree well with the ones found for a π-flux free fermion model (the parton mean-field ansatz of a U (1) Dirac spin liquid). These findings suggest that the ground state of the KAH is a gapless DSL, instead of a gapped QSL such as the Z 2 topological phase. This is more in line with several recent numerical [5,10,[14][15][16][17] and analytical [36] results obtained by methods other than DMRG.
However, we also have to make some cautious remarks. While we found indications of gapless features in the system sizes that we can access, we cannot draw a definite conclusion for the thermodynamic limit. Diverging correlation lengths require extremely large bond dimensions for proper convergence of the gap, making it impossible to rule out the existence of a small but finite gap with a large correlation length. In addition, fixing θ = 0 the spin gap should eventually decrease as v F /L y , which has not yet been observed. Within a DSL scenario, this could be a finite-size effect due to observed strong renormalization of the intra-valley backward scattering, or it could be a numerical artifact of finite-bond dimension. Finally, the present analysis does not reveal the nature of the gauge field, which could be either U (1) or Z 2 . Our result therefore is certainly not the final answer to the long-standing question on the KAH. Nevertheless, our results strongly suggest that the ground state of the KAH is a DSL, which has many theoretical and experimental implications. Numerically, we wrap a 2D lattice on a cylinder with one direction compactified with a small number of sites, and the other direction infinite. We use the MPS to cover the cylinder in the fashion of a snake chain, as graphically shown in Fig. 6(a). Along the compactified direction (a 1 ), the snake covering does not maintain the translational symmetry explicitly, thus one needs distinct MPS for each site. On the other hand, the MPS is translation invariant and repeating along the direction a 2 . Then one can use the MPS (of the smallest repeating unit cell) to define the transfer matrix (TM), as shown in Fig. 6(b). With the TM, one can further calculate its eigenvalues λ q,j = e ikq,j −ξ −1 q,j , which has a real part, corresponding to (e) (c) a correlation length ξ q,j , and an imaginary part, corresponding to a momentum k q,j . q is the quantum numbers, from which we can distinguish for example triplet excitation from singlet excitation.
Due to the snake covering, the MPS does not have translational invariance along the compactified direction (a 1 ). However, the Hamiltonian still has translational invariance along a 1 , hence the momentum k 1 along that direction can still be extracted. The way to calculate k 1 is similar as calculating a global quantum number q, for which one needs to obtain the momentum k 1 of each Schmidt basis. Technically, we implement the translational operation T 1 along the direction a 1 , and then obtains a mixed TM-T T1 , as shown in Fig. 6(c). The dominant eigenvector of T T1 will be V α,β = δ α,β e ikα , and k α gives the momentum k 1 of each Schmidt basis.
We want to remark that for certain special geometries, for instance the YC2n-2, the momentum along the compactified direction cannot be extracted using the method discussed above. The reason is because for those geometries, the momentum k 1 and k 2 are intertwined together. A consequence is that one cannot define the mixed TM for the translation T 1 . For example we can consider the -2 geometry for the square lattice Fig. 6(d). One can see that under the T 1 translation, the sites on one column doesn't go back to itself, instead some of the sites will go to the neighboring column. This is different from the normal geometry (see the left panel of Fig. 6(a)), for which T 1 translation maps the sites on one column back to itself, hence the mixed TM-T T1 can be well defined. A benefit of the -2 geometry, however, is that the snake-fashion MPS is actually translational invariant under 2 sites, no matter how large the circumference of the cylinder is. Then one can actually use a MPS with 2site structure (for the kagome lattice, it is 6 sites) to do the iDMRG simulation. With the momentum k from the TM's eigenvalue, one can then obtain the momentum Here L y is the width of the cylinder, θ is the twist boundary condition. Therefore, for the -2 geometry, one can still get k 1 and k 2 , but k 1 has a π ambiguity. Before closing this section, it is worth to make a few remarks on fermionic systems. Usually to simulate a fermionic system, we first do a Jordan-Wigner transformation to obtain the corresponding bosonic model, and then simulate the bosonic model directly. Therefore, to calculate the spectrum of the fermionic (e.g. single particle) excitations, one should consider the mixed TM-T f , with a Jordan-Wigner string inserted, Fig. 6(e).

Algorithm of calculating the spin gap
We use an algorithm that combines iDMRG and finite DMRG to calculate the spin gap. Firstly, we use iDMRG to obtain a converged wave-function of an infinite cylinder, then we cut the infinite cylinder into two halves, insert a 3 × L y × L y sites into the system. The left (L) and right (R) semi-infinite cylinder can be considered as environment (boundary conditions), and we calculate the energy of ground state E 0 (S z = 0), spin-1 sector E 0 (S z = 1) within the inserted 3 × L y × L y cylinder. Finally, we obtain the spin gap ∆ S=1 = E 0 (S z = 1) − E 0 (S z = 0). This algorithm is similar as the one used in finite DMRG, where one obtains the spin (triplet) and singlet gap by sweeping in the middle of a finite cylinder to minimize the boundary effect. The only difference is that the boundary environment we use comes from the iDMRG simulation, while finite DMRG uses environment from finite cylinder simulation.

The kagome Heisenberg model
a. Spin-stiffness Fig. 9 shows the response of the ground state energy and entanglement entropy under the twist boundary conditions before the failure of adiabaticity. Practically, the adiabaticity can be checked by looking at the wave-function overlap(≈ 0.99) between each two adjacent twist angles. During the (adiabatic) twist process, the system remains in a spin liquid phase that preserves all the lattice symmetries. Both the energy and entanglement entropy increase under the twist; the increase in entanglement entropy is very significant, and may be underestimated by finite m. This behavior is consistent with the DSL. It is worth noting that, for YC8-0, the twist trajectory is not symmetric about θ = π, which is a clear signature of fractionalization; the state at θ = π is not time-reversal-invariant as well (also true for YC8-2).

b. Transfer matrix
We provide more data for the transfer matrix spectrum of the KSL. We will show that by changing the system sizes and parameter point (J 2 interaction), the feature of a Dirac cone structure always exists. For example, the S z = 1 excitation shows a clear Dirac cone structure. As before, the S z = 0 excitation is always much lower than the S z = 1 excitation.
As we described in the main text, the behavior of the system can be grouped into two types, YC2n-(4k+2) cylinder and YC2n-4k cylinder. In the following, we will look at the two different classes separately. We plot the S z = 1 transfer matrix spectrum of J 2 = 0.05 Fig. 10, J 2 = 0.1 Fig. 11 and J 2 = 0 Fig.  12 for the YC6-2, YC8-2, YC10-2, YC12-2, YC8-6 cylinders. Clearly, the Dirac cone is independent of the J 2 interaction or system sizes. For J 2 = 0 in Fig. 12, we find that it is more difficult to main the adiabaticity of the twist around θ = π. Specifically, for the YC6-2 and YC8-2 cylinder, we can only adiabatically twist the system to θ = 5π/6; for θ = π, we always end up with an ordered state with a sudden jump. As we argued in the following section, once a Dirac spin liquid is put on a small cylinder, it might have the instability by spontaneously generating a mass gap. Such finite size effect will be gone in the pure 2+1D limit. This is consistent with our observation that, for a larger system size (i.e. YC10-2, YC12-2), the adiabatic twist can be maintained even for J 2 = 0.
Comparing the excitation spectrum of different system sizes, it seems that the larger the system is, the higher the Dirac mode (red circles) is. However, this is probably an artifact of the finite entanglement effect (truncation error) in our DMRG simulation. Since the Dirac mode is gapless, it suffers stronger finite entanglement effect than a gapped mode. That is to say, the larger the trun- cation error is, the more we will underestimate the Dirac mode. This artifact can be seen in Fig. 13, where clearly the Dirac mode becomes lower as the bond dimension increases. On the other hand, to achieve same accuracy (truncation error), the required bond dimension increases exponentially with the circumference of the cylinder. Table I shows the truncation error for different bond dimensions, system sizes and twist angles. For the large system size (YC10-2, YC12-2), m = 6000 roughly gives comparable accuracy as m = 1000 for YC8-2. Therefore, more care should be taken if one wants to compare the results between different system sizes.
Besides the S z = 1 transfer matrix spectrum, it is also interesting to look at the S z = 0 transfer matrix spectrum, Fig. 14. As we discussed in the main text, those S z = 0 excitations are more critical than the S z = 1 excitations. For a large system size, the singlet excitation  I. The truncation error for J2 = 0.05. The truncation error increases as the system sizes. It is also worth noting that the θ = π has much larger truncation error than the θ = 0 due to the fact that θ = π is likely to be gapless. is not gapless at the Dirac point θ = π. We think this is again the artifact of finite bond dimension, as one can see that the singlet excitation keeps on going down as the bond dimension m increases, as shown in Fig. 14(b).

d. YC2n-4k cylinder
At last, let's look at the YC2n-4k cylinder. As we discussed in the main text, this class of cylinder behaves very different for YC2n-(4k+2) that is discussed above. For the π-flux Dirac spin liquid, the YC2n-(4k+2) cylinder will hit the gapless Dirac point at θ = π, while YC2n-4k cylinder will hit the gapless Dirac point at θ = 2π. Our simulation on the YC2n-4k cylinder is also consistent with this scenario, namely the adiabaticity of the twist can be maintained after θ = π until θ ≈ 4π/3, after which the system collapses to the other topological sector.
When the π-flux DSL hits the Dirac points at θ = 2π on the YC2n-4k cylinder, the four Dirac fermions will be simultaneously gapless. Those four gapless Dirac fermions could then form different gapless fermion bilinears giving rise to gapless triplet (singlet excitations). This is again sharply distinct from YC2n-(4k+2) cylinder, where only two Dirac fermions are gapless when the system hits the Dirac points (at θ = π). Therefore for the π-flux DSL, we expect that gapless triplet excitation at all the three M points will be gapless on the YC2n-4k cylinder. However, we are working at a small cylinder, such that some lattice symmetry (e.g. C 3 ) is explicitly broken. Therefore, it is possible that the gapless triplet excitation at certain M point will be pushed to a higher energy level. Fig. 15 shows the S z = 1 transfer matrix spectrum of the cylinder YC8-0 with J 2 = 0. We find the lowest modes behaves like the Dirac modes. Similar as the YC2n-(4k+2) cylinder, the lowest modes show a linear dependence with the twist angle θ. These "Dirac modes" are actually two-fold degenerate, and they have the same k 1 but distinct k 2 . Since our simulation cannot adiabatically twist to the Dirac points, we cannot unambiguously determine the momentum of the Dirac points. But there are several indications that, the two Dirac modes correspond to the M 1 = (π, 0) and M 3 = (π, π) points (labeled by (k 1 , k 2 )). First, the YC8-0 cylinder has a reflection symmetry (with the reflection axis perpendicular to a 1 ), under which M 1 transforms to M 3 . Second, by doing a simple linear extrapolation for momentum k 1 (Fig. 15(b)) and k 2 (Fig. 15(c)), the momentum is consistent with M 1 and M 3 . We note that the extrapolation for k 2 gives k 2 ≈ ±0.1π and k 2 ≈ π ± 0.1π, which has considerable discrepancy from k 2 = 0 and k 2 = π. This discrepancy might come from the finite size effect, for example there is scattering (momentum transfer) between two modes. For the YC8-2 cylinder, once we tune the flux to θ = π in order to exactly hit the Dirac points, the kagome spin liquid is unstable to an ordered state. Such ordered state breaks spin flip symmetry P x = (2S x ) and the lattice symmetry (e.g. T a1 and C 6 ), its order pattern is shown in Fig. 16. Interestingly, such ordered state ac- tually preserves certain symmetries, which are i) the reflection symmetry R y , ii) translational symmetry T a2 , iii) the combination of spin flip and translational symmetry P x T a1 , iv) XY spin rotation symmetry.

Spontaneous mass generation of Dirac fermions in 1+1 dimension
As we discussed in Sec. IV, when the kagome spin liquid is tuned to exactly hit the Dirac points on a small cylinder (YC2n-(4k+2) geometry), we numerically find an instability toward an ordered state by spontaneously generating a mass term iΨσ 3 µ 3 Ψ. This immediately raises the question that, whether this results imply that the U(1) DSL (N f = 4 QED3) is unstable to the spontaneously chiral symmetry breaking (CSB). The CSB of QED3 is still an open issue, and it is unclear whether N f = 4 QED3 will eventually flow to an interacting conformal fixed point or not [60]. Our numerical simulation, on the other hand, was carried on a quasi-1D cylindrical geometry, for which the issue of spontaneous mass generation is different from the 2+1D limit. The differences are two-fold: i) The U(1) gauge field is more gentle in quasi-1D, it simply reduces the N f flavors of 1D Dirac fermions to N f − 1 coupled Tomonaga-Luttinger liquids (TLL) [58]. ii) The effect of four-fermion interactions on Dirac fermions are more drastic in 1D than that in 2D, namely in 2D all four-fermion interactions are irrelevant for (free) Dirac fermions; while in 1D, four-fermion interactions might be relevant or marginally relevant.
For the U(1) DSL on the YC2n-(4k+2) cylinder with θ = π-flux, two Dirac fermions are gapless. Then the dynamical U (1) gauge field reduces the system to a TLL with central charge c = 1, described by K is the Luttinger parameter. Generically, there are perturbations to that can potentially gap out the system. Assuming that the U(1) symmetry is preserved, the most relevant perturbation is λ cos 2φ, which is irrelevant when K > 1/2 [61]. For K < 1/2, the perturbation is relevant, and will gap out the Luttinger liquid, which corresponds to the mass generation of the Dirac fermions. For K = 1/2, the perturbation is marginally relevant or marginally irrelevant, depending on the sign of λ. Microscopically, it is difficult to extract the Luttinger parameter K. However, for the c = 1 Luttinger liquid with SU (2) symmetry, the K is fixed to K = 1/2. Our system is very close to this limit (the twist boundary condition slightly breaks SU (2) symmetry), hence we expect K ≈ 1/2. Consequently, it is reasonable that in our numerical simulation on a small YC2n-(4k+2) (e.g. YC8-2) cylinder, the DSL is unstable to an ordered state with a mass term spontaneously generated.