Dirac Spin Liquid on the Spin-1/2 Triangular Heisenberg Antiferromagnet

TL;DR

This study uses DMRG simulations with flux insertion to provide strong evidence that the spin-1/2 triangular Heisenberg antiferromagnet hosts a gapless U(1) Dirac spin liquid with characteristic excitations and universal entanglement properties.

Contribution

It demonstrates the presence of a gapless U(1) Dirac spin liquid in the model using flux insertion and excitation spectrum analysis, advancing understanding of quantum spin liquids.

Findings

Evidence for gapless U(1) Dirac spin liquid behavior.

Observation of Dirac cone structures in excitation spectrum.

Universal entanglement entropy response under flux insertion.

Abstract

We study the spin liquid candidate of the spin-$1/2$ $J_1$-$J_2$ Heisenberg antiferromagnet on the triangular lattice by means of density matrix renormalization group (DMRG) simulations. By applying an external Aharonov-Bohm flux insertion in an infinitely long cylinder, we find unambiguous evidence for gapless $U(1)$ Dirac spin liquid behavior. The flux insertion overcomes the finite size restriction for energy gaps and clearly shows gapless behavior at the expected wave-vectors. Using the DMRG transfer matrix, the low-lying excitation spectrum can be extracted, which shows characteristic Dirac cone structures of both spinon-bilinear and monopole excitations. Finally, we confirm that the entanglement entropy follows the predicted universal response under the flux insertion.

Figures

Tables

Table 1. Table 1: Discrepancy of bond strengths (subtracted by the average value ∼ − 0.18 similar-to absent 0.18 \sim-0.18 ) in the lowest-energy odd-sector state for various YC L y subscript 𝐿 𝑦 L_{y} - 0 0 cylinders by setting J 2 / J 1 = 0.12 subscript 𝐽 2 subscript 𝐽 1 0.12 J_{2}/J_{1}=0.12 and θ = 0 𝜃 0 \theta=0 . They are shown in three inequal directions 𝒂 2 subscript 𝒂 2 \bm{a}_{2} , 𝒂 1 subscript 𝒂 1 \bm{a}_{1} and 𝒂 1 − 𝒂 2 subscript 𝒂 1 subscript 𝒂 2 \bm{a}_{1}-\bm{a}_{2} . The data is accurate to 4 4 4 decimal places.

Geometry	$m$	${𝒂}_2$	${𝒂}_1$	${𝒂}_1-{𝒂}_2$
YC$6$-$0$	$1024$	$0.0002$	$-0.0004$	$0.0002$
	$2048$	$-0.0004$	$0.0007$	$-0.0004$
	$4096$	$-0.0005$	$0.0011$	$-0.0005$
	$6144$	$-0.0006$	$0.0012$	$-0.0006$
YC$8$-$0$	$1024$	$0.0286$	$-0.0576$	$0.0290$
	$2048$	$0.0251$	$-0.0503$	$0.0252$
	$4096$	$0.0227$	$-0.0455$	$0.0227$
	$6144$	$0.0218$	$-0.0437$	$0.0219$
	$8192$	$0.0214$	$-0.0428$	$0.0214$
	$12288$	$0.0210$	$-0.0419$	$0.0210$
YC$10$-$0$	$1024$	$0.0517$	$-0.1072$	$0.0555$
	$2048$	$0.0459$	$-0.0927$	$0.0468$
	$4096$	$0.0404$	$-0.0818$	$0.0413$
	$6144$	$0.0381$	$-0.0765$	$0.0384$
	$8192$	$0.0365$	$-0.0732$	$0.0368$
YC$12$-$0$	$1024$	$0.0682$	$-0.1481$	$0.0798$
	$2048$	$0.0623$	$-0.1319$	$0.0696$
	$4096$	$0.0573$	$-0.1171$	$0.0598$
	$6144$	$0.0534$	$-0.1093$	$0.0559$
	$8192$	$0.0509$	$-0.1042$	$0.0533$

Table 2. Table 2: Truncation error for various YC 2 ​ n 2 𝑛 2n - 1 1 1 cylinder by setting J 2 / J 1 = 0.12 subscript 𝐽 2 subscript 𝐽 1 0.12 J_{2}/J_{1}=0.12 . θ = π 𝜃 𝜋 \theta=\pi has much bigger truncation error than θ = 0 𝜃 0 \theta=0 .

	$m=1024$	$m=6144$	$m=12288$
YC$8$-$1$   $\left(\theta =0\right)$	$5.3\times {10}^{-5}$	$5.3\times {10}^{-6}$	$1.7\times {10}^{-6}$
YC$8$-$1$   $\left(\theta =\pi \right)$	$6.5\times {10}^{-5}$	$1.1\times {10}^{-5}$	$4.5\times {10}^{-6}$
YC$10$-$1$ $\left(\theta =0\right)$	$1.3\times {10}^{-4}$	$2.5\times {10}^{-5}$	$1.3\times {10}^{-5}$
YC$10$-$1$ $\left(\theta =\pi \right)$	$1.5\times {10}^{-4}$	$3.7\times {10}^{-5}$	$2.1\times {10}^{-5}$
YC$12$-$1$ $\left(\theta =0\right)$		$5.8\times {10}^{-5}$	$3.4\times {10}^{-5}$
YC$12$-$1$ $\left(\theta =\pi \right)$		$7.1\times {10}^{-5}$	$4.4\times {10}^{-5}$

Equations

