Spinon Fermi surface in a cluster Mott insulator model on a triangular lattice and possible application to 1T-TaS$_2$

TL;DR

This paper models 1T-TaS2 as a spin liquid with a spinon Fermi surface on a triangular lattice, using a derived XXZ model and DMRG simulations, suggesting experimental detection methods.

Contribution

It introduces a novel effective spin model for 1T-TaS2 and demonstrates the existence of a gapless spin liquid with a spinon Fermi surface through large-scale simulations.

Findings

Identification of a gapless spin liquid phase with spinon Fermi surface

Static spin structure factor peaks on the '2k_F' surface

Proposal for experimental detection of the spinon Fermi surface

Abstract

1T-TaS$_2$ is a cluster Mott insulator on the triangular lattice with 13 Ta atoms forming a star of David cluster as the unit cell. We derive a two dimensional XXZ spin-1/2 model with four-spin ring exchange term to describe the effective low energy physics of a monolayer 1T-TaS$_2$, where the effective spin-1/2 degrees of freedom arises from the Kramers degenerate spin-orbital states on each star of David. A large scale density matrix renormalization group simulation is further performed on this effective model and we find a gapless spin liquid phase with spinon Fermi surface at moderate to large strength region of four-spin ring exchange term. All peaks in the static spin structure factor are found to be located on the "$2k_F$" surface of half-filled spinon on the triangular lattice. Experiments to detect the spinon Fermi surface phase in 1T-TaS$_2$ are discussed.

Tables

Table 1. Table S1: The tight binding parameters in the unit of eV.

${ϵ}_{z^2}$	${ϵ}_{xy}$	${ϵ}_{x^2-y^2}$	$t_0$	$t_{11}$	$t_{22}$	$t_2$	${\lambda }_0$
1.4052	1.3440	1.3440	-0.1046	0.2406	-0.5320	-0.3701	-0.3130

Table 2. Table S2: The basis for the neighboring star of David clusters

$n$	${\Psi }_n\in \mathscr{H}$
$1\sim 4$	$|i\alpha \uparrow ,j\alpha \uparrow ⟩$	$|i\beta \downarrow ,j\alpha \uparrow ⟩$	$|i\alpha \uparrow ,j\beta \downarrow ⟩$	$|i\beta \downarrow ,j\beta \downarrow ⟩$
$5\sim 8$	$|i\alpha \uparrow \beta \downarrow ,0⟩$	$|0,j\alpha \uparrow \beta \downarrow ⟩$	$|0,j\alpha \uparrow \alpha \downarrow ⟩$	$|i\alpha \uparrow \alpha \downarrow ,0⟩$
$9\sim 12$	$|i\beta \uparrow \beta \downarrow ,0⟩$	$|0,j\beta \uparrow \beta \downarrow ⟩$	$|0,j\alpha \uparrow \beta \uparrow ⟩$	$|i\alpha \uparrow \beta \uparrow ,0⟩$
$13\sim 14$	$|i\alpha \downarrow \beta \downarrow ,0⟩$	$|0,j\alpha \downarrow \beta \downarrow ⟩$

Table 3. Table S3: The necessary conditions for detecting various spin liquid (SL) phase in DMRG simulation.

type of SL		structure factor	real space decay	static spin susceptibility
gapless SL	spinon Fermi surface	$2k_F$ peaks	power law	finite
	Z2 Dirac SL	-	power law	linear in $T$
	U(1) Dirac SL	-	universal power law	linear in $T$
gapped SL	Z2 SL	-	exponential	exponentially small
	U(1) SL*	-	exponential	exponentially small
	chiral SL	chiral order	exponential	exponentially small
*not allowed in 2D

Equations

$${\tilde{\mathcal{H}}}_{\text{eff}}$$

$$+\left(S_j^x{S}_k^x+S_j^y{S}_k^y+(1+\gamma )S_j^z{S}_k^z\right)\left(S_i^x{S}_l^x+S_i^y{S}_l^y+(1+\gamma )S_i^z{S}_l^z\right)-\left({\mathbf{S}}_i\cdot {\mathbf{S}}_k\right)\left({\mathbf{S}}_j\cdot {\mathbf{S}}_l\right)].$$

$$H_0={\sigma }_0\otimes \left(\begin{array}{ccc}{f}_0& f_1& f_2\\ f_1& f_{11}& f_{12}\\ f_2& f_{12}& f_{22}\end{array}\right)+{\lambda }_0{\sigma }_z\otimes \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& i\\ 0& -i& 0\end{array}\right),$$

$$f_0$$

$$f_{11}$$

$$f_{22}$$

$$f_1$$

$$f_2$$

$$f_{12}$$

$${\Psi }_{\alpha }^{\uparrow }$$

