Superconductivity in the doped quantum spin liquid on the triangular lattice

Broad interest in quantum spin liquid (QSL) phases was triggered by the notion that they can be viewed as insulating phases with preexisting electron pairs, such that upon light doping they might automatically yield high temperature superconductivity. Yet despite intense experimental and numerical efforts, definitive evidence showing that doping QSLs leads to superconductivity has been lacking. We address the problem of a lightly doped QSL through a large-scale density-matrix renormalization group study of the t-J model on finite-circumference triangular cylinders with a small but nonzero concentration of doped holes. We provide direct evidences that doping QSL can naturally give rise to d-wave superconductivity. Specifically, we find power-law superconducting correlations with a Luttinger exponent, Ksc ≈ 1, which is consistent with a strongly diverging superconducting susceptibility, χsc~T−(2−Ksc)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{sc} \,\sim\, {T}^{-(2\,-\,{K}_{sc})}$$\end{document} as the temperature T → 0. The spin–spin correlations—as in the undoped QSL state—fall exponentially which suggests that the superconducting pair-pair correlations evolve smoothly from the insulating parent state.

The spin-1/2 antiferromagnet on the triangular lattice with nearest-neighbor (NN) and next-nearest-neighbor (NNN) Heisenberg interactions is highly frustrated. A number of numerical simulations have provided strong evidence that its ground state is a QSL in the region of 0.07 ≤ J 2 /J 1 ≤ 0.15, which is sandwiched by a 120°magnetic order phase and a stripe magnetic order phase [25][26][27][28][29][30][31][32] . Therefore, it can serve as an ideal platform to investigate the consequence of doping the QSL. Although it is still under debate whether this QSL is gapped or gapless in the 2D limit, it is consistent with a gapped (possibly "Z 2 " 33 ) spin liquid on finite width cylinders. This is evidenced by the observed nonzero spingap and exponentially falling spin-spin correlations where the correlation length is significantly shorter than the width (L y in Fig. 1) of the cylinders, and short-range dimer-dimer and chiral-chiral correlations [25][26][27][28][29] . Therefore, these leave little doubt that doping this state corresponds to doping a gapped QSL.
As the QSL has been realized on the spin-1/2 antiferromagnetic J 1 -J 2 Heisenberg model on the triangular lattice, a natural question is whether doping this QSL will give rise to superconductivity? Quasi-one-dimensional (1D) systems such as cylinders have become an important starting point to resolve this problem which has essential degrees of freedom that allow for two-dimensional characteristics to emerge. Moreover, it can be accurately studied using the density-matrix renormalization group (DMRG) [34][35][36] . According to the Mermin-Wagner theorem, a superconducting (SC) state that can be realized in quasi-1D systems such as cylinders has quasi-long-range SC correlation. The Luther-Emery (LE) liquid is an example of this 37 , which has one gapless charge mode with quasi-long-range SC and chargedensity-wave (CDW) correlations, but the spin-spin correlations fall exponentially [38][39][40][41][42][43][44] .
In this paper, we directly address the question of doping the QSL on the triangular lattice by studying the t-J model defined in Eq. (6). On the cylinders we study, the undoped system is a fully gapped QSL. Upon lightly doping we find that the system still shows exponentially falling spin-spin correlations. Most importantly, we find that even at the smallest doping level δ and on our largest 6-leg cylinders, the SC correlations are strong and decay with a slow power law, ΦðrÞ $ jrj ÀKsc , with K sc ≈ 1. This slow decay implies a SC susceptibility that diverges as χ sc $ T Àð2 À KscÞ as the temperature T → 0. Moreover, the SC correlations dominate over the CDW correlations in the sense that in all cases K c > K sc . This is suggestive that doping the QSL state on the triangular lattice could naturally give rise to superconductivity.