$$H=J_1\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{⟨⟨i,j⟩⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$\mathcal{S}=S_0(L_y)-B\sum _{n=1}^{N_f}\ln \left|2\sin \left[\frac{s}{2}\left(\theta -{\theta }_n^c\right)\right]\right|,$$

$$H_{MF}=-\sum _{⟨ij⟩}\sum _{\sigma =\uparrow /\downarrow }[(-1{)}^{s_{ij}}{f}_{i,\sigma }^{†}{f}_{j,\sigma }+h.c.].$$

$$2k_1=k+2\theta /L_y,k_2=k{L}_y/2.$$

Full text

Introduction.—Quantum spin liquids (QSLs) are exotic phases of matter which remain disordered due to quantum fluctuations, which in turn give rise to remarkable properties of fundamental importance, such as fractionalizations, gauge fluctuations, topology, and unconventional superconductivity Savary and Balents (2016); Zhou et al. (2017); Lee et al. (2006); Broholm et al. (2019). However, despite of a long-running quest, theoretical and experimentally relevant models for enigmatic QSLs are still limited and rare.

Historically, it has been proposed that geometric frustrations on the spin-$1/2$ triangular antiferromagnetic Heisenberg model (TAFM) could lead to a spin disordered ground state Anderson (1973). Although the nearest neighbor TAFM turns out to exhibit a $120^{\circ}$ magnetic order Sachdev (1992); Bernu et al. (1992); Capriotti et al. (1999); Zheng et al. (2006); White and Chernyshev (2007), the possibility of increasing the frustration by adding next-nearest-neighbor (NNN) interactions has captured much interest in the literature Manuel and Ceccatto (1999); Mishmash et al. (2013); Kaneko et al. (2014); Iqbal et al. (2016); Li et al. (2015); Hu et al. (2015); Zhu and White (2015); Saadatmand and McCulloch (2016); Szasz et al. (2018); Wietek and Läuchli (2017); Gong et al. (2017); Bauer and Fjærestad (2017); Yuste et al. (2018); Ferrari and Becca (2019); Jolicoeur et al. (1990); Chubukov and Jolicoeur (1992) for the $J_{1}$-$J_{2}$ TAFM

[TABLE]

where $\langle i,\,j\rangle$ and $\langle\langle i,\,j\rangle\rangle$ respectively denote NN and NNN bonds. So far, the general consensus is that an intermediate region ($0.07\lesssim J_{2}/J_{1}\lesssim 0.15$) without magnetic ordering Manuel and Ceccatto (1999); Mishmash et al. (2013); Kaneko et al. (2014); Iqbal et al. (2016); Li et al. (2015); Hu et al. (2015); Zhu and White (2015); Saadatmand and McCulloch (2016); Szasz et al. (2018); Wietek and Läuchli (2017); Gong et al. (2017); Bauer and Fjærestad (2017); Yuste et al. (2018); Ferrari and Becca (2019) is sandwiched between a stripe ordered phase ($J_{2}/J_{1}\gtrsim 0.15$) Jolicoeur et al. (1990); Chubukov and Jolicoeur (1992) and a $120^{\circ}$ magnetically ordered phase ($0.0\leq J_{2}/J_{1}\lesssim 0.07$) Sachdev (1992); Bernu et al. (1992); Capriotti et al. (1999); Zheng et al. (2006); White and Chernyshev (2007). However, the underlying physics and precise nature of this intermediate phase is under an intense debate. For instance, variational Monte Carlo simulations suggest a gapless $U(1)$ Dirac QSL Iqbal et al. (2016) as candidates for this intermediate phase. Density-matrix renormalization group (DMRG) calculations Hu et al. (2015); Zhu and White (2015); Saadatmand and McCulloch (2016); Szasz et al. (2018) found an indication of a gapped QSL as the nonmagnetic phase, while its internal structure (e.g. Z2, chiral) has yet to be determined. In addition, extensive exact diagonalization calculations fail to find evidence in support of either theory in the accessible system sizes Wietek and Läuchli (2017). Taken as whole, although a possible QSL phase has been identified on TAFM, the exact nature of this intermediate phase remains elusive.

It was shown that other experimental-relevant spin models on the triangular lattice also show spin liquid behavior which is continuously connected to the spin-liquid phase of $J_{1}$-$J_{2}$ TAFM model Zhu et al. (2018); Maksimov et al. (2019). Thus understanding the underlying physics in the $J_{1}$-$J_{2}$ TAFM, will give deep insight into a whole class of new triangular materials, for example, the recent synthesized Na-based chalcogenides Baenitz et al. (2018); Liu et al. (2018); Ding et al. (2019); Ranjith et al. (2019); Bordelon et al. (2019). In particular, the spin dynamics of NaYbO2 shows low-energy spectral weight accumulating at the $K$-point of the Brillouin zone Ding et al. (2019). So far it is unclear if these findings can be interpreted within spin liquid picture Song et al. (2019, 2018), which demonstrates the need for detailed theoretical predictions.

In this paper, we unveil the QSL nature of the triangular $J_{1}$-$J_{2}$ model by using large-scale DMRG simulations armed with recently developed state-of-the-art transfer matrix analysis He et al. (2017); Zauner et al. (2014). We find smoking-gun signatures of the $U(1)$ Dirac QSL (DSL), which consistently appear in 16 different geometries and/or system sizes (see Fig. 1 (a) for details). These signatures include: 1) momentum-dependent “excitation spectra”, extracted from the DMRG transfer matrix Zauner et al. (2014); He et al. (2017), which reveals gapless modes of the Dirac spin liquid showing recently predicted behavior of both fermion bilinear excitations as well as intricate monopoles Song et al. (2019, 2018); 2) strong dependence of the energy gap on twisted boundary conditions He et al. (2017); 3) universal entanglement entropy response under flux insertion Zhu et al. (2018). These evidences unambiguously show that the intermediate phase in TAFM is a gapless $U(1)$ DSL.

