Gapless spin liquid and pair density wave of the Hubbard model on three-leg triangular cylinders

We study the ground state properties of the Hubbard model on three-leg triangular cylinders using large-scale density-matrix renormalization group simulations. At half-filling, we identify an intermediate gapless spin liquid phase between a metallic phase at weak coupling and Mott insulating dimer phase at strong interaction, which has one gapless spin mode and algebraic spin-spin correlations but exponential decay scalar chiral-chiral correlations. Upon light doping the gapless spin liquid, the system exhibits power-law charge-density-wave (CDW) correlations but short-range single-particle, spin-spin, and chiral-chiral correlations. Similar to CDW correlations, the superconducting correlations are also quasi-long-ranged but oscillate in sign as a function of distance, which is consistent with the striped pair-density wave. When further doping the gapless spin liquid phase or doping the dimer order phase, another phase takes over, which has similar CDW correlations but all other correlations decay exponentially.

Pair density wave (PDW) is a superconducting (SC) state in which Cooper pairs have finite momentum and the order parameter varies periodically in space in such a way that its spatial average vanishes. [1][2][3]. The first example of PDW is the Fulde-Ferrell-Larkin-Ovchinnikov state [4,5] which arises in a conventional s-wave superconductor in response to a small degree of spinpolarization so that the Fermi surface is spin-split. Increasing interest of PDW state has emerged as a mechanism to understand recent discoveries in underdoped cuprate superconductors, where direct observation of PDW has been made experimentally via local Cooper pair tunneling and scanning tunneling microscopy in Bi 2 Sr 2 CaCu 2 O 8+x [6][7][8] as well as the dynamical interlayer decoupling observed in 1/8 hole-doped La 2 BaCuO 4 [9,10]. Although theoretically much is known about the properties of the PDW state, its realization in microscopic models remains still very few. [11][12][13][14][15][16] These include the one-dimensional (1D) Kondo-Heisenberg model [11], extended two-leg Hubbard-Heisenberg model [12] and an extended Hubbard model with a staggered spindependent magnetic flux per plaquette on a three-leg triangular lattice. [13]. The signature of the PDW ordering was also observed in t-J model with four-spin ring exchange interaction on the four-leg triangular cylinder [14] and the t-J-like extension of the Kitaev model on the three-leg honeycomb cylinder. [17] However, there is no evidence of the PDW state found in the standard Hubbard model on systems wider than a two-leg ladder.
Aside from the QSL, a closely related question is the physics of doping it. Intuitively, QSL can be viewed as an insulating phase with preformed electron pairs such that it might produce superconductivity upon light doping. [21,[53][54][55][56][57][58] This idea was supported by recent large-scale DMRG studies that nematic d-wave, [59] and topological d ± id-wave superconductivity, [60] were observed on the lightly-doped time-reversal symmetric QSL and CSL, respectively. As for the doped Hubbard model on the triangular lattice, a number of SC states are proposed, including the d-wave, d ± id-wave, and pwave superconductivity [13,[61][62][63], however, these were arXiv:2103.07998v1 [cond-mat.str-el] 14 Mar 2021 challenged by the recent DMRG study, which reported the absence of superconductivity in doping the Hubbard model on both three-and four-leg cylinders. [64].
Principal results: In this paper, we address the above questions by studying the Hubbard model on three-leg triangular cylinders of length up to L x = 128 using large-scale DMRG simulations. Our main results are summarized in the ground state phase diagram in Fig.1. At half-filling, an intermediate time-reversal symmetric gapless spin liquid phase separates the metallic phase at weak coupling U < U c1 = 7.0 ± 0.5t and the Mott insulating dimer phase at strong coupling U > U c2 = 12.0 ± 0.5t. Distinct with the gapped CSL on four-and six-leg cylinders, [49,52] we find that the spin liquid phase on three-leg cylinders is gapless, manifest as gapless spin mode and quasi-long-range spin-spin correlations. With the chiral-chiral correlations decay exponentially at long distances, this phase preserves the time-reversal symmetry. Upon light-doping, this gapless spin liquid evolves into a state consistent with that of the striped PDW [3]. Both the SC correlations and charge-density-wave (CDW) decay as a power-law and oscillate in distance. While other correlations (singleparticle, spin-spin, and scalar chiral-chiral) are all shortrange, all these correlations are intertwined and mutually commensurate in terms of the wavevector. In contrast, a CDW phase is identified When further doping the gapless spin liquid phase with δ 10% or doping the dimer order phase.  Model and method: We employ DMRG [51] to study the ground-state properties of the Hubbard model on the triangular lattice, whose Hamiltonian is defined as iσĉ iσ is the electron number operator. t denotes the electron hopping amplitude between the nearest-neighbor (NN) sites ij , and U is the on-site Coulomb repulsion. The lattice geometry used in our simulations is depicted in the inset of Fig.1, with open (periodic) boundary condition along the e 1 (e 2 ) direction. We focus on three-leg triangular cylinders with width L y = 3 and length up to L x = 128, where L y and L x are the number of sites along the e 2 and e 1 directions, respectively. The doped hole concentration is defined as δ = N h /N , where N = 3L x is the total number of lattice sites and N h is the number of doped holes. We set t = 1 as an energy unit and consider 6t ≤ U ≤ 18t in the present study. We perform up to 69 sweeps and keep up to m = 25000 number of states with a typical truncation error ∼ 5 × 10 −7 . Further details of the numerical simulation are provided in the Supplemental Material (SM).
Gapless spin liquid: At half-filling, we identify three distinct phases (see Fig.1 and Fig.2) separated by two phase transitions at U c1 = 7.0 ± 0.5t and U c2 = 12.0 ± 0.5t. These phases are determined by various energy gaps including the single-particle gap ∆ p , charge gap ∆ c and  spin-triplet gap ∆ s defined as Here E N ↑ ,N ↓ is the ground state energy of the system with N ↑ spin-up and N ↓ spin-down electrons. Our calculations identify a metallic phase at U < U c1 where all three gaps vanish in the thermodynamic limit, consistent with previous studies [39,44,48,49]. At large U > U c2 , the ground state of the system can be mapped onto the spin-1/2 antiferromagnetic Heisenberg model. It has a dimerized ground state on three-leg cylinders [65] where all three gaps are expected to be finite in the thermodynamic limit. This is indeed consistent with our results as shown in Fig.2a-c including the dimer pattern in the inset of Fig.2c. Independently, the phase boundaries can also be determined by n d U 2 , with the double occupancy n d = 1 N i n i,↑ni,↓ [39], which exhibits peak and kink at the two phase boundaries (see Fig.2d).
We focus on the intermediate phase among these three phases. Distinct with four-and six-leg cylinders, this intermediate phase on three-leg cylinders is consistent with a gapless spin liquid, where both ∆ p and ∆ c remain finite but ∆ s vanishes in the thermodynamic limit as shown in Fig.2a-c. To further support this, we consider U = 10t as an example (deeply in the intermediate phase) and investigate the scaling behavior.
We first calculate the spin-spin correlation where S i is the S = 1/2 spin operator on site i and (x 0 , y) is the reference site with x 0 ∼ L x /4 and r is the distance between two sites in the e 1 direction. As shown in Fig.3a, it is clear that F (r) decays with a power-law at long distances which can be well fitted by F (r) ∼ r −Ks with corresponding Luttinger exponent K s = 1.1 (1). As a further test, a key feature of the gapless spin liquid is its finite gapless spin mode characterized by the central charge c. It can be obtained from fitting the von Neumann entanglement entropy, where ρ x is the reduced density matrix of a (quasi-) 1D subsystem with length x [66,67]. For critical (quasi-) 1D systems, it has been established [66,67] that S(x) = c 6 ln[ Lx π sin( πx Lx )] + const. Examples are shown in Fig.3b for cylinders of length L x = 48 and L x = 72, the extracted central charge is c = 1.0(1) suggesting that the intermediate phase has one gapless mode.
In contrast to the spin channel, a finite single-particle gap in the intermediate phase suggests that the singleparticle correlation should decay exponentially as G σ (r) ∼ e −r/ξ G with a correlation length ξ G . This is indeed the case as shown in Fig.3c, where G σ (r) decays exponentially and the extracted correlation length is ξ G = 1.1 (2). To test the possibility of time-reversal symmetry breaking, we have also calculated the scalar chiral-chiral correlation X(r), which is defined as Here χ i = S i · (S j × S k ) is the scalar chiral operator, where i, j and k label clockwise vertices of a triangle. On three-leg cylinders, we find that X(r) decays exponentially as X(r) ∼ e −r/ξχ at long distances with the correlation length ξ χ = 2.2(1). Therefore, we conclude that the intermediate gapless spin liquid phase on three-leg cylinders preserves time-reversal symmetry, in stark contrast to the gapped CSL on four-and six-leg cylinders. [49] Lightly doped gapless spin liquid: Upon light doping the gapless spin liquid, a state which is consistent with that of the striped PDW emerges where the CDW and SC pair-field correlations decay spatially in a powerlaw at long distances. We provide two detailed examples (U = 9t, δ = 1/18 and U = 10t, δ = 1/24) in Fig.4, while the conclusion holds for all parameter in the PDW+CDW phase of Fig. 1. In this paper, we have studied a sizeable system with length up to L x = 128 to suppress the finitesize effect. As shown below, the oscillation period of SC correlations is rather large, which results in the absence of the PDW correlation in previous study. [64]  Pair density wave. To test the possibility of superconductivity, we calculate the equal-time SC pair-field correlations. As the ground state with an even number of electrons always have total spin 0, we focus on spin-singlet SC correlation, which is defined as ,↓ĉ † (x,y)+α,↑ ] is spin-singlet pair creation operator living on bond α =a, b and c (see Fig.1 inset). (x 0 , y) is the reference site with x 0 ∼ L x /4 and r is the distance between two bonds in the e 1 direction. The spatial distribution of SC correlations Φ aa (r) and Φ cc (r) for the two examples are shown in Fig.4: Φ(r) exhibits clear spatial oscillation which can be well fitted by Φ(r) ∼ f (r)φ(r) for a large region of r, where f (r) sets envelope and φ(r) determines the oscillation, as discussed below. At long distances, the envelope function f (r) is consistent with a power-law decay f (r) ∼ r −Ksc . The extracted exponent is K sc = 3.6(2) for Φ aa (r) and K sc = 3.9(3) for Φ cc (r), respectively. We have also calculated the spin-triplet SC correlations, which however are much weaker than the spin-singlet SC correlations.
The spatial oscillation of the SC correlations Φ(r) is characterized by the normalized function φ(r) = (−1) r Φ(r)/f (r) as mentioned above. Examples of φ aa (r) and φ cc (r) are shown in Fig.4c-d, both of which oscillate periodically in real space and can be well fitted by φ(r) ∼ sin(Qr + θ) for φ aa (r) when r 8 and φ cc (r) when r 24. This is consistent with the striped PDW state with vanishing spatial average of φ(r). Q = 3πδ is the corresponding PDW ordering wavevector which corresponds to the wavelength λ sc = 2/3δ, i.e., λ sc = 12 for δ = 1/18 and λ sc = 16 for δ = 1/24. As we will see below, our results clearly show the relationship λ sc = λ s = 2λ c = λ χ , which is expected for the striped PDW state. Here λ s , λ c and λ χ are wavelengths of the spin-spin, CDW and scalar chiral-chiral correlations.
Charge density wave. To measure the charge order, we define the local rung density operator asn(x) = Ly y=1n (x, y) and its expectation value as n(x) = n(x) . Fig.5a shows the charge density profile n(x) on cylinders of length L x = 108 at δ = 1/18 and L x = 128 at δ = 1/24. The system forms 1/3-filled charge stripes of wavelength λ c = 1/3δ, which is the spacing between two adjacent charge stripe along the cylinder. This corresponds to an ordering wavevector K = 6πδ = 2Q with 1/3 doped hole per CDW unit cell.
At long distances, the spatial decay of CDW correlations is dominated by a power-law with the Luttinger exponent K c , which can be obtained by fitting the charge density oscillations (Friedel oscillations) induced by the boundaries of the cylinder [68,69] n(x) = n 0 + δn cos(K * x + θ)x −Kc/2 .
Here n 0 denotes the background electron density, δn and θ are model-dependent constants. Note that a few data points (Fig.5a, in gray) are excluded to minimize the boundary effect for a more reliable fit. The extracted exponent K c = 1.6(1) is shown in the inset of Fig.5a. Alternatively, K c can also be obtained from the charge density-density correlation, which gives consistent results (see SM for details).
Other correlations. To further characterize the PDW phase, we have also calculated other correlations including F (r), X(r) and G σ (r) as shown in Fig.5. Contrary to CDW and SC correlations, we find that they decay exponentially at long distances as F (r) ∼ e −r/ξs , X(r) ∼ e −r/ξχ and G σ (r) ∼ e −r/ξ G , where the corresponding correlation lengths ξ s , ξ χ and ξ G are given in Table I. It may be worth mentioning that while F (r) decays exponentially at long distances, its correlation length is fairly long ξ s ∼ 22 (2). This can be attributed to the fact that the lightly doped case is very close to the gapless spin liquid at half-filling, which has divergent correlation length. Interestingly, we find that both F (r) and X(r) exhibit clear spatial oscillation as shown in the insets of Fig.5b-c with wavelengths λ s and λ χ that are the same as that of the SC correlation, i.e., λ s = λ χ = λ sc . This gives the same ordering wavevector Q as the SC correlation. These features further support the striped PDW state in the lightly doped system.
Conclusion: In summary, we have studied the ground state properties of the Hubbard model on sizeable threeleg triangular cylinders. Based on our results, we con-clude that the exact ground state of the system has the following properties: (1) At half-filling, there is an intermediate gapless spin liquid phase which is characterized by one gapless spin mode and power-law spin-spin correlation but a gap to all charge excitations. (2) Light doping (δ 10%) the gapless spin liquid phase can give rise to a striped PDW state with power-law SC correlations with moderate exponent K sc ∼ 4 and an ordering wavevector Q. (3) There are power-law CDW correlations with an ordering wavevector K = 2Q. (4) While both spin-spin and scalar chiral-chiral correlations are short-ranged, they are mutually commensurate with both CDW and SC correlations with an ordering wavevector Q. To the best of our knowledge, this is the first numerical observation of power-law PDW correlation in the standard Hubbard model on a system wider than the 2leg ladder.

Supplemental Material
More results for lightly doped gapless spin liquid We provide more results here for the charge densitydensity correlation D(r) measured from U = 10t at δ = 1/24 in the lightly doped gapless spin liquid phase. Alternatively, the Luttinger exponent K c , which is extracted using Eq.7, can also be extracted from the densitydensity correlation D(r), defined as D(r) = 1 L y Ly y=1 | (n (x0,y) − ρ)(n (x0+r,y) − ρ) |, (S1) where ρ = 1 − δ is the electron density. For the system with quasi-long-range CDW order, D(r) should also decay with a power-law D(r) ∼ r −Kc where the exponent K c is identical within the error bar to the one extracted from Eq.7 in the thermodynamic limit. Fig.S1 shows the extracted K c from both methods L x = 128 cylinder. In Fig.S1a, the extracted K c ( ) from the Friedel oscillation as a function of the finite truncation error is K c = 1.6(1). For comparison, Fig.S1b shows the density-density correlation D(r) of the same cylinder, where the extracted exponent is also K c = 1.6(1), which is as expected consistent with that obtained by the Friedel oscillation the within the error bar. We can estimate the central charge c for (quasi-) 1D system of length L x using a more precise formula [66,67], Here A,c and q are model dependent constants. As shown in Fig.S2, the estimated central charge is c = 1.3(1) for U = 9t at δ = 1/18, and 1.2(1) for U = 10t at δ = 1/24, respectively. This is fairly close to c = 1 which suggests that there is probably single gapless charge mode.
Lightly doped Dimer phase To characterize the ground-state properties of lightly doped dimer phase, we consider U = 18t as an example which is deep inside the dimer phase and consider a typical doping concentration δ = 1/18. As shown in Fig.S3a, the charge density distribution n(x) has a well-defined ordering wavevector K = 6πδ in the e 1 direction and can be fitted by the Freidel oscillation (Eq.7) with K c = 1.5(1) in the = 0 limit. Alternatively, K c can also be obtained from the algebraic fitting of D(r) ∼ r −Kc , as shown in Fig.S3b, with the extracted exponents K c = 1.5(1), which is also consistent with that obtained from the Freidel oscillation.
We have also calculated various other correlation functions as shown in Fig.S4 and find that they all decay exponentially at long distances. These include the singleparticle correlation G σ (r) ∼ e −r/ξ G , the superconducting correlation Φ αβ (r) ∼ e −r/ξsc , the spin-spin correlation F (r) ∼ e −r/ξs , and the scalar chiral-chiral correlation X(r) ∼ e −r/ξχ . The corresponding correlation lengths are summarized in the main text in Table.I. Aside from the correlation functions, we have also calculated the central charge c. For the critical (quasi-) 1D system, the von Neumann entanglement entropy of the subsystem follows the formula S(x) = c 6 ln[ Lx π sin( πx Lx )] + const, with the central charge c equals to 1. As shown in Fig.S4d, the central charges extracted is c = 1.04 (7). This is consistent with single gapless charge mode.
Compare results with complex DMRG code Previous study [64] suggests lightly doping the intermediate phase could lead to a chiral metal phase, which spontaneously breaks the time-reversal symmetry. As a result, an important numerical check which needs to be done is whether real-value and complex-value DMRG simulations provide the qualitatively same results. In this section, we provide direct evidences that both DMRG simulations indeed gives us the similar results, both qualitatively and quantitatively. As examples, we have calcu- lated both the spin-spin F (r) and scalar chiral-chiral correlations X(r) on L x = 36 ∼ 108 systems by keeping up to m = 25000 number of states in the real-value DMRG simulation and up to m = 16000 number of states in the complex-value DMRG calculation. As shown in Fig.S5, while the results obtained from complex-value DMRG simulations suffer from a significantly larger boundary effects for relatively small systems (such as N = 36 × 3), they are clearly consistent with the results obtained from real-value DMRG simulations, including both F (r) and X(r). For both DMRG simulations, it is clear that F (r) and X(r) decay exponentially and can be wellfitted by an exponential function as F (r) ∼ e −r/ξs and X(r) ∼ e −r/ξχ . Moreover, we have checked and calculated other correlations, including single particle and SC correlations, where both real-and complex-value DMRG simulations give us the same results.