$$+{\alpha }_3\left({\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{2\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{3\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{4\uparrow }+{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{5\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{6\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{7\uparrow }\right)$$

$$+{\alpha }_4\left({\psi }_{d_{z^2}}^{2\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{3\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{4\uparrow }+{\psi }_{d_{z^2}}^{5\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{6\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{7\uparrow }\right)$$

$$+{\alpha }_5\left({\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{8\uparrow }+{\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{9\uparrow }+{\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{10\uparrow }+{\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{11\uparrow }+{\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{12\uparrow }+{\psi }_{d_{x^2-y^2}+i{d}_{xy}}^{13\uparrow }\right)$$

$$+{\alpha }_6\left({\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{8\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{9\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{10\uparrow }+{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{11\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{12\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{x^2-y^2}-i{d}_{xy}}^{13\uparrow }\right)$$

$$+{\alpha }_7\left({\psi }_{d_{z^2}}^{8\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{9\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{10\uparrow }+{\psi }_{d_{z^2}}^{11\uparrow }+e^{-i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{12\uparrow }+e^{i\frac{2\pi }{3}}{\psi }_{d_{z^2}}^{13\uparrow }\right),$$

$$H=\left(\begin{array}{cccc}2tC+\mu & 2t_1{C}_{-}& 0& 0\\ 2t_1^{\ast }{C}_{+}& 2tC+\mu & 0& 0\\ 0& 0& 2tC+\mu & 2t_1{C}_{-}\\ 0& 0& 2t_1^{\ast }{C}_{+}& 2tC+\mu \end{array}\right)-\frac{1}{2}\lambda {\sigma }_z\otimes {\tau }_z,$$

$$C$$

$$C_{+}$$

$$C_{-}$$

$$\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_1,$$

$${\mathcal{H}}_0$$

$$+U_1{c}_{i\alpha \downarrow }^{†}{c}_{i\alpha \downarrow }{c}_{i\beta \downarrow }^{†}{c}_{i\beta \downarrow })+\sum _i\lambda \left(c_{i\beta \uparrow }^{†}{c}_{i\beta \uparrow }+c_{i\alpha \downarrow }^{†}{c}_{i\alpha \downarrow }\right),$$

$${\mathcal{H}}_1=\sum _{⟨i,j⟩,\sigma }\left(t{c}_{i\alpha \sigma }^{†}{c}_{j\alpha \sigma }+t{c}_{i\beta \sigma }^{†}{c}_{j\beta \sigma }+t_1{c}_{i\alpha \sigma }^{†}{c}_{j\beta \sigma }+t_1^{\ast }{c}_{i\beta \sigma }^{†}{c}_{j\alpha \sigma }\right)+h.c.,$$

$${\mathcal{H}}^{\left(2\right)}=\frac{4t^2}{U_1}\left(S_i^x{S}_j^y+S_i^y{S}_j^y\right)+4\left(\frac{t^2}{U_1}+\frac{t_1{t}_1^{\ast }}{U_0+\lambda }-\frac{t_1{t}_1^{\ast }}{U_1+\lambda }\right)S_i^z{S}_j^z-4\left(\frac{t^2}{U_1}+\frac{t_1{t}_1^{\ast }}{U_0+\lambda }+\frac{t_1{t}_1^{\ast }}{U_1+\lambda }\right),$$

$${\mathcal{H}}^{\left(2\right)}=\left(\begin{array}{cccc}-\frac{2t^2}{U_1}+\frac{2t_1{t}_1^{\ast }}{U_1+\lambda }-\frac{2t_1{t}_1^{\ast }}{U_0+\lambda }& \frac{2t^2}{U_1}& 0& 0\\ \frac{2t^2}{U_1}& -\frac{2t^2}{U_1}+\frac{2t_1{t}_1^{\ast }}{U_1+\lambda }-\frac{2t_1{t}_1^{\ast }}{U_0+\lambda }& 0& 0\\ 0& 0& \frac{2t_1{t}_1^{\ast }}{U_1+\lambda }+\frac{2t^2}{U_1}+U_1& \frac{2t^2}{U_1}\\ 0& 0& \frac{2t^2}{U_1}& \frac{2t_1{t}_1^{\ast }}{U_1+\lambda }+\frac{2t^2}{U_1}+U_1\end{array}\right).$$

$${\mathcal{H}}^{\left(2\right)}$$

$$\approx \tilde{\mathcal{H}}+\left[S,\tilde{\mathcal{H}}\right]+\frac{1}{2}\left[S,\left[S,{\tilde{\mathcal{H}}}_{\text{diag}}\right]\right],$$

$$\tilde{\mathcal{H}}=\left(\begin{array}{cccc}0& 0& t_{\alpha }& t_{\beta }\\ 0& 0& -t_{\beta }& -t_{\alpha }\\ t_{\alpha }& -t_{\beta }& U& 0\\ t_{\beta }& -t_{\alpha }& 0& U\end{array}\right),S=\frac{1}{U}\left(\begin{array}{cccc}0& 0& -t_{\alpha }& -t_{\beta }\\ 0& 0& t_{\beta }& t_{\alpha }\\ t_{\alpha }& -t_{\beta }& 0& 0\\ t_{\beta }& -t_{\alpha }& 0& 0\end{array}\right),$$

$$U$$

$$t_{\alpha }$$

$$t_{\beta }$$

$$\tilde{\mathcal{H}}=\sum _{⟨i,j⟩}\left(t_{\alpha }{\tilde{c}}_{i\alpha \uparrow }^{†}{\tilde{c}}_{j\alpha \uparrow }+t_{\beta }{\tilde{c}}_{i\beta \downarrow }^{†}{\tilde{c}}_{j\beta \downarrow }\right)+\sum _{i}U{\tilde{c}}_{i\alpha \uparrow }^{†}{\tilde{c}}_{i\alpha \uparrow }{\tilde{c}}_{i\beta \downarrow }^{†}{\tilde{c}}_{i\beta \downarrow }+h.c.,$$

$${\tilde{\mathcal{H}}}^{\left(4\right)}$$

$$+\left[t_{\alpha }{t}_{\beta }\left(S_j^x{S}_k^x+S_j^y{S}_k^y\right)+\frac{t_{\alpha }^2+t_{\beta }^2}{2}{S}_j^z{S}_k^z\right]\left[t_{\alpha }{t}_{\beta }\left(S_i^x{S}_l^x+S_i^y{S}_l^y\right)+\frac{t_{\alpha }^2+t_{\beta }^2}{2}{S}_i^z{S}_l^z\right]-t_{\alpha }^2{t}_{\beta }^2\left({\mathbf{S}}_i\cdot {\mathbf{S}}_k\right)\left({\mathbf{S}}_j\cdot {\mathbf{S}}_l\right)\}.$$

$${\tilde{\mathcal{H}}}_{\text{eff}}$$

$$+\left(S_j^x{S}_k^x+S_j^y{S}_k^y+(1+\gamma )S_j^z{S}_k^z\right)\left(S_i^x{S}_l^x+S_i^y{S}_l^y+(1+\gamma )S_i^z{S}_l^z\right)-\left({\mathbf{S}}_i\cdot {\mathbf{S}}_k\right)\left({\mathbf{S}}_j\cdot {\mathbf{S}}_l\right)].$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_{i^{\prime }}=\frac{1}{2}\left(S_i^{+}{S}_{i^{\prime }}^{-}+S_i^{-}{S}_{i^{\prime }}^{+}\right)+S_i^z{S}_{i^{\prime }}^z$$

$$D_b(i)D_b(i^{\prime })=\left(\frac{1}{2}\left(S_i^{+}{S}_{i+\delta }^{-}+S_i^{-}{S}_{i+\delta }^{+}\right)+S_i^z{S}_{i+\delta }^z\right)\left(\frac{1}{2}\left(S_{i^{\prime }}^{+}{S}_{i^{\prime }+\delta }^{-}+S_{i^{\prime }}^{-}{S}_{i^{\prime }+\delta }^{+}\right)+S_{i^{\prime }}^z{S}_{i^{\prime }+\delta }^z\right)$$

$$X_{△}=\frac{i}{2}\left(S_i^{+}{S}_j^{-}-S_i^{-}{S}_j^{+}\right)S_k^z+\frac{i}{2}\left(S_k^{+}{S}_i^{-}-S_k^{-}{S}_i^{+}\right)S_j^z+\frac{i}{2}\left(S_j^{+}{S}_k^{-}-S_j^{-}{S}_k^{+}\right)S_i^z$$

$$X_{△}{X}_{{△}^{\prime }}$$

Full text

Quantum spin liquid (QSL) was first proposed by P. W. Anderson in 1973 Anderson (1973). He argued that the ground state of spin-1/2 Heisenberg antiferromagnet on the triangular lattice is a random quantum superposition of singlets, the so called resonating valence bonds (RVB). Although the RVB state is not the true ground state of the triangular lattice spin-1/2 Heisenberg model, Anderson’s proposal has inspired a great deal of study this new "quantum liquid" state in frustrated magnetic systems. QSL is a highly entangled states and is very difficult to realize and characterize in experiments due to the lacking of an obvious order parameter and symmetry breaking. During the past forty years people only find a few QSL candidates, such as organic compounds $\kappa$-(BEDT-TTF)2Cu2(CN)3 Shimizu et al. (2003) and EtMe3Sb[Pd(dmit)2]2 Itou et al. (2007), herbertsmithite (ZnCu3(OH)6Cl2) Shores et al. (2005), Na4Ir3O8 Okamoto et al. (2007), YbMgGaO4 Li et al. (2015a, b); Shen et al. (2016); Paddison et al. (2017) and recently proposed 1T-TaS2 Law and Lee (2017); Klanjšek et al. (2017); Yu et al. (2017); Ribak et al. (2017), but still with many controversies in details. The ongoing efforts are either in the direction to explore new QSL candidate materials Feng et al. (2017); Ding et al. (2018), or push our theoretical and numerical understanding further.

It is well-known that, the geometrically frustration on kagome, pyrochlore and triangular lattices, or spin anisotropy such as Kitaev type interaction on a honeycomb lattice Kitaev (2006), play an important role to stabilize a QSL phase Balents (2010). On a kagome lattice, the isotropic nearest neighbor antiferromagnetic Heisenberg interaction is probably enough to result in QSL phase based on density matrix renormalization group (DMRG) Yan et al. (2011); Depenbrock et al. (2012); Jiang et al. (2012); Liao et al. (2017) or variational Monte Carlo (VMC) Iqbal et al. (2013) calculations, while on a triangular lattice, it is not the case. The ground state of Heisenberg model on the triangular lattice is the $120^{\circ}$-AFM state Miyake (1992); Chubukov et al. (1994); White and Chernyshev (2007). Thus to stabilize QSLs, more frustration, such as next neighbor frustrations Li et al. (2015c); Zhu and White (2015); Hu et al. (2015); Kaneko et al. (2014); Iqbal et al. (2016), anisotropic 111there is spatial anisotropy in organic compounds $\kappa$-(BEDT-TTF)2Cu2(CN)3 Shimizu et al. (2003) and EtMe3Sb[Pd(dmit)2]2 Itou et al. (2007) or high order exchange interactions is needed. The ring exchange terms become important for systems close to the insulating side of the Mott transition and it is suggested that the organics belongs to this case Motrunich (2005); Lee and Lee (2005). Exact diagonalization and variational study of the triangular lattice spin model with ring exchanges find a gapless QSL ground state with a spinon Fermi surface Motrunich (2005). Later DMRG simulation on two and four spin ladders and Gutzwiller variational wave functions calculation also find a similar QSL phase Sheng et al. (2009); Block et al. (2011).

1T-TaS2 was proposed to be a QSL candidate by two of us Law and Lee (2017). It has quasi-2D structure and each layer is made up of a triangular lattice with Ta atoms. It is recognized that 13-site clusters are formed with very narrow band near Fermi surface due to spin-orbit coupling (SOC) Fazekas and Tosatti (1979); Rossnagel and Smith (2006). A weak residual repulsion interaction is enough to open a Mott gap. Charge fluctuations induce high order exchange processes for the local moments if the system is close to the Mott transition (a weak Mott insulator). There are good reasons to expect this to be the case for 1T-TaS2 because it is the only insulator among all the CDW compounds and a related material 1T-TaSe2 is metallic. This motivates us to derive an effective spin model that include the effect of SOC and ring exchange. The geometric frustration and high order exchange interaction and spin anisotropy are new ingredients for the possible QSL physics in 1T-TaS2. In this paper we first derive a microscopic effective spin model including the anisotropy modified ring exchange interaction for this kind of material, and then perform the state-of-art large-scale DMRG simulation to explore ground state over quite a large range of parameter space. Our work will not only shed a new light in the understanding of ground state of 1T-TaS2, but will also push the limit of DMRG results for XXZ model with ring exchange on the triangular lattice, which is relevant to many other materials.

Effective spin model of 1T-TaS2 —  In 1T-TaS2, the Ta atoms form a planar triangular lattice sandwiched by S atoms in an octahedral coordination. The Ta layer and S layers have the ABC type stacking, which restores the global inversion symmetry for the crystal structure. As the temperature is lowered, 1T-TaS2 undergoes a series of charge-density wave phase (CDW) transition and eventually entering the commensurate CDW phase around 180K. This is the Mott insulating state Sipos et al. (2008) where the lattice is deformed into a superlattice with the unit cell of star of David, formed by 13 Ta clustered atoms seen from Fig. S1 (b) of supplementary materials (SM) sup . In the $\sqrt{13}\times\sqrt{13}$ star of David unit cell, the outer twelve Ta atoms have displacement toward the centered Ta atom, which strengthens the interatomic bonds inside the star of David and weakens others. As the early first principles bulk band structure calculation for 1T-TaS2 indicates that the Ta $5d$ orbitals are dominant in the conduction and valence bands Mattheiss (1973); Smith et al. (1985); Rossnagel and Smith (2006); Darancet et al. (2014); Qiao et al. (2017), the atomic SOC from $d_{x^{2}-y^{2}}$ and $d_{xy}$ orbitals is expected to modify the reconstructed band structure in the commensurate CDW phase. Importantly, the joint effect of lattice deformation and atomic SOC gives rise to the well isolated narrow band at the Fermi level, as is shown in Fig. S1 (a) of SM sup . As a result, in the presence of weak repulsive interaction, the 1T-TaS2 is susceptible to the Mott-Hubbard transition and turns out to be a Mott insulator.

In order to describe the Mott state in the 1T-TaS2, we consider a single star of David unit cell as a super-site, which is described by the intra-cluster tight binding Hamiltonian. Through numerical diagonalization, the Wannier orbitals localized inside the star of David with corresponding eigen-energies can be obtained in terms of the linear combination of atomic orbitals from the 13 Ta atoms. At the energy of the narrow band, it is found that the Wannier orbitals $\Psi^{\uparrow}_{\alpha}$ and $\Psi^{\downarrow}_{\beta}$ form the Kramers doublet while the Wannier orbitals $\Psi^{\downarrow}_{\alpha}$ and $\Psi^{\uparrow}_{\beta}$ are lifted in the energy due to the atomic SOC. Here the expressions for the two Wannier orbitals can be found in SM sup . Taking the two Wannier orbitals as the basis, we can construct a two-orbital Hubbard model with both inter-orbital and intra-orbital interactions for 1T-TaS2 sup . Since each star of David unit cell occupied with the single state $\Psi^{\uparrow}_{\alpha}$ or $\Psi^{\downarrow}_{\beta}$ would have the lowest energy, all other occupation states can be perturbatively dealt with through the Schrieffer-Wolff transformation. As a result, the effective XXZ spin model with the anisotropy modified ring exchange terms can be obtained as

[TABLE]

In the effective spin model, the $J$-term is a XXZ type nearest neighbor interaction, where $\gamma$ denotes spin anisotropy, which arises as the ratio between the inter-orbital and intra-orbital interaction deviates from one sup . Due to the atomic SOC, the effective spin model does not have the SU(2) spin rotational symmetry but preserves the U(1) rotation around the $z$ direction. Eq. (1) is of general interest as an effective spin Hamiltonian including SOC. Therefore, given the large SOC in 1T-TaS2, the smallness of $\gamma$ was not obvious a priori and required a demonstration. However in practice, it turns out that for 1T-TaS2 when the inter-orbital and intra-orbital interactions are in the same order, the anisotropy $\gamma$ remains smaller than 0.1 sup . In the large limit of atomic SOC, the anisotropy $\gamma$ will be further suppressed sup . In the rest of the paper we will mostly treat the case $\gamma$=0. The $K$-term is the four spin ring exchange term and is modified by the spin anisotropy. In general, the strength of $K/J$ depends on the ratio between the effective in-plane hopping and interaction. In the weak Mott insulating regime, the effective hopping and interaction are at the same scale, which is verified in several first principle calculations of 1T-TaS2 Darancet et al. (2014); Qiao et al. (2017), and then the strength of $K/J$ is of order one. The details on the derivation of the effective spin model and a comparing of parameters definition with earlier studied ring exchange model Motrunich (2005); Sheng et al. (2009); Block et al. (2011) can be found in SM sup .

For the spin model in Eq. (1), there are some well-known limits. (i) $K/J=0$, $\gamma=0$ case. In this case, we have a pure Heisenberg model on the triangular lattice and the ground state is the famous $120^{\circ}$-AFM state Miyake (1992); Chubukov et al. (1994); White and Chernyshev (2007). (ii) $K/J=0$,$\gamma\rightarrow\infty$ case. When $\gamma=\infty$, we have a pure Ising model on the triangular lattice. Due to the geometry frustration, the Ising spin does not order at zero temperature. As this paramagnetic state is highly degenerate, small perturbation may drive it to an ordered state via the order by disorder Villain et al. (1980). (iii) $K/J=\infty$, $\gamma=0$ case. In this case, we only have isotropic four-spin exchange terms. The ground state in the classical limit has been discussed in Ref. Kubo and Momoi (1997). As in real materials, $K$ is usually in the same order of $J$ or smaller, this case is less relevant.

Results —  For general values of $K/J$ and $\gamma$, the ground states are not known. To identify all possible ground states over a wide range of parameter space, we use DMRG to solve the effective spin model (1). The matrix product states (MPS) representation is used in our DMRG simulation. Due to the $U(1)$ spin rotational symmetry, the model has a total $S^{z}$ conservation, and all results are obtained in $S_{\text{tot}}^{z}=0$ sector 222Except the calculation with magnetic fields, which is performed with general DMRG. We use the cylindrical geometry with open boundary condition in the $x$ direction (See Fig. S4 of SM sup ). We use a bond dimension up to $m=5120$ and all results are obtained with a truncation error in or less than the order of $10^{-5}$.

To detect possible orders, we measure the spin-spin correlation $\langle\mathbf{S}_{i}\cdot\mathbf{S}_{i^{\prime}}\rangle$, the dimer-dimer correlation $\langle D_{b}(i)D_{b}(i^{\prime})\rangle-\langle D_{b}\rangle^{2}$ where $D_{b}(i)=\mathbf{S}_{i}\cdot\mathbf{S}_{i+\delta}$ (where $b=x,y,xy$ denotes $\delta=\hat{x}$, $\hat{y}$, or, $\hat{x}+\hat{y}$), the chirality-chirality correlation $\langle X_{\triangle}X_{\triangle^{\prime}}\rangle$ where $X_{\triangle}=\mathbf{S}_{i}\cdot\left(\mathbf{S}_{j}\times\mathbf{S}_{k}\right)$ (with $i,j,k\in\triangle$). The phase diagram in Fig. 1 shows the results that we obtain mainly from six wide systems ($L_{y}=6$) and confirm in eight wide systems ($L_{y}=8$). Our work extends previous results done on four wide systems  Sheng et al. (2009); Block et al. (2011). The details of finite size scaling to obtain the phase diagram can be found in SM sup . It turns out the small anisotropy ($\gamma\lesssim 0.1$) only shifts phase transition point slightly, thus we mainly focus on the isotropic case ($\gamma=0$). At small $K/J$, the ground state is the $120^{\circ}$-AFM state. In the intermediate value of $K/J$, a staggered valence bond solid (VBS) phase emerges. When $K/J>0.3$, we enter a QSL state with a spinon Fermi surface(SFS). This state was called a spin Bose metal in some earlier literature  Sheng et al. (2009); Block et al. (2011)

We here focus on the SFS phase. In the SFS phase, all structure factors, including the spin, dimer and chirality, shows no features of ordering. Thus we can rule out all spin, dimer and chirality orders in this phase. Second, the real space correlations of spin-spin and dimer-dimer show long range correlations and generally can be well fitted with power laws, as showed in Fig 2. The long range correlation is consistent with a gapless phase rather than a gapped one. In addition to a global power law decay in all these correlation functions, there are modulations and sign changes superposed over them. In the following, we will see the modulations in the spin correlator are actually the consequence of "$2k_{F}$" singularity due to the existence of a spinon Fermi surface. The modulation of the dimer correlator is discussed further in the SM sup . Thus the scatter of the data in Fig. 2 is due to these modulations and not due to numerical noise. As expected for a six-spin correlator, the chirality-chirality correlator shows much more rapid decay and can be fitted either with a large power law or an exponential. In addition, it exhibits sign changes. Thus we cannot interpret the chirality correlator in terms of the gauge flux correlator which is expected to have a power law decay with no sign changes in the asymptotic long distance limit. Apparently the system size is too small to reach that limit.

Another evidence for a gapless spin liquid comes from the finite static spin susceptibility. As shown in the SM sup , we apply a small magnetic filed $B$ along $z$-direction, measure magnetic moment density $M$ along the same direction and calculate the static spin susceptibility by $\chi=\frac{\partial M_{z}}{\partial B_{z}}$. To reduce the finite size effect, several twisted boundary conditions are considered. We find that the magnetic moment density and magnetic fields obey a linear behavior and we get a finite static spin susceptibility ($\chi\approx 0.22$ for $K/J=0.8$ where $J$ and $g\mu_{B}$ are set to one). Based on the value of the finite static spin susceptibility and using a half-filled free spinon band theory (static spin susceptibility is predicted to be $(g\mu_{B})^{2}N(0)/4$, where density of states at Fermi surface is about $N(0)\approx\frac{1}{1.65\pi t}$ with $t$ the spinon hopping parameter), we can estimate the spinon hopping $t$ to be $\sim 0.22J$. This estimate is not precise because we expect Fermi liquid corrections to the free spinon expression for the susceptibility. The important point is that the finite spin susceptibility is consistent with a gapless spin liquid with a spinon Fermi surface.

In order to confirm the SFS state, we look for $2k_{F}$ peaks in the spin structure factor $S(\mathbf{k})$. To simplify the discussion, we write $\mathbf{k}$ in the basis of reciprocal lattice primitive vectors $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$, namely, $\mathbf{k}=k_{1}\mathbf{b}_{1}+k_{2}\mathbf{b}_{2}$. We analyze the spin structure factor $S(\mathbf{k})$ for each fixed $k_{2}$ line, and pick out all peaks, as is shown in Fig 3. The positions of all those peaks are plotted in Fig. 4. Here results of both $L_{y}=6$ and $L_{y}=8$ wide systems are plotted together. All the points are located on the $2k_{F}$ surface of a half-filled spinon on the triangular lattice, within finite size resolution. This strong signature of the existence of spinon Fermi surface is a definitive evidence for SFS phase.

We point out that there is an unexpected feature of the spin structure factor, namely there is no peak at the $2k_{F}$ position along the $\mathbf{b}_{1}$ line. This is seen in Fig. 3(d) ($k_{2}=0$ curve) and agrees with earlier results on four wide systems Block et al. (2011). A possible explanation is that there are residual antiferromagnetic spin fluctuations associated with the $120^{\circ}$-spin order located at the K point $(\frac{1}{3},\frac{1}{3})$. These may connect the Fermi surface crossings near $(\frac{1}{3},0)$ and $(0,-\frac{1}{3}$) and suppress the density of states there. An even more interesting possibility is that a gap is open along the $\Gamma$ (0,0) to M ($\frac{1}{2},0$) directions due to spinon pairing. The Amperean scenario proposed by Lee et al Lee et al. (2007) will create a set of gaps separated by spinon Fermi arcs. In this case the spin liquid will belong to the class of Z2 QSL with spinon Fermi surface. Further discussion of this possibility and its relation to the structure of the dimer correlator can be found in the SM sup .

The partially filled spinon Fermi surface have an infinite number of gapless modes, but on a finite width cylinder, there are only finite numbers of them. A fitting of central charge $c$ has been done by analyzing the entanglement entropy calculated with DMRG, which measures the number of gapless modes and is found to increase roughly linearly with width of cylinder, consistent with SFS phase. As discussed in the SM sup ,the central charge can in principle distinguish between the U(1) and Z2 QSL. Unfortunately in the SFS phase the entanglement entropy converges slowly when we increase the number of states in DMRG, and it is not possible to extract the central charge unambiguously. On the other hand the positions of peaks as well as the spin structure factor converge faster, as showed in Fig 3.

Another interesting finding is the VBS phase. As shown in the phase diagram (Fig. 1), at intermediate $K/J$, there is a VBS phase that is revealed by the dimer-dimer correlation. We clearly observe y-direction dimer structure with sharp peak at $M$ point ( $(\frac{1}{2},\frac{1}{2})$ in the basis of reciprocal lattice primitive vectors). Detailed results can be found in SM sup . From the continuous dependence of energy on $K/J$ in our DMRG simulation on finite sizes, it is very likely that the transition from AFM to VBS or VBS to SFS is a continuous one, where the AFM to VBS transition may generate a deconfined quantum critical point (DQCP) Senthil et al. (2004) and worth a further study.

As we have previously mentioned, the anisotropy in 1T-TaS2 is actually small. We have explored the effect of this small anisotropy to the phase diagram, and found that the small anisotropy moderately suppresses the region of VBS phase and thus helps stabilize the QSL phase.

Discussion and Conclusions —  We have derived the effective spin model for 1T-TaS2. The effective model is essentially an XXZ model with a four-spin ring exchange term on the triangular lattice. We have demonstrated that the effective model supports a SFS ground state at moderate to large $K/J$ regions. The SFS that we find is a gapless QSL phase with a spinon Fermi surface. We clearly observe the singular wavevectors in the spin structure factor from the "$2k_{F}$" surface that can be described by half-filled spinon with a uniform hopping on the triangular lattice.

Finally we remark on the applicability of these results to 1T-TaS2. Our monolayer model is directly applicable to mono-layer 1T-TaS2 which can be grown by MBE. For bulk samples, if the interlayer coupling is weaker than the intra-layer exchange $J$, it is possible that the ground state is made up of weakly coupled layers of SFS. It will be of great interest to use neutron scattering to look for the $2k_{F}$ peaks in the static spin structure factor $S(\mathbf{k})$. Furthermore, the SFS state is expected to have low energy excitations concentrated around $2k_{F}$. Gapless excitations are also expected for $k<2k_{F}$ and to extend to an energy scale which is a fraction of the spinon bandwidth. This is seen in calculations based on the free spinon model Li et al. (2017). The SFS is expected to have finite spin susceptibility. Experimental there is indeed a residual temperature independent susceptibility Klanjšek et al. (2017) but it is not known how much of it is spin or orbital in origin. We also expect a linear term $\gamma$ in the specific heat. This is seen experimentally but it is suppressed by a magnetic field Ribak et al. (2017). This contribution has been explained as mainly due to local moment Kimchi et al. (2018). If we use the large magnetic field limit to extract an intrinsic $\gamma$, we find a value of about $0.1\text{mJ/K}^{2}$ per mole or $1.3\text{mJ/K}^{2}$ per mole of star of David cluster. This is about a factor of 10 smaller than what is observed in the organics Shimizu et al. (2003); Itou et al. (2007) , suggesting that the exchange energy scale is a factor of 10 larger. Finally the definitive evidence for a spinon Fermi surface will be the observation of a linear $T$ term in the thermal conductivity. While a null result was found earlier Yu et al. (2017), such a linear T term has recently been reported Murayama et al. (2018).

An exciting future avenue is the doping of the weak Mott insulator. Doping of the spinon Fermi surface spin liquid will likely lead to a correlation driven superconductor Anderson (1987). As discussed earlier Law and Lee (2017), the interesting regime is a doping concentration of about $10\%$ per cluster which is less than $1\%$ per Ta. Carrier localization is a serious challenge at such a low doping level and gating of atomically thin samples may be the preferred method. The recent success of inducing superconductivity by doping the weak Mott insulator in twisted bilayer graphene is certainly encouraging Cao et al. (2018). The 1T-TaS2 system has the advantage of having a much higher temperature scale compared with graphene and the organics.

Acknowledgments —  X.Y. Xu thanks the discussion with E.M. Stoudenmire and Donna Sheng. The calculations are performed using the ITensor C++ library (version 2.1.1). W.-Y. He, X.Y. Xu and K.T.L. thank the support of HKRGC through HKUST3/CRF/13G and C6026-16W. W.-Y. He also acknowledges the support of Hong Kong PhD Fellowship. G.C. is supported by the Ministry of Science and Technology of China with the Grant No.2016YFA0301001 and the Thousand-Youth-Talent Program of China. K.T.L. is further supported by the Croucher Foundation and the Dr Tai-chin Lo Foundation. P.A.L. acknowledges support by the US Department of Energy, Basic Energy Sciences under grant DE-FG02-03-ER46076. He also thanks the hospitality of the Institute for Advanced Studies at the Hong Kong University of Science and Technology. The simulation is performed at Tianhe-2 platform at the National Supercomputer Center in Guangzhou.

I Effective tight binding Hamiltonian for the star of David superlattice

In TaS2, at the Fermi level, electrons are mainly from the 5d orbitals of the Ta atoms. Before lattice deformation, the $d_{z^{2}}$, $d_{xy}$ and $d_{x^{2}-y^{2}}$ orbitals from the Ta atoms are considered to construct the tight binding Hamiltonian as follows

[TABLE]

where the second term represents the atomic SOC. The hopping terms are

[TABLE]

with $a$ the lattice constant for the triangular lattice before deformation. All the tight binding parameters can be found in Table.S1 Smith et al. (1985); Rossnagel and Smith (2006). In the low temperature, the 1T-TaS2 undergoes the commensurate charge density wave transition and form the triangular superlattice with $\sqrt{13}\times\sqrt{13}$ star of David unit cell, as is seen in Fig. S1 (b). In the star of David unit cell, as the outer twelve Ta atoms have displacement toward the centered Ta atom, the interatomic distance from the inner surrounding Ta atoms (labeled by 2-7) and outer surrounding Ta atoms (labeled by 8-13) to the centered Ta atom 1 shrink by $6.4\%$ and $4.4\%$ respectively. Since the hopping integral is scaled as $R^{-5}$ with $R$ the Ta-Ta distance Heine (1967), the atomic bonds within the star of David are therefore strengthened in the intra-cluster Hamiltonian $H_{\textrm{intra}}$ while other bonds are weakened in the inter-cluster Hamiltonian $H_{\textrm{inter}}$. As a result, after the lattice deformation, it is found in the reconstructed band structure that a narrow band is well isolated at the Fermi level, as is shown in Fig. S1 (a).

We first regard a single star of David cluster as a super-atom and consider only the intra-cluster tight binding Hamiltonian $H_{\textrm{intra}}$. Through numerical diagonalization, at the energy of the narrow band the Kramers doublet $\Psi^{\uparrow}_{\alpha}$, $\Psi^{\downarrow}_{\beta}$ are found to be the Wannier function that compose the narrow band. The Wannier orbitals can be expressed in terms of the linear combination of atomic orbitals from the 13 Ta atoms with the coefficients $\alpha_{1},\alpha_{2},...,\alpha_{7}$

[TABLE]

and $\Psi^{\downarrow}_{\beta}=i\sigma_{y}\Psi^{\uparrow\ast}_{\alpha}$ due to the time reversal symmetry, while the Wannier orbitals $\Psi^{\downarrow}_{\alpha}$, $\Psi^{\uparrow}_{\beta}$ are lifted in energy due to the atomic SOC. With the Wannier orbitals for the star of David super-atom, in the basis of $\left[\Psi^{\uparrow}_{\alpha},\Psi^{\uparrow}_{\beta},\Psi^{\downarrow}_{\alpha},\Psi^{\downarrow}_{\beta}\right]^{\textrm{T}}$ the minimum two orbital Hamiltonian can be constructed as

[TABLE]

with $t=-0.0072\textrm{eV},\ t_{1}=\left(-0.015+0.0068i\right)\textrm{eV},\ \mu=0.2897\textrm{eV},\ \lambda=0.5622\textrm{eV}$, and $\sigma$, $\tau$ for the spin, orbital space respectively. The matrix element for $H$ is obtained through $H_{i,j}=\langle\Psi_{i}\rvert H_{\textrm{intra}}+H_{\textrm{inter}}\lvert\Psi_{j}\rangle$. The energy dispersion for the narrow band from the minimum two orbital model is plotted as the green dashed line in Fig.S1 for comparison. Here the basis functions in Eq. S9 read

[TABLE]

II XXZ model from the two-orbital Hubbard model

With the two Wannier orbitals for the star of David super-atom, considering both the intra-orbital interaction $U$ and inter-orbital interaction $U_{1}$, we can construct the two-orbital Hubbard model as

[TABLE]

with

[TABLE]

and

[TABLE]

where $c^{\left(\dagger\right)}_{i,\alpha\left(\beta\right)\sigma}$ is the annihilation (creation) operator that annihilates (creates) an $\alpha$ ($\beta$) Wannier orbital state with spin state $\sigma$ at the site $i$. In order to get the effective spin model, we first analyze the neighboring star of David clusters at the site $i$, $j$. At quarter filling for the two-orbital Hubbard model, the states with single occupancy, listed as the first row in the table.S2, form the subspace with low energy. Through kinetic hopping, these states can transfer to the excited states listed in the remaining part of table.S2. As a result, writing the two-orbital Hubbard model in the basis from table.S2, we can perturbatively deal with the excited states through the second order Schrieffer-Wolff transformation and obtain the effective Hamiltonian in the low energy subspace

[TABLE]

which is the XXZ model. For the neighboring star of David clusters, the spin exchange mechanism is schematically shown in Fig. S2. As the two Wannier orbitals $\Psi_{\alpha}$ and $\Psi_{\beta}$ are both involved due to the atomic SOC, the spin exchange can go through either the intra-orbital hopping or the inter-orbital hopping. In the process in Fig. S2 (a), the exchange through the intra-orbital hopping gives rise to the Heisenberg model as usual. In the process in Fig. S2 (b) and (c), the virtual hopping between the low energy state and the excited state is inter-orbital type and the low energy state cannot be changed, so the inter-orbital type virtual hopping results in the $S^{z}$ terms in the spin model. Since the intermediate states in the process (b) and (c) have intra-orbital interaction and inter-orbital interaction as the excitation energy respectively, the two process cannot be canceled out and eventually generate the anisotropy in the $S^{z}$ term. As a result, considering the neighboring star of David clusters we obtain the XXZ spin model.

III Extrapolation to the anisotropy modified ring exchange term

In the neighboring two clusters analysis, as the anisotropic term in the XXZ model arises from the interplay of inter-orbital hopping, inter-orbital interaction and intra-orbital interaction, we would like to explore the possibility to reduce the two-orbital complexity to an effective one orbital Hubbard model with pseudo-spin states and simplify the whole system. We notice that including both the low energy states and excited states, the effective Hamiltonian $\mathcal{H}^{\left(2\right)}$ from the second order Schrieffer-Wolff transformation in the basis of $c^{\dagger}_{j\alpha\uparrow}c^{\dagger}_{i\beta\downarrow}\lvert 0\rangle$, $c^{\dagger}_{j\beta\downarrow}c^{\dagger}_{i\alpha\uparrow}\lvert 0\rangle$, $c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\beta\downarrow}\lvert 0\rangle$, $c^{\dagger}_{j\alpha\uparrow}c^{\dagger}_{j\beta\downarrow}\lvert 0\rangle$ reads

[TABLE]

We take $c^{\dagger}_{\alpha\uparrow}\lvert 0\rangle$ and $c^{\dagger}_{\beta\downarrow}\lvert 0\rangle$ as the pseudo-spin states and construct an effective one orbital Hubbard model $\tilde{\mathcal{H}}$ so that under Schrieffer-Wolff transformation the effective Hamiltonian $\mathcal{H}^{\left(2\right)}$ is recovered

[TABLE]

where $\tilde{\mathcal{H}}_{\textrm{diag}}$ is the diagonal part of the Hamiltonian matrix $\tilde{\mathcal{H}}$ and

[TABLE]

with the modified on-site interaction and kinetic hopping

[TABLE]

We write the effective one orbital Hubbard model $\tilde{\mathcal{H}}$ in the form of annihilation and creation operator

[TABLE]

and then apply the Schrieffer-Wolff transformation to generate the effective spin model MacDonald et al. (1988). For the neighboring two star of David clusters, we recover the XXZ model $\tilde{\mathcal{H}}^{\left(2\right)}=\mathcal{H}^{\left(2\right)}$. For a diamond cluster of star of Davids, we get the anisotropy modified ring exchange term

[TABLE]

Eventually the general form for the effective spin model of the 1T-TaS2 is derived as

[TABLE]

Here $\gamma=\frac{(t_{\alpha}-t_{\beta})^{2}}{2t_{\alpha}t_{\beta}}$ denotes the anisotropy of the spin model, and the anisotropy $\gamma$ can be estimated with Eq. S5 and Eq. S6 as is shown in Fig. S3 (a). In the estimation, we fix the $U_{1}+U_{0}=200\textrm{meV}$, which is about the Mott gap measured in the tunneling spectroscopy Qiao et al. (2017), and vary the ratio between the inter-orbital interaction $U_{1}$ and the intra-orbital interaction $U_{0}$ to obtain the evolution of the anisotropy $\gamma$. It can be seen that the anisotropy $\gamma$ increases slowly with the increase of $\frac{U_{1}}{U_{0}}$ and $\gamma$ reaches 0.1 when the $U_{1}=2U_{0}$. If we further increase the ratio $\frac{U_{1}}{U_{0}}\rightarrow\infty$, the anisotropy will reach a constant value around 0.4 and converge. In the limit of $\frac{U_{1}}{U_{0}}\rightarrow\infty$, the intra-orbital interaction is eliminated and the anisotropy is only affected by the atomic SOC $\lambda$. In the two orbital Hubbard model, as the gap between the two bands increases with the increase of atomic SOC, it regress to the single orbital Hubbard model and the anisotropy $\gamma$ will approach [math] in large atomic SOC, as is shown in Fig. S3 (b).

IV DMRG simulation setup

As we already discussed in the main text, several limits of effective spin Hamiltonian (S23) have been widely explored, but for general parameter region, especially in the situation strongly related to real materials 1T-TaS2 where $\gamma\lesssim 0.1$ and with a meaningful large of $K/J$, the physics is still unknown. To identify possible ground state around this region, we use DMRG to simulate the effective spin Hamiltonian (S23). The MPS representation is used in our DMRG simulation and calculations are performed using the ITensor C++ library (version 2.1.1), http://itensor.org/. The model has total $S^{z}$ conservation, and all results are got in $S_{\text{tot}}^{z}=0$ sector if it is not specified. We use the cylinder geometry with open boundary condition in $x$ direction, as showed in Fig. S4.

V Determine phase diagram

V.1 Identify different phases with various correlation functions

We will mainly focus on $\gamma=0$ as in 1T-TaS2, anisotropy $\gamma$ is small enough to be neglected. To identify different phases, various correlation functions are measured, including spin correlation, dimer correlation and chirality correlation.

Spin correlation $C_{S}(i,i^{\prime})=\langle\mathbf{S}_{i}\cdot\mathbf{S}_{i^{\prime}}\rangle$ with

[TABLE]

Dimer correlation $C_{D_{b}}(i,i^{\prime})=\langle D_{b}(i)D_{b}(i^{\prime})\rangle-\langle D_{b}\rangle^{2}$, where the dimer $D_{b}=\mathbf{S}_{i}\cdot\mathbf{S}_{i+\delta}$ can have three types labeled by $b=x,y,xy$ due to the three possible directions of bonds, namely $\delta=\hat{x}$, $\delta=\hat{y}$ or $\delta=\hat{x}+\hat{y}$.

[TABLE]

Chirality correlation $C_{X}(\triangle,\triangle^{\prime})=\langle X_{\triangle}X_{\triangle^{\prime}}\rangle$, where the spin chirality $X_{\triangle}$ is defined as $X_{\triangle}=\mathbf{S}_{i}\cdot\left(\mathbf{S}_{j}\times\mathbf{S}_{k}\right)$ (with $i,j,k\in\triangle$).

[TABLE]

[TABLE]

After the DMRG energy optimization, we will get the ground state wavefunction. Measuring above correlations is tedious and slow for large system size and large bond dimension. We need to take advantage of the mixed-canonical matrix product state (matrices to the left are left-normalized, to the right are right-normalized.) and classify the correlations with different site indexes into several classes and measure them separately. With this trick, the measurement is highly accelerated and thus we are able to measure all the correlation mentioned above with limited computation resources we have.

When we get the real space correlations, a fast way to check whether there is any order is to first perform the Fourier transformation to the momentum space and get the structure factor, and then check whether there is any peak in the structure factor. This first step usually will help us identify whether there exists certain ordered phase, like spin order, bond order, etc. If we do not see any peak (ignore finite size broadening), then we will roughly rule out any orders which result in peaks in the corresponding structure factor and it will be more like a QSL phase, although more evidences are needed.

One easily obtainable evidence is the functional form of the real space decay of different correlations. The real space decay of the spin correlation function will reflect the spin gap information, a short correlation length is consistent with a spin gap, while a decay of universal power law $\sim r^{-4}$ is consistent with (algebraic) U(1) Dirac spin liquid which is gapless. Actually, for algebraic spin liquid as a critical phase, the dimer correlation will also follow the same power law decay. For other correlations, we can also get information from the real space decay. For example, the short correlation length of dimer correlations is consistent with a singlet gap, while the long range chiral correlation is consistent with time-reversal broken chiral spin liquid phase. In addition, detailed analysis of momentum space structure will also provide further evidences. As discussed in Refs Sheng et al. (2009); Block et al. (2011), a detailed analysis of the $2k_{F}$ peaks in structure factor should help to identify the gapless QSL with spinon Fermi surface.

V.2 Distinguish gapped and gapless phase with static spin susceptibility

Another way to distinguish gapped and gapless phase is to calculate static spin susceptibility. Here we can apply small magnetic fields along $z$-direction and measure the magnetic moment density $M$. As the $z$-direction magnetic field will not break the total $S^{z}$ conservation, the changing of magnetic field does not make the common eigenstates of field term and zero-field Hamiltonian evolve Savary and Balents (2017), so the static spin susceptibility is directly related to the density of excitations. The finite static spin susceptibility will correspond to a gapless phase while a exponentially small static spin susceptibility will correspond to a gapped phase. For gapless Dirac phase, the static spin susceptibility will obey a linear $T$ dependence and will not be finite at zero temperature.

In Table S3, we list almost all possible types of spin liquids and their corresponding necessary conditions which can be checked in DMRG simulation as we discussed in above.

VI phase diagram

The phase diagram is already showed in Fig. 2 of main text.

VI.1 AFM phase

The AFM phase has the in-plane $120^{\circ}$-spin order with $Q$ vector at the corner $K$ point ( $(\frac{1}{3},\frac{1}{3})$ in the basis of reciprocal lattice primitive vectors). Fig. S5 shows spin structure factor at $K/J=0.08$ with $L_{x}=12$ and $L_{y}=6$. The finite size scaling of this spin order is showed in Fig. S6 only for $L_{y}=6$ systems, while for 2D finite size scaling with fixing $L_{x}$ and $L_{y}$ ratio is not possible, as we need both $L_{x}$ and $L_{y}$ to be multiples of 6 and to have three points in $y$-direction $L_{y}=6$, $12$ and $18$ is out of current limit of computation power. The $120^{\circ}$-spin order is prominent in small $K/J$ region (for example $K/J=0.08$ in the Fig. S6), and absent in larger $K/J$ region (for example $K/J=0.24$ and $0.40$ in the Fig. S6).

VI.2 VBS phase

In the intermediate $K/J$ region we get VBS phase, which forms orders of bond dimer, also called dimerized phase. As the original model has the threefold rotational symmetry, there is no difference among $\hat{x}$, $\hat{y}$ or $\hat{x}+\hat{y}$ direction dimerized phase. However, in finite size DMRG calculation with cylinder geometry we have broken this symmetry and lift the degeneracy, that is why we only find $y$-direction dimerized phase survive as the dominant. Fig. S7 shows the structure factor of $y$-direction dimer correlation for $K/J$=0.24 with $L_{x}=12$, $L_{y}=6$. The finite size scaling with fixing $L_{x}$ and $L_{y}$ ratio $L_{x}/L_{y}=3$, is showed in Fig. S8. It is obvious that at $K/J=0.24$, we have dimer order in thermodynamic limit.

VI.3 SFS phase

In $K/J>0.3$ region, we have SFS phase. We have showed real space decay of various correlations in the main text. In SFS phase, there is no spin or bond order from different structure factors and the real space decay of all correlations indicate long range correlation. In addition, we measured static spin susceptibility by applying $z$-direction magnetic fields with a coupling $-B_{z}\sum_{i}S_{i}^{z}$ and computed $M_{z}=\frac{1}{N}\sum_{i}\langle S_{i}^{z}\rangle$. As we know $\chi=\frac{\partial M_{z}}{\partial B_{z}}$, the linear fitting of the magnetic moment density with magnetic fields will give us an estimation of the static spin susceptibility. Fig S9 is such an estimation, and gives a finite value of static spin susceptibility $\chi=0.22(2)$ (we have set $J=1$ as the energy unit). Restoring units, we have $\chi=0.22(g\mu_{B})^{2}/J$. According to a band theory of half-filled free spinon on the triangular lattice, the static spin susceptibility is $(g\mu_{B})^{2}N(0)/4$, where density of states at Fermi surface is about $N(0)\approx\frac{1}{1.65\pi t}$ with $t$ the spinon hopping parameter. In combination, we can also make an estimation of the hopping strength $t$ to be $0.22J$. The definitive evidence of SFS phase comes from the observation of $2k_{F}$ surface in the spin structure factor, which we have discussed in detail in the main text.

In addition to the spin correlation function discussed in the main text, it is interesting to consider the dimer correlation function $C_{D_{b}}$ defined after Eq. (S24). In Fig. S10, we show the Fourier transform of this correlator for $b=y$ (Fig. S10(a)) and $b=xy$ (Fig. S10(b)). Note that the $D_{y}$ correlator shows a broad peak at the same M point as in Fig. S7 and looks like a broadened version of it. The $D_{xy}$ correlator should be equivalent by 60 degrees rotation but it appears more smeared in Fig. S10(b) due to finite size effect. It is thus very interesting that the SFS state retains strong remnants of the VBS phase in the dimer correlators. Perhaps the SFS phase can be thought of as the quantum melted version of the VBS phase. We also point out that if pair density wave (PDW) order is responsible for a gap on the Fermi surface along the $\Gamma$ to M direction as we suggested in the main text, the wave vector of the PDW is $2k_{F}$ in that direction. We expect to form a charge density wave or bond density wave with wave vector $4k_{F}$ in the same direction as a composite order parameter which may exhibit a slow power law decay. It turns out that up to Umklapp, this wave-vector is very close to the M point. So the peak structure in the dimer correlator shown in Fig. S10 can be regarded as a possible signature of PDW order or a Z2 spin liquid. This is only a speculative idea at this point and much more work will clearly be needed in the future.

Finally, another numerical signature of entering SFS phase is that the simulation become harder and harder to converge. In Fig. S11, we show the subsystem entanglement entropy for $L_{x}=24$, $L_{y}=6$ systems in SFS phase. The convergence is very slow even we increase bond dimension exponentially. At $m=5120$, the truncation error is already at the level of $10^{-5}$, but the entanglement entropy still increase slowly, especially at the center region of the system. Although the full system is not conformally invariant, we will proceed with the following analysis. For one-dimensional gapless state with conformally invariant correlation functions in space-time, the subsystem entanglement entropy of cylinder geometry system is Calabrese and Cardy (2004)

[TABLE]

where $l$ is the subsystem length, $N$ is the total length, $A$ is a constant and $c$ is the effective central charge which measures the number of gapless modes directly. From Fig S11, we find for small subsystem size, the entanglement is already converged, if we use this part to fit the central charge, we will get $c\approx 8.6$. The black dashed line shows the function in Eq.(S28) with this central charge value. If we do similar fitting for systems with different widths, we will find the central charge increase roughly linearly with the width of the cylinder, that is to say, the number of gapless mode increases linearly with the width of cylinder, which is consistent with the SFS phase. We note that in the U(1) SFS in a six wide system, there are five bands that cross the Fermi surface giving ten gapless modes. However, one mode should be removed due to coupling to the U(1) gauge field which enforces a constraint and we expect $c=9$. On the other hand, for a Z2 SFS, the gauge photon is gapped due to spinon pairing and if a pairing gap exists along the $\Gamma$ to M direction as we suggested in the main text, we expect four bands to cross the Fermi surface and $c=8$. Unfortunately due to the slow convergence for subsystems near the middle of the sample, we cannot determine the central charge with sufficient accuracy to distinguish between the two cases.

Although from the entanglement entropy point of view, the convergence is slow, the $2k_{F}$ peaks positions converge rather quickly, as does the spin structure factor itself when increasing bond dimensions as showed in Fig. 4 of main text.

VI.4 Phase transition

Being limited by the computation resources, we only focus on different phases here and leave the study of details of phase transitions to future works. However, from the total energy data of different $K/J$ values we have, it seems both phase transitions of AFM to VBS and VBS to SFS are continuous. Fig. S12 shows the energy density versus $K/J$ for $L_{x}=12$, $L_{y}=6$ with $\gamma=0$ systems. There is not any singularity of first derivative, although confidence of the information is limited as we need to also consider finite size broadening.

VII The role of anisotropy

The $\gamma$ in 1T-TaS2 is very close to zero and consider $\gamma=0$ is enough, but for a complete model study, we will also discuss the small anisotropy case. For a small anisotropy the overall phase diagram is the same as isotropic ($\gamma=0$) case. From the available data we have, small anisotropy may shrink the VBS phase a little bit. In Fig. S13, we show the y-direction dimer correlation structure factor at $M$ point versus $K/J$ for $L_{x}=12$, $L_{y}=6$ with $\gamma=0$ and $\gamma=0.1$. Where a suppressing of y-direction dimer order in small anisotropy case ($\gamma=0.1$) is found.

VIII Comparing with earlier studied model

Comparing our model with four-spin ring exchange model studied before Motrunich (2005); Sheng et al. (2009); Block et al. (2011), we find if we expand the permutation operators in their model explicitly to spin operators, then the relations between our model (considering $\gamma=0$ isotropic case, with free parameters $J$ and $K$) and their model (with parameters $J_{2}$ and $J_{4}$ ) will be $J=2J_{2}+5J_{4}$, $K=4J_{4}$, and an extra next nearest neighbor (nnn) two spin exchange term (with exchange strength $J_{4}$) is in their model while absent in ours. Our results will be a good reference to discuss the role of nnn term and four spin ring exchange term. On the one hand, we find SFS phase without nnn term while both VMC and DMRG simulations find SFS phase with nnn term. That means nnn term is not essential to make SFS phase, but four spin ring exchange terms should be. On the other hand nnn term may be important to make other QSL phase. For example, without four-spin ring-exchange term and considering $J_{1}$-$J_{2}$ model on the triangular lattice, several groups, with coupled-cluster Li et al. (2015c), DMRG Zhu and White (2015); Hu et al. (2015), or Gutzwiller-projected methods Kaneko et al. (2014); Iqbal et al. (2016), find QSL phase in medium $J_{2}/J_{1}$ region (for example in Ref. Zhu and White (2015); Hu et al. (2015); Iqbal et al. (2016), $0.08<J_{2}/J_{1}<0.16$), although controversies still exist in whether it is gapped or nodal QSL.