Properties of $U(1)$ DSL.—Let us begin with a brief review of properties of the $U(1)$ DSL on the triangular lattice Zhou and Wen (2002); Iqbal et al. (2016); Song et al. (2019, 2018). We begin with rewriting spin operator in terms of fractional fermionic spinons $\bm{f}=(f_{\uparrow},f_{\downarrow})^{T}$, $\bm{S}=\bm{f}^{\dagger}\bm{\sigma}\bm{f}$, where the partons $\bm{f}$ are coupled to a $U(1)$ dynamic gauge field due to the $U(1)$ redundancy. The $U(1)$ DSL can then be realized by putting spinons in a staggered $\pi$ flux mean-field ansatz, whose band structure will have two Dirac cones located at the $\pm Q$ points (valley, black dots in Fig. 1 (b)) of the Brillouin zone Iqbal et al. (2016); sm . The low energy physics of the $U(1)$ DSL is captured by $N_{f}=4$ QED3, namely there are four Dirac fermions (two from spins $\uparrow$/$\downarrow$ and two from valleys) coupled to a dynamic $U(1)$ gauge field. This $N_{f}=4$ QED3 theory may flow into a $2$+$1$D conformal field theory (CFT) in the infrared, therefore the $U(1)$ DSL is a critical/conformal phase Hermele et al. (2004, 2005, 2008), which is a close analog to the familiar spin-$1/2$ Heisenberg chain in 1+1D Giamarchi (2003); Brézin and Zinn-Justin (1989); Eggert and Affleck (1992). One effective way to detect the $U(1)$ DSL is to measure its gapless modes. It has been shown that the $U(1)$ DSL has two types of fundamental gapless modes, namely fermion (spinon) bilinears and monopoles (of the $U(1)$ gauge field) Song et al. (2019, 2018); Hermele et al. (2008); Alicea (2008). The fermion bilinears are “particle-hole” excitations of four Dirac fermions, while monopoles are instantons of the $U(1)$ gauge field. It is worth emphasizing that both fermion bilinears and monopoles are gauge invariant, which correspond to local operators such as spin $\bm{S}$, dimer operators $\bm{S}_{i}\cdot\bm{S}_{j}$, etc. Moreover, these critical operators have distinct quantum numbers (spins, momentum, angular momentum, etc.), enabling us to detect them directly.

There are in total $16$ fermion bilinears Hermele et al. (2005), which can be grouped into $1\oplus 15$, namely $SU(4)$ singlet and adjoint. They are distributed at different momenta in the $1^{\text{s}t}$ Brillouin zone (BZ) of the triangular lattice (seen in Fig. 1 (b)). The singlet bilinear is a spin singlet with zero momentum ($\Gamma$ point, black $\Circle$). The $15$ adjoint bilinears can be classified into three types Song et al. (2019): Type B1) $3$ time-reversal-even spin singlets with momenta located at three $M$ points of BZ: $M_{1}$ (violet $\varhexagon$), $M_{2}$ (blue $\varhexagon$) and $M_{3}$ (dark-cyan $\varhexagon$). Type B2) $3$ time-reversal-even spin triplets with zero momentum located at the $\Gamma$ points (black $\Circle$) of BZ. Type B3) $9$ time-reversal-odd spin triplets with momenta located at three $M$ points ($\varhexagon$). The quantum numbers of the monopoles in the $U(1)$ DSL remained elusive for decades until recently solved in Ref. Song et al. (2019, 2018). There are $6$ monopole operators (which are complex) of two types: Type M1) $3$ time-reversal-odd spin-triplets with momenta located at $K_{\pm}$ (red $\triangleleft$ and magenta $\triangleright$). Type M2) $3$ time-reversal-even spin-singlets with momenta located at the $X_{\pm}$ points (orange $\triangleleft$ and gold $\triangleright$). Physically, the condensation of spin-triplet monopoles will give the familiar $120^{\circ}$ non-colinear magnetic order, while the condensation of spin-singlet monopoles leads to valence bond solid such as $\sqrt{12}\times\sqrt{12}$ state Song et al. (2019). Later we will show that signatures of both the fermion bilinear and monopole operators have been measured in our DMRG simulations.

Method.—We use infinite-DMRG White (1992, 1993); McCulloch (2008) to simulate the $J_{1}$-$J_{2}$ TAFM wrapped on infinitely long cylinders. The evidence for DSL excitations is based on two major improvements in this study. First, we study different types of cylindrical geometries which correspond to different ways of wrapping a cylinder. As shown in Fig. 1(a), we define the YC$L_{y}$-$n$ cylinders by identifying the site $\bm{r}$ with the site $\bm{r}+L_{y}\bm{a}_{1}-n\bm{a}_{2}$ sm . For instance, the notation YC$8$-$1$ denotes a cylinder (C) with circumference of eight lattice spacing and a shift of one column in the y direction (Y) when connected periodically. Here, $\bm{a}_{1/2}$ are the triangular Bravais lattice primitive vectors, $L_{y}$ is the “circumference” of the cylinder, and $n$ amounts to a shift along the cylindrical direction. Simulating different geometries not only proves the observed DSL signatures are robust against finite size effect, but also serves as a nontrivial check as DSL on different geometries show qualitatively different behaviors (see Supplemental Materials sm ). Second, we carry out a numerical Aharonov-Bohm experiment by inserting flux $\theta$ in the cylinder, see Fig. 1 (a)). This is implemented by twisted boundary conditions, which modifies interactions crossing the boundary by a phase factor, e.g. $S^{+}_{i}S^{-}_{j}e^{i\theta}+S^{+}_{j}S^{-}_{i}e^{-i\theta}$ with a flux angle $\theta$ He et al. (2014). With the flux insertion, we can fully scan the momentum points in the Brillouin zone on a given geometry and are therefore not limited by finite-size energy gaps. Furthermore, certain physical quantities in DSL, such as the entanglement entropy, have a nontrivial response under flux insertion Zhu et al. (2018).

In the simulation it is important to keep the ground-state evolving adiabatically under the flux insertion. In most cases, the adiabatic flux insertion can be maintained, except very close to the Dirac cone (large flux $\theta$), where accurate infinite-DMRG simulation becomes very challenging due to the small gap and large entanglement of the state. Once adiabatic flux insertion fails at large $\theta$, the infinite-DMRG simulation may suddenly collapse to a competing state in other super-selection sectors of the ground-state sm or a symmetry broken state due to the instability of the gapless state He et al. (2017). We will not present the data of flux $\theta$ for which adiabatic flux insertion fails, as they do not reveal any direct information of spin-liquid ground state at zero flux.