Principal results
We find that the system exhibits power-law CDW correlations corresponding to a local pattern of "partial-filled" charge stripes.
doped hole per unit cell (Fig. 2b). Similar with the 4-leg cylinders, this again corresponds to two holes per 1D unit cell.
Because L x ≫ L y , our system can be thought of as quasi-1D, and we find that the ground state of the system is consistent with that of a LE liquid 37 , which is characterized by power-law decaying CDW (Fig. 2) and SC correlations (Figs. 3, 4), but exponentially decaying spin-spin correlations (Fig. 5), with one gapless charge mode ( Fig. 6). At long distance, the spatial decay of the CDW correlation is dominated by a power-law with the Luttinger exponent K c . The exponent K c can be obtained by fitting the charge-density oscillations (Friedel oscillations) induced by the boundaries of the cylinder 45 Here nðxÞ ¼ hnðxÞi, wherenðxÞ ¼ P Ly y ¼ 1n ðx; yÞ=L y is the local rung density operator and x is the rung index of the cylinder. δn is a nonuniversal amplitude, ϕ is a phase shift, n 0 is the background density, and k F is the Fermi wavevector. Note that a few data points (Fig. 2a, b, gray color) are excluded to minimize the boundary effect and improve the fitting quality. The extracted exponent K c is similar on both width L y = 4 and L y = 6 cylinders, where we find K c > 1 at all doping levels. Similarly, K c can also be obtained from the charge density-density correlation which gives consistent results (see Supplementary Fig. 4).
At long distance, the SC pair-field correlation Φ(r), defined in Eq. (3), is also characterized by a power-law with the appropriate Luttinger exponent K sc defined by where r is the distance between two Cooper pairs along the e x direction. As shown in Fig. 3d, the extracted exponent in the long cylinder limit L x → ∞ is K sc < 1 on L y = 4 cylinders. This also holds true for L y = 6 cylinders, again we find that K sc < 1 when L x → ∞.
For both cases, this implies a divergent SC susceptibility  H.-C. Jiang χ sc $ T Àð2 À KscÞ when the temperature T → 0. As expected theoretically from the LE liquid 37 , we find that the relation K c K sc = 1 holds within the numerical uncertainty and finite-size effect (Figs. 3d, 4d, insets) for L y = 4 cylinders for doping levels δ ≤ 1/12 and L y = 6 cylinders for doping levels δ ≤ 1/18. However, a notable deviation can be seen for the doping level δ = 1/12 on L y = 6 cylinders ( Fig. 4d inset), we suspect that it might be attributed to the finite-size effect even after finite-size (L x ) extrapolation (Figs. 2d, 4d). This is also consistent with the central charge c shown in Fig. 6d, which also suffers from visible finite-size effect, although an apparent trend that c may approach to the expected value c = 1 in the long cylinder limit L x → ∞. In any case, the fact that K c > K sc (Figs. 3d, 4d) demonstrates the dominance of the SC correlations Φ(r) on both L y = 4 and L y = 6 cylinders.

Charge-density wave order
To describe the charge-density properties of the ground state of the system, we have calculated the charge-density profile n(x). Figure 2a shows the charge-density distribution n(x) for width L y = 4 cylinders, which is consistent with the "half-filled" charge stripe of wavelength λ = 1/2δ, i.e., spacings between two adjacent stripes, for all doping concentration δ with half a doped hole per unit cell. The charge-density profile n(x) for width L y = 6 cylinders is shown in Fig. 2b, which is however consistent with "third-filled" charge stripes for all δ with only one third of a doped hole per unit cell. As a result, it has a shorter wavelength λ = 1/3δ on L y = 6 cylinders than the "half-filled" stripes on L y = 4 cylinders at the same doping concentration δ.
At long distance, we find that the CDW correlations decay with a power-law with the Luttinger exponent K c . The exponent K c , which was obtained by fitting the data points using Eq. (1), is shown in Fig. 2c for width L y = 4 cylinders and Fig. 2d for width L y = 6 cylinders as a function of L x and δ. For all cases, we find that K c increases with L x , which is consistent with weaker CDW correlations in the longer cylinders. A notable difference between the L y = 4 and L y = 6 cylinders is the doping dependence of K c as shown in Figs. 3d, 4d. On the L y = 4 cylinders, the value of K c decreases with the increase of δ on width L y = 4 cylinder. On the contrary, K c increases with the increase of δ on width L y = 6 cylinders. However, in any case, our results show that K c > 1, or more precisely K c > 1.2 in the long cylinder limit L x → ∞ (Fig. 2c, d), for all doping levels on both L y = 4 and L y = 6 cylinders. K c can also be obtained directly from the charge-density-density correlation, which produces similar results (see Supplementary Fig. 4).

Superconducting correlation
To test the possibility of superconductivity, we have also calculated the equal-time SC pair-field correlations. As the ground state of the system with even number of doped holes is always found to have spin 0, we focus on spin-singlet pairing. A diagnostic of the SC order is the pair-field correlator, defined as Here Δ y α ðx; yÞ ¼ 1 ffiffi 2 p ½c y ðx;yÞ;" c y ðx;yÞ þ α;# À c y ðx;yÞ;# ; c y ðx;yÞ þ α;" is the spin-singlet pair-field creation operator, where the bond orientations are designated α = a, b, and c ( Fig. 1), (x 0 , y) is the reference bond with x 0~Lx /4, and r is the distance between two bonds in the e x direction.
Due to the presence of CDW modulations (Fig. 2), SC correlations Φ αβ (r) exhibit similar spatial oscillations with n(x). For a given cylinder of length L x , we determine the decay of SC correlations with reference bond located at the peak position around x 0~Lx /4 of the charge-density distribution n(x) to minimize the boundary effect 42,43 . This gives accurate values of Φ αβ (r) for reliable finite-size scaling. Examples of Φ aa (r) are given in Figs. 3, 4 for width L y = 4 and L y = 6 cylinders, respectively. Further details are provided in the Supplementary Fig. 2.
As the ground state of the system is a spin-singlet state, theoretically, the pair symmetry of the SC correlations will be consistent with either s-wave, d ± id or d-wave 46,47 . In the following we show that the pair symmetry is only consistent with the d-wave. First of all, our results show that the ground state wavefunction of the system is purely real and the imaginary part of all the SC correlations Φ αβ (r) is zero. As a result, the pair symmetry is inconsistent with d ± id-wave. Second, the SC correlations change sign between different bonds, i.e., Φ ab~Φac −0.4Φ aa (see Supplementary Fig. 3), so the pair symmetry is also inconsistent with s-wave. Therefore, the only option left is the dwave, which is expected to survive in two dimensions 47 . This is indeed consistent with our results, where we find that the SC correlations are stronger on a bonds but weaker on other bonds, i.e., Φ bb~Φcc~0 .2Φ aa on both the L y = 4 and L y = 6 cylinders. Taken together theory and numerical results, we therefore conclude that the pairing symmetry of the SC correlations is consistent with the d-wave.
Similar to the CDW correlations, the SC correlations Φ(r) of the system on both width L y = 4 (Fig. 3) and L y = 6 cylinders (Fig. 4) also decay with a power law, whose exponent K sc was obtained by fitting the results using Eq. (2) [42][43][44] . The fact that the observed K sc < 1 holds for all the cases when L x ≫ L y undoubtedly demonstrates the dominance of the SC correlations over the CDW correlations since K c > 1 > K sc . These findings are in stark contrast to the case of doping QSL on the Kagome lattice 22,24 , where the SC correlations decay exponentially but the CDW correlations have 2D long-range order. This indicates that the lattice geometry may play a critical role in giving rise to superconductivity in doping a QSL.

Spin-spin correlation
To describe the magnetic properties of the ground state, we calculate the spin-spin correlation functions defined as where S ! x;y is the spin operator on site i = (x, y). (x 0 , y) is the reference site with x 0~Lx /4 and r is the distance between two sites in the e x direction. It is worth noting that the oscillations have different wavelength between small r region (short distance) and large r region (long distance), and there are no antiphase domain walls in the spin-spin correlation h S ! x0;y Á S ! x0 þ r;y i. Following similar procedure as for n(x) and Φ(r), we first extrapolate F(r) to the limit ϵ → 0 for a given cylinder of length L x and then perform finitesize scaling as a function of L x . As shown in Fig. 5a, F(r) decays exponentially as FðrÞ $ e Àr=ξ s at long-distances, with a correlation length ξ s = 1 ∼ 8 lattice spacings for the L y = 4 cylinders and ξ s = 2~7 lattice spacings for the L y = 6 cylinders (Fig. 5b) for doping concentration 0 ≤ δ ≤ 1/12. Therefore, the spin-spin correlations are also short-ranged which is similar with the QSL at half-filling, where the spin-spin correlations also decay exponentially.