Excitation gap.—Previous DMRG studies have found a considerably large spin gap in the $J_{1}$-$J_{2}$ TAFM Zhu and White (2015). However, this is not sufficient to exclude a DSL since on a cylinder the momentum is discrete, so the gapless Dirac point may be missed. The flux insertion, which effectively changes the quantized momentum of spinons, can make spinons hit Dirac point at specific values of flux $\theta$ He et al. (2017). By carefully studying the DSL ansatz incorporating the effect of emergent gauge fields sm , we find that DSL on different cylinder geometries YC$L_{y}$-$n$ have distinct $\theta$ dependence. If both $L_{y}$ and $n$ are even, spinons are gapless when $\theta=2\pi$ (Since spinons are fractional particles, the flux insertion has $4\pi$ periodicity). For all other three cases, spinons will be gapless at $\theta=\pi$ or $3\pi$.

Figure 2 shows the energy gap as a function of flux $\theta$. Although the gap is large at $\theta=0$ not , we find it significantly decreases as $\theta$ increases. The sensitivity of the energy gap is an indication of the gapless DSL: 1) for a gapped spin liquid the spin gap should have a small dependence (exponentially in $L_{y}$) on the flux; 2) finite flux drags the momentum lines toward the Dirac points, thus the gap monotonically reduces. Due to the small gap when Dirac points are approached, we are not able to maintain the adiabatic flux insertion when $\theta\sim 1.5\pi$ for YC$2n$-$2m$ cylinder, and $\theta\sim\pi$ for all other cylinders. There are also truncation effects from the finite bond dimensions $m$ in infinite-DMRG, which may explain that the YC$10$-$1$ gap appears larger than the YC$8$-$1$ gap in Fig. 3. We discuss results for different $m$ in the Supplemental Material sm . We also remark that the gap we measured may come from monopoles (rather than spinons), whose finite size effect is more subtle to analyze. The important message to take is, in all cases the gap systematically decreases as a function of $\theta$, and it is consistent with the theoretical expectation that the finite-size gap of spinons vanishes at i) $\theta=\pi$ for YC$8$-$1$ and YC$10$-$1$, ii) $\theta=2\pi$ for YC$8$-[math] and YC$10$-[math].

Correlation-length Spectrum.—While the energy gap is an important indication, the Dirac cone structure of the energy-momentum resolved spectrum will be a much stronger evidence of a DSL. So far the study of a large number of excited states has been very challenging, but fortunately recent seminal works Zauner et al. (2014); He et al. (2017) have uncovered a relationship between the energy spectrum and the spectrum of the transfer matrix in tensor-network formulation, which opens a window to the current problem.

The essence of this technique simply relies on a familiar fact: the information of excitations is encoded in the ground state, which can be decoded by measuring correlation functions of various operators. In infinite-DMRG simulations, the information of correlation functions of all operators can be straight-forwardly obtained through the eigenvalues of the transfer matrix Schollwöck (2011). Each eigenvalue takes the form $\lambda=e^{{ik}-1/\xi}$, where $\xi$ is the corresponding correlation length and $k$ is the momentum along the infinite-DMRG direction. The momentum around the cylinder can also be calculated from a revised transfer matrix He et al. (2017). The correlation lengths $\xi$ set an upper bound for excitation gaps $\Delta$ (up to a non-universal factor), and for a Lorentz invariant systems it holds that $\Delta\propto 1/\xi$. One can make this statement precise by an exact mapping from the infinite-DMRG transfer matrix to the partition function in the Euclidean path integral Zauner et al. (2014). In other words, if Lorentz (space-time rotation) symmetry is emergent in the system, the correlation-length spectrum precisely corresponds to the excitation spectrum of the Hamiltonian.

Fig. 3 shows the $S^{z}=1$ correlation-length spectrum of the $J_{1}$-$J_{2}$ TAFM at $J_{2}/J_{1}=0.12$. The left column shows the spectrum as a function of flux $\theta$. Since $\theta$ effectively changes the quantization of the momenta, we can then obtain the full dispersion relation as a function of $k_{1}$ and $k_{2}$ in the two right columns. For the cylinders YC$8$/$10$-[math] in Fig. 3 (a)/(b) the Dirac cones at the $M_{2}$ point $(k_{1},k_{2})=(0,\pi)$ can clearly be identified (blue $\varhexagon$), corresponding to fermion bilinear excitations of type B3 discussed above. In addition there are low lying monopole excitations close to the $K_{\pm}$ points $(k_{1},k_{2})=(-2\pi/3,2\pi/3)$ and $(2\pi/3,-2\pi/3)$ (red $\triangleleft$ and magenta $\triangleright$ respectively) of type M1. For the cylinders YC$8$/$10$-$1$ in Fig. 3 (c)/(d), we again find low lying excitations at the $M$ point $(2k_{1},k_{2})=(0,\pi)$ (blue $\varhexagon$) and $K$-points $(2k_{1},k_{2})=(2\pi/3,2\pi/3)$, $(-2\pi/3,-2\pi/3)$ (red $\triangleleft$ and magenta $\triangleright$). These observation of low lying excitations of fermion bilinear and monopole operators is clear evidence for a $U(1)$ DSL. We note that the lattice rotation symmetry $C_{6}$ is broken on the cylinder geometry, so its corresponding degeneracy is naturally split. In the Supplemental Material a total of 16 different cylinder geometries are analyzed, which consistently show the predicted $U(1)$ DSL excitations sm .

Entanglement Entropy.—Gapless spin liquids have nontrivial long-ranged quantum entanglement, in contrast to Landau ordered phases. We therefore also consider the bipartite entanglement entropy, ${\cal S}=-{\rm Tr}_{\rm sys}(\rho_{\rm sys}\ln\rho_{\rm sys})$, where the reduced density matrix $\rho_{\rm sys}={\rm Tr}_{\rm env}(|\Psi\rangle\langle\Psi|)$ for the half-cylinder “system” is constructed by ground-state wave function $|\Psi\rangle$ and traced over the degrees of freedom in the other half-cylinder “environment.” It was recently proposed that the entanglement entropy of $2$+$1$D CFT may have a universal response to an external Aharonov-Bohm flux Chen et al. (2017). In particular, for DSL Zhu et al. (2018), we expect

[TABLE]

where $S_{0}(L_{y})$ represents the area law part of entropy and $B$ is a prefactor which may or may not be universal. Other parameters are universal and can be determined by the underlying theory: $N_{f}$ accounts for the number of flavors of different Dirac spinons, $s=1/2$ is the fractional spin carried by Dirac spinons, and $\theta^{c}_{n}$ corresponds to the flux value at which the $n^{\textrm{th}}$ Dirac spinon becomes gapless. This scaling function (Eq. (2)) has been successful applied to identify the emergent DSL of the kagome antiferromagnet Zhu et al. (2018).

Fig. 4 shows the flux dependence of the entanglement entropy $\cal S$ at $J_{2}/J_{1}=0.12$, which has a strong dependence on flux $\theta$. This is a hallmark of low energy gapless excitations. In contrast, a fully gapped state would be largely insensitive to $\theta$. Moreover, as shown in Fig. 4, the dependence of $\cal S$ on $\theta$ can be fitted by the scaling function Eq. (2) with parameters $N_{f}=4$, $s=1/2$ and $\theta^{c}_{n}=2\pi$ for YC$2n$-[math] and $\theta^{c}_{n}=\pm\pi$ for YC$2n$-$1$. This agrees well with our theoretical expectation.

Summary and Discussion.—By combining large-scale DMRG simulations and recent analytical predictions, we study the intermediate spin liquid phase on the $J_{1}$-$J_{2}$ triangular antiferromagnetic Heisenberg model. Using flux insertion on different cylinder geometries we demonstrate that the energy gap of the spin liquid closes, and more importantly we find the low energy excitations of fermion bilinears and monopoles of Dirac spin liquid. The simultaneous appearance of fermion bilinears and monopoles is in favor of a Dirac spin liquid scenario, as opposed to the scenario of proximity to an ordered phase. Moreover, the entanglement entropy response under flux insertion agrees with a universal scaling law of the Dirac spin liquid. These findings strongly suggest that the intermediate phase of the $J_{1}$-$J_{2}$ TAFM is a gapless Dirac spin liquid.

Appendix A $U(1)$ Dirac spin liquid on

the triangular lattice

We fractionalize spin operators into fermionic partons, $\bm{S}=\bm{f}^{\dagger}\bm{\sigma}f$, $\bm{f}^{\dagger}=(f_{\uparrow}^{\dagger},f_{\downarrow}^{\dagger})$. The $U(1)$ Dirac spin liquid (DSL) can be described by a mean-field ansatz,

[TABLE]

$(-1)^{s_{ij}}$ is chosen to give a $\pi$/[math] flux on the up/down triangles. One can find a band structure of two Dirac cones in the $1^{\text{s}t}$ Brillouin zone. Therefore, the low-energy theory of this state is described by $N_{f}=4$ QED3.

Once we put $U(1)$ DSL on a finite-width and infinite-length cylinder, there are several subtle finite size effects, which can be understood in the mean-field theory. First, the cylinder may trap an extra flux $\phi=0$ or $\pi$ from the emergent $U(1)$ gauge field He et al. (2017). The value of $\phi$ will be energetically determined. Second, we note that there is an even-odd effect for the mean-field ansatz. As shown in Fig. AS1, cylinders with even and odd circumference $L_{y}$ should have different gauge configurations. Next, we solve band structure with an additional phase $e^{i\sigma\theta/2+i\phi}$ for the hoppings across the boundary in axis $\bm{a}_{1}$, where $\theta$ is the external Aharonov-Bohm flux inserted in the cylinder and the factor $\sigma/2$ is due to that $f_{\uparrow/\downarrow}$ carries $\pm 1/2$ spin. So on the YC$L_{y}$-$n$ cylinders, the properties of the band structure can be classified according to the even or oddness of $L_{y}$ and $n$. If both $L_{y}$ and $n$ are even, $\phi=\pi$ will be energetically favored, and the fermions become gapless when $\theta=2\pi$. For the other cases, the system behaves qualitatively similar no matter $\phi=0$ or $\pi$, and the fermions will become gapless at $\theta=\pi$ or $3\pi$, where only two (instead of four in the thermodynamic limit) Dirac fermions become gapless.

Appendix B Transfer matrix technique

In this section, we provide more details of the transfer matrix technique. In infinite-DMRG, we wrap a 2D lattice on a thin cylinder with a finite circumference $L_{y}$ but an infinite length $L_{x}$ and use the snake-chain matrix product state (MPS) to cover the cylinder as shown in Fig. BS2 (b) (for simplicity we view triangular lattice as a square lattice with a diagonal bond). Generically the snake-chain MPS geometry breaks the lattice translation symmetry along the $a_{1}$ direction, so we need to use $2n$ distinct matrices $A_{l}$ ($l=1$, $\cdots$, $2n$) in the MPS of YC$2n$-[math] cylinders (see Fig. BS2 (b), left panel), while $4n+2$ matrices for YC$(2n+1)$-[math] cylinder. For a special geometry, i.e. YC$L_{y}$-$1$, MPS recover two-site repeating structure, namely $A_{1}$-$A_{2}$, independent if $L_{y}$ is even or odd (see Fig. BS2 (b), right panel). On the other hand, the MPS is translation invariant and repeating along the direction $a_{2}$. Then, one can use the smallest repeating unit cell to define the transfer matrix (TM), as shown in Fig. BS2 (c). The eigenvalues $\lambda_{j}=e^{ik_{j}-\xi^{-1}_{j}}$ of TM–${\mathbf{T}}$ contain information of the correlation functions of all operators, which are further related to excitations of the system Zauner et al. (2014). Physically, each eigenvalue corresponds to one excitation mode of the system: $\xi_{j}$ gives the correlation length (or equivalently the inverse of gap), while $k_{j}$ gives the momentum along the infinite direction of the cylinder.