Central charge
A key feature of the LE liquid is that it has a single gapless charge mode with central charge c = 1, which can be obtained by calculating the von Neumann entropy S ¼ ÀTrρlnρ, where ρ is the reduced density matrix of a subsystem with length x. For critical systems in 1 + 1 dimensions described by a conformal field theory, it has been established 48,49 that for an open system of length L x , j sin kF j þS; whereÃ andS are model dependent fitting parameters, and k F is the Fermi momentum. Examples of S(x) are shown in Fig. 6a, b for the L y = 4 and L y = 6 cylinders, respectively. For a given cylinder of length L x , a few data points in blue close to the ends are excluded in the fitting to minimize the boundary effect. The extracted central charge c is shown in Fig. 6c, d as a function of L x and δ.
Although the value of c deviates notably from c = 1 on short cylinders, apparently it approaches c = 1 with the increase of L x by taking into account the numerical uncertainty and finite-size effect. The result is consistent with one gapless charge mode with c = 1, which provides additional evidence for the presence of the LE liquid in doping the QSL on the triangular lattice.

DISCUSSION
Taken together, our results show that the ground state of doping QSL on the triangular lattice is consistent with a LE liquid with a finite gap in the spin sector. Although both CDW and SC correlations decay as a power-law, the fact that K c > 1 > K sc on both the 4-leg and 6-leg cylinders when L x ≫ L y at all doping levels undoubtedly demonstrates the dominance of the SC correlations with divergent SC susceptibility. To the best of our knowledge, this is the direct numerical realization of dominant superconductivity in doping a time-reversal symmetric QSL, which provides compelling evidence that doping a QSL could naturally give rise to superconductivity [6][7][8][9] . In this paper, we have focused on the lightly doped case, it will also be interesting to study the higher doping case as well as the consequence of doping different phases on the triangular lattice such as the magnetic ordered phase 50,51 . Answering these questions may lead to better understanding of the mechanism of high temperature superconductivity.

Model Hamiltonian
We employ DMRG 34 to study the ground state properties of the holedoped t-J model on the triangular lattice, defined by the Hamiltonian Hereĉ þ iσ (ĉ iσ ) is the electron creation (annihilation) operator on site i = (x i , y i ) with spin σ. S ! i is the spin operator andn i ¼ P σĉ þ iσĉ iσ is the electron number operator. The electron hopping amplitude t ij is equal to t 1 (t 2 ) if i and j are NN (NNN) sites as shown in Fig. 1. J 1 and J 2 are the spin superexchange interactions between NN and NNN sites, respectively. The Hilbert space is constrained by the no-double occupancy condition, n i ≤ 1. At half-filling, i.e., n i = 1, Eq. (6) reduces to the spin-1/2 antiferromagnetic J 1 -J 2 Heisenberg model.

Numerical details
The lattice geometry used in our simulations is depicted in Fig. 1, where e x = (1, 0) and e y ¼ ð1=2; ffiffi ffi 3 p =2Þ denote the two basis vectors. We take the lattice geometry to be cylindrical with periodic (open) boundary condition in the e y (e x ) direction. Here, we focus on cylinders with width L y and length L x , where L x and L y are number of sites along the e x and e y directions, respectively. The total number of sites is N = L x × L y with N e = N electrons at half-filling. The doping level of the system is defined as δ = N h / N, where N h = N − N e is the number of doped holes.
For the present study, we focus on the lightly doped case with δ ≤ 1/12 on width L y = 4 cylinders of length up to L x = 144 and width L y = 6 cylinders of length up to L x = 72. We set J 1 = 1 as an energy unit and J 2 = 0.11 such that the system is deep in the QSL phase at half-filling [25][26][27][28][29][30] . We consider t 1 = 3 and t 2 ¼ t 1 ffiffiffiffiffiffiffiffiffiffi ffi J 2 =J 1 p ¼ 0:995 to make a connection to the corresponding Hubbard model 52 . The results also hold for other t 1 and t 2 . We perform up to 100 sweeps and keep up to m = 12,000 states for width L y = 4 cylinders with a typical truncation error ϵ < 3 × 10 −7 , and up to m = 30000 states for width L y = 6 cylinders with a typical truncation error ϵ < 4 × 10 −6 . This leads to excellent convergence for our results when extrapolated to m = ∞ limit. Further details of the numerical simulation are provided in the Supplementary Information sections I, II, and III.

DATA AVAILABILITY
The author declares that the main data supporting the findings of this study are available within the article and its Supplementary Information file. Extra data are available from the corresponding author upon reasonable request.