We can further extract the conserved quantum number of each excitation mode from TM. We note that the Schmidt basis (virtual index) of MPS has a well defined conserved quantum number ($Q^{\alpha}$), hence each eigenvector of TM–$\mathbf{T}$ shall have a quantum number $Q=Q^{\alpha}-Q^{\alpha^{\prime}}$ (Fig. BS2(c)). From this calculation we can get the correlation-length spectrum of different $S^{z}$ sectors. Below we discuss the method of calculating momentum $(k_{1},k_{2})$ of the correlation length spectrum.

For the YC$L_{y}$-0 geometry, $k_{2}$ is nothing but $k$ from the eigenvalues ($e^{ik-\xi^{-1}}$) of TM–$\mathbf{T}$, $k_{1}$ on the other hand requires a bit more work. Due to the snake covering the MPS does not have translational invariance along the compactified direction ($\bf a_{1}$). However, the Hamiltonian still has translational invariance along $\bf a_{1}$, hence the momentum $k_{1}$ is a conserved quantum number. In other words, $k_{1}$ is similar to other conserved quantum number (e.g. $S^{z}$) even though it is not encoded in MPS explicitly. To extract $k_{1}$ of each mode one needs to first obtain the conserved momentum $k_{1}$ of each Schmidt basis of the MPS using the mixed transfer matrix TM–$\mathbf{T}[T^{(y)}_{1}]$. The mixed TM–$\mathbf{T}[T^{(y)}_{1}]$ is pictorially defined in Fig. BS2(d), namely it is defined by translating the MPS by one site along the $\bf a_{1}$ direction. Due to the translation invariance (along $\bf a_{1}$) the dominant eigenvector of ${\mathbf{T}}[T^{(y)}_{1}]$ should be $\mathbf{V}_{\alpha,\alpha^{\prime}}=\delta_{\alpha,\alpha^{\prime}}e^{i\bf k^{\alpha}}$, with $\bf k^{\alpha}$ being the conserved momentum $k_{1}$ of each Schmidt basis of the MPS. At last, each mode has the momentum $k_{1}=\bf k^{\alpha}-\bf k^{\alpha^{\prime}}$, where $\alpha$ and $\alpha^{\prime}$ are the Schmidt basis of eigenvectors of TM–$\bf T$. In the simulation we inserted a finite flux $\theta$ in the cylinder, so the aforementioned momentum is further modified $(k_{1}+\theta/L_{y},k_{2})$ to obtain the final correlation spectrum presented in the paper.

For the special geometry YC$L_{y}$-$1$, the momenta $k_{1}$ and $k_{2}$ are intertwined together. A benefit is that the MPS with snake-geometry is actually invariant under two-sites translation. And the momentum $(k_{1},k_{2})$ can be obtained directly form the TM’s eigenvalues $e^{ik-1/\xi}$,

[TABLE]

Therefore, for the $-1$ geometry, one can still get $k_{1}$ and $k_{2}$, but $k_{1}$ has a $\pi$ ambiguity.

Appendix C Even/odd-sector states

It is well known that a topological ordered state has a number of topological degenerate groundstates once it is placed on a torus or an infinite cylinder (e.g. see He et al. (2014)). Physically all topological sectors can be understood as a state with gauge fluxes and/or gauge charge lines (i.e. anyon lines) threaded in the torus or cylinder. Similarly the $U(1)$ DSL also has different groundstate sectors, which we will call as super-selection sectors instead of topological sectors. The properties of super-selection sectors of $U(1)$ DSL are not well understood theoretically, and it will not be pursued here either. In this section we will try to clarify some confusion regarding the topological/super-selection sectors obtained in DMRG simulations.

For a spin liquid Hamiltonian, DMRG simulations may yield several different “groundstates” (They are local minima of the Hamiltonian). However, there is no a-priori knowledge that these “groundstates” are different topological/super-selection sectors of a spin liquid phase. One has to conduct a thorough study on these “groundstates”. For example, for a gapped topological spin liquid one needs to check if the modular matrix calculated from these “groundstates” agrees with the theoretical expectation for a topologically ordered phase. Otherwise one cannot exclude the possibility that different “groundstates” are the groundstates of different competing phases.

For the $J_{1}$-$J_{2}$ triangular spin liquid, one can obtain two “groundstates” on the YC$2n$-0 geometry. Previous works Hu et al. (2015); Zhu and White (2015); Saadatmand and McCulloch (2016); Gong et al. (2017) call these two “groundstates” the even and odd sectors, which pictorially corresponds to the number of valence bonds which cross a cut through the system. Their entanglement spectrum is one-fold or two-fold degenerate respectively, related to whether the $SO(3)$ spin rotation symmetry is realized projectively (two-fold) or not (one-fold). Numerically, we can get even/odd-sector states in the model without/with a pair of blank sites center-symmetrically located at two edge-columns during infinite-DMRG “warming-up” steps. For three cylindrical geometries (YC$2n$-$(2k+1)$, YC$(2n+1)$-$2k$, YC$(2n+1)$-$(2k+1)$) there is no distinction between the even and odd sector, as these two sectors are simply related by a translation symmetry. In contrast, on the YC$2n$-$2k$ cylinders the odd sector and even sector are different states, and we mostly focus on the odd sector for this geometry. The even sector on the YC$2n$-$2k$ cylinders has a smaller gap, which also shows signatures of $U(1)$ DSL in the correlation length spectrum.

We remark that a parallel work Gong et al. (2019) reported a large central charge $c=5$ for the even sector state on a finite YC$8$-[math] cylinder with a small length $L_{x}$. However, our data of the even sector state on YC$8$-[math] cylinder does not agree with their observation. We use $SU(2)$-infinite-DMRG and have simulated bond dimension up to $m^{\star}=8192$ (equivalent to $U(1)$-DMRG with $m=32014$). In Fig. CS3, we fit the central using three different methods Pollmann et al. (2009), namely the scaling of entropy with (a) the correlation length, (b) the $SU(2)$-DMRG bond dimension, and (c) the equivalent $U(1)$-DMRG bond dimension. The fitted central charge is much smaller than $c=5$, and it keeps on decreasing as bond dimension increases. The central charge is likely to eventually go to $c=0$ at the infinite bond dimension. So the observation of $c=5$ in Ref. Gong et al. (2019) might be a finite size effect: They only simulated a small cylinder ($L_{x}=16$), but the central charge is only well defined in the 1D limit with $L_{x}\gg L_{y}=8$.

We also note that there is a weak ($C_{6}$ breaking) nematic ordering in the spin liquid groundstates, namely a small difference of bond strengths in three inequal directions. The cylindrical geometry of DMRG breaks $C_{6}$ symmetry explicitly, so it is not surprising that a gapless $U(1)$ Dirac spin liquid weakly breaks $C_{6}$ Zhu et al. (2018). Second, the $C_{6}$ breaking is weak and it becomes even weaker as we increase the bond dimension (Tab. 1). The YC$6$-[math] cylinder has extremely tiny nematicity, which may come from the fact that DMRG simulation is fully converged for such a small $L_{y}$.

Appendix D Additional numerical data

D.1 Correlation length spectrum

We have in total simulated $16$ different geometries and/or system sizes. Similar to the $4$ clusters shown in main text, the fermion bilinear and monopole operators also show up in correlation length spectrum of all the other $12$ clusters. Here provide these data in this appendix.

D.1.1 YC$2n$-$(2k+1)$

Let us first look at the YC$2n$-$(2k+1)$ cylinder. We have already shown results of YC$8$-$1$ and YC$10$-$1$ cylinders in the main text, but for comparison we still plot them together with YC$6$-$1$ and YC$12$-$1$ ones in Fig. DS4. For YC$6$-$1$ in the panel (a), we can insert the flux adiabatically until the flux $\theta>0.3\pi$. Once a Dirac spin liquid is put on a small cylinder, it might have an instability by spontaneously generating massive terms during flux insertion. Therefore sometimes the numerical adiabatic change breaks down, but such finite size effect will be gone in the pure $2$+$1$D limit, so only data is shown where the ground state remains in the spin liquid state. This is consistent with our observation that, for a larger system size (i.e. YC$8$-$1$, YC$10$-$1$ and YC$12$-$1$), the adiabatic twist can be maintained even when $\theta=\pi$. The trend of the spectrum is similar to larger system sizes, although the information of $\theta\in(0,0.3\pi)$ is not enough for us to mark the type of its excitations. For other three cylinders, we also find one Dirac cone (low-lying excitations) is at a $M$ point $(2k_{1},k_{2})=(0,\pi)$ ($M_{2}$ or $M_{3}$ point) and the other is close to the $K_{\pm}$ point $(2k_{1},k_{2})=\pm(2\pi/3,2\pi/3)$. The former correspond to Fermion bilinears and the latter correspond to monopole operators.

Comparing the correlation length spectrum of various $L_{y}=8$, $10$ and $12$, it seems that the larger the system is, the higher the Dirac mode (blue $\varhexagon$) is. We think it is an artifact from the finite entanglement effect (bond dimensions) in infinite-DMRG simulations. For larger system sizes ($L_{y}$), the gapless modes suffers more severe truncation error, yielding a smaller correlation length at a given bond dimension. This artifact can be seen in Fig. DS5, where clearly the Dirac mode becomes lower as the bond dimension increases. On the other hand, to achieve the same accuracy (truncation error), the required bond dimension increases exponentially with the circumference of the cylinder. Table 2 shows the truncation error for different bond dimensions, system sizes and flux. For the large system size (YC$10$-$1$, YC$12$-$1$), $m=12288$ roughly gives comparable accuracy as $m=1024$ for YC$8$-$1$. Therefore, more care should be taken if one wants to compare the results between different system sizes.

D.1.2 YC$2n$-$2k$

Next we look at the YC$2n$-$2k$ cylinder. As we discussed previously, this class of cylinder behaves very different than YC$2n$-$(2k+1)$ discussed above. For the $U(1)$ Dirac spin liquid, the YC$2n$-$(2k+1)$ cylinder hits the gapless Dirac cone at $\theta=\pi$, while YC$2n$-$2k$ cylinder hits the gapless Dirac cone at $\theta=2\pi$. Our simulation on the YC$2n$-$2k$ cylinder is also consistent with this scenario. For YC$6$-[math], the adiabaticity of the twist can be maintained until $\theta=1.1\pi$, after which the system collapses to the other topological sector. For YC$8$-[math], YC$10$-[math] and YC$12$-[math] cylinders, adiabatic twist can persist until $\theta\approx 1.5\pi$.

Fig. DS6 shows the $S^{z}=1$ correlation length spectrum of the YC$6$-[math], YC$8$-[math], YC$10$-[math] and YC$12$-[math] cylinders at $J_{2}/J_{1}=0.12$. We find the lowest modes behave like the fermion bilinears and monopole operators of $U(1)$ DSL. Similar to the YC$2n$-$(2k+1)$ cylinder, the lowest modes show a linear dependence with the flux $\theta$. The lowest-lying fermion bilinear appears at $M_{2}=(0,\pi)$ point (labeled by $(k_{1},k_{2})$). We do not find spinon-pair excitations at $M_{1}$ and $M_{3}$, which can be understood by the explicitly broken rotation symmetry $C_{6}$ on a finite-width and infinite-length cylinder. Importantly, we also find that one additional branch appears nearby $K_{\pm}$ points $(k_{1},k_{2})=\pm(-2\pi/3,2\pi/3)$, signaling the expected monopole excitations.

D.1.3 YC$(2n+1)$-$(2k+1)$

Thirdly, the YC$(2n+1)$-$(2k+1)$ cylinder is basically the same as YC$2n$-$(2k+1)$. For the YC$(2n+1)$-$(2k+1)$ cylinder, we expect that spinons hit Dirac cones when $\theta=\pi$ or $3\pi$ independent of the emergent gauge flux $\phi=0$ or $\pi$. In Fig. DS7, we plot the $S^{z}=1$ correlation length spectrum as a function of the flux $\theta$, $2k_{1}$ and $2k_{2}$ respectively. Different from YC$2n$-$(2k+1)$ cylinder, the adiabatic twist cannot persist to $\pi$ due to the small gap and instability of the state as Dirac cone is approached. Furthermore, from the momentum-resolved spectrum we find Fermion bilinears appearing at a $M$ point $(2k_{1},2k_{2})=(0,0)$, and monopoles close to the $K_{\pm}$ points $(2k_{1},2k_{2})=\pm(2\pi/3,-2\pi/3)$.

D.1.4 YC$(2n+1)$-$2k$

In Fig. DS8, we plot the $S^{z}=1$ correlation length spectrum of YC$(2n+1)$-[math] cylinders as a function of the flux $\theta$, $k_{1}$ and $2k_{2}$ respectively. We find that the spinon-pair excitations appear at $(k_{1},2k_{2})=(\pi,0)$ ($M_{1}$ or $M_{3}$ point). Additionally, monopole excitations appears nearby $K_{\pm}$ $(k_{1},2k_{2})=\pm(-2\pi/3,-2\pi/3)$.

D.1.5 Other geometries

In previous studies, people use different geometries, such as “XC” geometry as defined in Ref. Zhu and White (2015) (left column of Fig. DS9 (a)) and “YC” geometry as defined in Ref. Hu et al. (2015) (left column of Fig. DS9 (b)). Therefore, we also analyze the adiabatic flux insertion for cylinders “XC”$8$-$1$ and “YC”$8$-$2$ as shown in Fig. DS9.

For the “XC”$8$-$1$ cylinder (equivalent to the YC$8$-$5$ cylinder), the spinon-pair excitations hit Dirac cones when $\theta=\pi$ independent of the emergent gauge flux $\phi=0$ or $\pi$. In Fig. DS9 (a), we find that its Dirac mode appears at a $M$ point $(4k_{1},k_{2})=(0,\pi)$ modulo $2\pi$. Monopole excitations are very close to $K_{\pm}$ points $(4k_{1},k_{2})=\pm(-2\pi/3,2\pi/3)$.

“YC” cylinders can be transformed to “XC” by a rotation $\pi/2$ in the $xy$-plane. For the “YC”$8$-$2$ cylinder (equivalent to the YC$9$-$2$ cylinder), the spinon-pair excitations hit a Dirac cone when $\theta=\pi$ independent of the emergent gauge flux $\phi=0$ or $\pi$ too. In Fig. DS9 (b), we find that its Dirac mode appears at a $M$ point $(2k_{1},4k_{2})=(0,0)$. Monopole excitations are very close to $K_{\pm}$ points $(2k_{1},4k_{2})=\pm(2\pi/3,2\pi/3)$ too.

D.2 Scaling behavior of entanglement entropy

In addition to Fig. 4, we show more data of entanglement entropy $\cal S$ for the cylinder YC$L_{y}$-[math] (Fig. DS10) and YC$L_{y}$-$1$ (Fig. DS11). Firstly, we find that Eq. (2) accurately fits the data around the minimal value of $\cal S$ for all the geometries and bond dimensions. Secondly, we notice that the fitting parameter $B$ shows strong dependence on bond dimension, system geometry, etc.. The strong dependence on bond dimension makes it hard to draw a conclusion on the question whether $B$ is a universal quantity or not. However, it is clear that the entanglement entropy always follows the universal scaling law conjectured for the $U(1)$ DSL.

D.3 Gaps

D.3.1 Gap measurement

We use an algorithm that combines infinite-DMRG and finite DMRG to calculate the spin gap in Fig. DS12. We first obtain a converged wave-function of an infinitely-long cylinder using infinite-DMRG “warm-up” steps. Then we insert a sector of cylinder consisting of $L_{y}\times L_{y}$ sites (red cylinder) into the middle of two half-chains. The left (L) and right (R) semi-infinite cylinder can be considered as environment (boundary conditions). We further do sweeps inside the small cylinder and get the lowest energy $E_{0}(S^{z}=1)$ in the $S^{z}=1$ sector and the energy of the $1^{\text{st}}$ excited state $E_{1}(S^{z}=0)$ in the $S^{z}=0$ sector. Finally, we obtain the gap $\Delta_{S^{z}=1}=E_{0}(S^{z}=1)-L^{2}_{y}e_{0}(S^{z}=0)$ and $\Delta_{S^{z}=0}=E_{1}(S^{z}=0)-L^{2}_{y}e_{0}(S^{z}=0)$ where we get the average energy per-site $e_{0}(S^{z}=0)$ during infinite-DMRG calculations.

D.3.2 Gap scaling

In addition to Fig. 2, we show more data for the gap in Fig. DS13. For large enough $L_{y}$, the spin gap of a DSL at $\theta=0$ decreases with the circumference size $L_{y}$ as $\Delta_{S^{z}=1}\sim v_{S^{z}=1}/L_{y}$, but numerically larger $L_{y}$ also has a larger trunctation error from finite bond dimension $m$, which tends to overestimate the spin gap in Fig. 2. A simple finite size scaling is therefore difficult, if $m$ sets a larger energy scale. We therefore demonstrate that the gap value becomes lower with increasing truncated bond dimension in Fig. DS13, which is consistent with a vanishing gap in the thermodynamic limit ($L_{y}\rightarrow\infty$).