Pair-Density-Wave in the Strong Coupling Limit of the Holstein-Hubbard model

A pair-density-wave (PDW) is a novel superconducting state with an oscillating order parameter. A microscopic mechanism that can give rise to it has been long sought but has not yet been established by any controlled calculation. Here we report a density-matrix renormalization group (DMRG) study of an effective $t$-$J$-$V$ model, which is equivalent to the Holstein-Hubbard model in a strong-coupling limit, on long two-, four- and six-leg triangular cylinders. While a state with long-range PDW order is precluded in one dimension, we find strong quasi-long-range PDW order with a divergent PDW susceptibility as well as the spontaneous breaking of time-reversal and inversion symmetries. Despite the strong interactions, the underlying Fermi surfaces and electron pockets around the $K$ and $K^\prime$ points in the Brillouin zone can be identified. We conclude that the state is valley-polarized and that the PDW arises from intra-pocket pairing with an incommensurate center of mass momentum. In the two-leg case, the exponential decay of spin correlations and the measured central charge $c\approx 1$ are consistent with an unusual realization of a Luther-Emery liquid.

A pair-density-wave (PDW) is a novel superconducting state with an oscillating order parameter. A microscopic mechanism that can give rise to it has been long sought but has not yet been established by any controlled calculation. Here we report a density-matrix renormalization group (DMRG) study of an effective t-J-V model, which is equivalent to the Holstein-Hubbard model in a strong-coupling limit, on long two-, four-and six-leg triangular cylinders. While a state with long-range PDW order is precluded in one dimension, we find strong quasi-long-range PDW order with a divergent PDW susceptibility as well as spontaneous breaking of time-reversal and inversion symmetries. Despite the strong interactions, the underlying Fermi surfaces and electron pockets around the K and K points in the BZ can be identified. We conclude that the state is valleypolarized and that the PDW arises from intra-pocket pairing with an incommensurate center of mass momentum. In the two-leg case, the exponential decay of spin correlations and the measured central charge c ≈ 1 are consistent with an unusual realization of a Luther-Emery liquid.
In a PDW state, the Cooper pairs have non-zero center of mass momentum [1]. Long-studied versions of such states are the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phases [2,3] that can arise in the context of BCS theory in the presence of a small degree of spin magnetization. However, it has been conjectured that a PDW could also arise from a strong-coupling mechanism and could be responsible for the dynamical layer decoupling phenomena observed in certain cuprate materials [4][5][6]. While PDW states have been shown to be energetically competitive in certain mean-field and variational calculations [7][8][9][10], and evidence has been sought in numerical studies of generalized t-J [4,[11][12][13][14] and Hubbard models [15,16], PDW long-range order has not been established in any controlled calculation in two or higher dimensions. In one-dimensional systems -which can be reliably treated using DMRG -the closest one can get to such order is PDW quasi-long-range order with a divergent susceptibility. However, even this has only been seen to date for the case of commensurability-two PDW order in Heisenberg-Kondo chains [17,18]; other reported DMRG sightings [14,19,20] have found PDW correlations that fall off fast enough to yield only finite superconducting susceptibility.
Previously [21], we analyzed a strong coupling limit of the Holstein-Hubbard model, where the on-site direct electron-electron repulsion, U e-e , is stronger than the phonon induced attraction, U e-ph . We showed that the low energy physics is captured by an effective model of a form similar to the t-J-V model derived from the strong coupling limit of the ordinary Hubbard model, but in an unfamiliar range of parameters. In the effective model, the electron hopping is exponentially suppressed by a Frank-Condon factor, e − U e-ph 2ω D , relative to its bare value, * Electronic address: kivelson@stanford.edu † Electronic address: yaohong@tsinghua.edu.cn t, where ω D is the optical phonon frequency. Meanwhile, the nearest-neighbor magnetic exchange coupling J and density repulsion V , are relatively unsuppressed. Thus one can readily access novel regimes in which the effective interactions are comparable to or larger than the electron kinetic energy. Specifically, we analyzed the model on the two-dimensional triangular lattice in the adiabatic limit, where the suppression is large. We found a PDW when the electrons are dilute and the bare hopping integral t is negative (i.e. the band maximum occurs at k = 0 while minima occur at the K and K points). Unfortunately, in the adiabatic limit, the PDW is a Bose-Einstein condensate (BEC) of real space pairs, in which the same Frank-Condon factor suppresses the condensation temperature, meaning that this order would be destroyed at exponentially small temperatures.
In the present work, we employ DMRG to study the effective model in a range of parameters suggested by the above discussed analysis, but now for the case of a finite Frank-Condon factor [28]. For all the cylinders we can treat in this way, we find a divergent PDW susceptibility, as shown in Fig. 2. The PDW state spontaneously breaks time-reversal and inversion symmetries, resulting in a local pattern of equilibrium currents as shown in the inset of Fig. 3(a). The state is far from the BEC limit as can be seen in Fig. 3(b) from the existence of a sharp drop in the electron occupation at well-defined Fermi points corresponding to electron pockets around the K and K points in the Brillouin Zone (BZ). However, the size of the pockets (i.e. the value of 2k F ) is a factor of two larger than what it would be for non-interacting electrons, which we will argue reflects spontaneous valley polarization. As shown in Figs. 2(a)&3(c), the PDW ordering wavevector, Q, is incommensurate and density dependent. We will show that the ordering vector takes the value expected for intra-pocket pairing.
Model and Methods. -The effective Hamiltonian (with implicit projection to states with no doubly occupied sites) is: where t 1 is the renormalized nearest-neighbor hopping, t 2 is a weak next-nearest-neighbor hopping term via an intermediate site m, i, m, j represents a triplet of sites such that m is a nearest-neighbor of two distinct sites i and j, and (τ + 2t 2 ) is a weak singlet hopping term whereŝ ij = (ĉ i,↑ĉj,↓ +ĉ j,↑ĉi,↓ )/ √ 2 is the annihilation operator of a singlet Cooper pair on bond ij . Explicit expressions are given in Ref. [21] for the values of these effective couplings as functions of the parameters in the original Holstein-Hubbard model. For the data shown in this study, we fix t = −1, ω D = 3, U e-e = 22, and U e-ph = 18, which leads to J ≈ 0.208, V /J ≈ 0.666, t 1 /J ≈ −0.240, t 2 /J ≈ 0.0329 and τ /J ≈ 0.0445 in the effective model. However, we have checked that in a broad range of parameters, or setting t 2 = τ = 0, the results we find do not change qualitatively. We include some results with different settings in the Supplemental Material.
We perform DMRG on triangular cylinders with N x × N y sites with open boundary conditions in the x-direction and periodic boundary conditions in the y-direction. The x-direction is aligned with a primitive vector, as shown in Fig. 1 2 ) are the three primitive vectors of the triangular lattice. The compactification of the lattice to a cylinder restricts us to N y even. The momentum values in the y-direction take values k y = m Ny 4π √ 3 , m ∈ Z, indicated in Fig. 1(b). We will primarily focus on the two-leg case but some key findings are also verified in the four-and six-leg cases. Note that two-leg cylinders can be flattened to be a purely onedimensional chain by simply neglecting the y coordinate of each site, so we plot the data for all sites in the two-leg case. All the DMRG data collected are obtained from the lowest energy state out of five trials with independently randomized initial states and all the results shown (unless otherwise stated) are extrapolated to zero truncation error, utilizing data collected with six truncation errors ranging from 3×10 −6 to 1×10 −7 . In the two-leg case, we have checked our results do not change significantly down to truncation error 7 × 10 −11 , corresponding to keeping bond dimensions up to m = 4320. The correlation functions shown are averaged over all legs of the ladder, taking N 0 = 5N y different reference sites r 0 centered around the middle of the system. All data involving sites within N x /4 to the open boundary are discarded, i.e. we only retain the data on the interval x ∈ [N x /4, 3N x /4], to reduce boundary effects. Thus, each structure factor we compute utilizes the data on N = N x N y /2 sites.
Pair density wave correlations. -We first examine the equal-time pair-pair correlations. The singlet pairing order parameter is defined on the nearestneighboring bonds to a given site r i as ∆ a ( r i ) ≡ c ri,↑ĉ ri+êa,↓ +ĉ ri+êa,↑ĉ ri,↓ . The pair-pair correlator is de- . Our principle finding is that P ab ( r) oscillates in sign and falls slowly in magnitude at large distance, i.e.
This can be seen in Fig. 2(a), which is a plot of P ab ( r)/A ab |x| −νSC at large distance for a two-leg cylinder, taking ν SC = 1.20. Without breaking of the reflection symmetry across the x direction,ê 2 andê 3 are equivalent, so we only show pairing correlations involving ∆ 1 and ∆ 2 . The fact that all the components of P ab oscillate in-phase with the same wavevector, Q, indicates the local s-wave character of the PDW. An illustration of the pair correlation in the four-leg case can be seen in the inset of Fig. 2(a), which matches our previous prediction in Ref. [21] at short distance. Correspondingly, the Fourier transform of the correlator, i.e. the structure factor S PDW has pronounced peaks at two non-zero momenta and vanishes at zero momentum, as shown in Fig. 2(b). Around each peak but outside a cutoff window of width δk x ∼ 2π Nx (to avoid finite size effects), the structure factor is well-fitted by a functional form A − B|k x − Q| νSC−1 as is expected, given the long-distance correlation behavior in Eq. (2). It is important to note that the wavevector of the PDW, |Q|, appears to be a smooth monotonic function of the electron density and is not locked to a multiple of the momenta at the K points, i.e. the PDW is incommensurate with the lattice.
Since 1D systems generically exhibit emergent Lorenz invariance, we can infer from the equal time correlator that the static susceptibility to PDW order should vary as as the system size tends to infinity and temperature to zero [22]. Thus, χ(Q) diverges at zero temperature for ν SC < 2, as we find is the case for a wide range of electron densities and model parameters. That the divergence of susceptibility is not restricted to the two-leg case can be seen in Fig. 2(c), which shows the structure factor and the extracted exponents at the same electron density n = 0.05 for the four and six-leg cases.
Spontaneously broken time-reversal and inversion symmetries. -Another prominent feature of the ground state is that it spontaneously breaks both time-reversal and inversion symmetries. This can be directly seen from the current-current correlation σ (ĉ † ri,σĉ ri+êa,σ −ĉ † ri+êa,σĉ ri,σ ) is the current on the bond directed in theê a direction. (The actual current operator that we compute also receives minor contributions from t 2 and τ , see Supplemental Materials for the full expression). As shown in Fig. 3(a), the current-current correlation oscillates around a non-zero value at long distance, and we also confirm in Supplemental Materials that the peak of the structure factor of the current-current correlation scales linearly with system size. These facts signify persisting currents in the ground states. The pattern of the current flows is shown in the inset, indicating that the ground state is an orbital anti-ferromagnet.
-We further identify the broken symmetry states by investigating the occupation number in momentum space, n( k) ≡ 1 N ri, rj ,σ e i k·( ri− rj ) c † ri,σ c rj ,σ , for the two-leg cylinder. As shown in Fig. 3(b), although the strong interactions shift a small fraction of occupation weight to the vicinity of zero momentum (far above the non-interacting Fermi surface), most of the weight is confined to narrow intervals of k about the two minima of the non-interacting bands, which occur at K = (4π/3, 0) and K = −K. However, the width of these intervals -labeled "2k F " in the figure -is twice as large as it would be for a noninteracting system. Rather, this observation is consistent with the supposition that the ground-state is valley polarized [29]. We thus conclude that the presence of electrons at both K and K is an artifact of the fact that, for a finite size system, the ground-state is a superposition of states with the two possible senses of valley polarization. This is also consistent with the observation that the maximum value of n k is somewhat less than 1, while for an unpolarized non-interacting system it should equal 2 for all states below the Fermi energy. The occupancy of only one valley naturally explains the observed broken symmetries and the loop currents order, since the persisting current can be estimated as J a ≈ ±2t 1 n sin( K ·ê a ), where ± corresponds to K or K valley is occupied. Note that the current pattern is translationally invariant, and correspondingly the pattern of flux breaks point group symmetry but not translation symmetry -it is not induced by the PDW order.
The valley polarization can be understood with a simple mean field theory. We propose a trial Hamiltonian to capture the essence of the physics. To simplify the problem, we neglect the constraint on no double occupancy. Further neglecting the weak τ and t 2 terms, we solve the mean field equationt i,a,σ ≈ −t 1 − J 2 ĉ † ri+êa,σĉ ri,σ − V ĉ † ri+êa,σĉ ri,σ , to the leading order in n. We find that the valley-polarized solution, with complex hopping elementst i,a,σ = t 1 e iθ and θ ≈ ± √ 3 4 J/2+V t1 n, is always energetically favored. The band structure is renormalized to ( k)/|t 1 | = 2 cos(k x + θ) + 4 cos The introduction of the complex phase θ thus energetically distinguishes the two valleys by 9 2 (J/2 + V )n, the amount of which is sufficiently large to fully valleypolarize the system while keeping the positions of the band minima. This mechanism is similar to that of Stoner magnetism, but here the density of states is divergent at the band bottom so the polarization always occurs in the dilute limit. (Conversely, a finite critical interaction strength is necessary in dimensions d ≥ 2.) Therefore, the pairing in each ground state with valley polarization must happen between the two Fermi points located at k x =K ± k F and the pair momentum should be twice the center momentum 2K. This intra-valley singlet pairing mechanism is distinct from the intra-valley triplet or the inter-valley pairing mechanism proposed for a spin-valley locked system [16,23], and is enforced by the broken symmetries. This weak-coupling, mean-field pairing mechanism is complementary to the strong coupling, BEC-type mechanism proposed in Ref. [21]. Furthermore, due to the asymmetry of the non-interacting band structure around its band minima,K can be calculated to the leading order in n: In Fig. 3(c), we see that the observed Q indeed matches 2K calculated in this way modulo a reciprocal lattice vector, G = (4π, 0), up to an error of order π Nx . This accounts for the incommensurate nature of the PDW.
Luther-Emery liquid. -All the long-distance behaviors we have observed are consistent with a Luther-Emery liquid [24], at least in the two-leg case. For instance, as shown in Fig. 4(a), for electron density n = 0.15, the oscillatory piece of the charge correlator C( r) ≡ 1 N0 r0 n r+ r0 n r0 has wave-vector 2k F = 2nπ and exhibits a power-law decay with exponent ν CDW = 0.90, such that the expected relation ν CDW · ν SC = 1.08 ≈ 1 is approximately satisfied. The spin correlator S( r) ≡ 1 N0 r0 S r+ r0 · S r0 is short-ranged corresponding to a correlation length ξ ≈ 4.36, as shown in Fig. 4(b). To extract the central charge, we computed the von Neumann entanglement entropy S E (x) ≡ −tr(ρ x ln ρ x ) where ρ x is the reduced density matrix of the subsystem on one side of a cut at x. For critical systems in 1 + 1 dimensions described by conformal field theory, it has been established [25,26] that for an open boundary system with length N x , where c is the central charge, A and B are modeldependent parameters, and q is an adjustable fitting parameter that should approach the Fermi momentum k F in the thermodynamic limit. We see in Fig. 4(c) that this formula fits well with the observed data. As shown in the inset, the central charge c ≈ 1 and q → k F = nπ when N x → ∞, which confirms that there is only one gapless mode resulting from filling only one of the two valleys. Extensions and speculations: -With increasing number of legs, the PDW correlations seem to weaken somewhat, as suggested by the growth of the exponent ν SC . We speculate that this is an artifact in the few-leg system: at low electron densities, the Fermi pockets are small, so the bands with non-zero k y and higher energies remain unoccupied for a range of n. The effect of adding legs only leads to an increase of the filling of the k y = 0 band. In other words, the true scaling process to the two-dimensional limit has not started until the allowed values of k y become sufficiently closely spaced that additional bands with non-zero k y cross the Fermi surface.
To have more than one band crossing the Fermi surface, the required number of legs scales as 4π/ √ 3 n , which is ∼ 12 for n = 0.05. Going to larger n could in principle help improve the situation, but this will generally place the system outside of the range of n the PDW occursat least for ladders with more than 2 legs. We have not found a range of parameters amenable to DMRG in which the PDW exists and there are multiple Fermi points. Indeed, the results obtained here in the few-leg cases are essentially still one-dimensional; the main insight we obtain from them that could help understanding the 2D case is that valley polarization physics accounts well for our DMRG results in the 1D case, and this physics should apply as well in two dimensions.
More generally, our findings suggest a new route to PDW order in two or higher dimensions -one that does not rest on complicated strong-coupling physics. Consider the situation in which, either due to explicit or spontaneous time-reversal symmetry breaking, a low density of weakly interacting quasiparticles form a Fermi pocket about a band extremum at a single point in the BZ (such as the K point in the present context) about which the dispersion is not symmetric. Then from this starting point, it is natural to consider a pairing instability involving electrons on opposite sides of this pocket. To the extent that the dispersion relation can be treated as quadratic (i.e. in the effective mass approximation), these states would be perfectly nested, so deviations from nesting (e.g. trigonal warping) are small in proportion to a power of the electron density. Consequently, the pairing instability occurs for correspondingly weak couplings, and the resulting state is a PDW with ordering vector or vectors that are likewise continuously varying functions of n, as in Eq. 5. This is, in a sense, an orbital version of the original FFLO mechanism.
Gordon and Betty Moore Foundations EPiQS through Grant No. GBMF4302 at Stanford.
Competing interest statement. Author S.A.K is Editor-in-Chief of npj Quantum Materials. The authors declare no other competing interests.
Author Contribution. KSH and ZH contributed equally to this work. KSH performed the numerical investigation with the help from ZH. All authors analysed the data and contributed to the writing of the paper.

Data availability.
The authors declare that the data supporting the findings of this study are available within the paper and its supplementary information files.

Code availability.
The codes implementing the calculations of this study are available from the corresponding author upon request. B. The key observables without t2 and τ terms In Supplementary Figure 2, we show some key data for the same parameters as in main text but setting t 2 = τ = 0, i.e. the pure t-J-V model. We confirm that, without the t 2 and τ terms, the PDW correlation weakens in comparison to the original effective model, but the exponent ν SC = 1.90 < 2 still yields a divergent susceptibility, as shown in Supplementary Figure 2(a). The occupation number and thus the valley polarization picture is basically unaffected, except that there is no longer a small fraction of occupation weight being shifted to the vicinity of zero momentum, as shown in Supplementary Figure 2  The current operator can be derived from the time derivative of the density operator: where the last two terms will contribute to the current on both bond im and mj . Therefore, the current operator on bond ab can be written as: The z component of the spin current operator on bond ab can be written asĴ s b→a =Ĵ ↑ b→a −Ĵ ↓ b→a with: The singlet hopping term does not contribute to the spin current. In Supplementary Figure 3(a) we plot the magnitude of the zero-momentum current-current structure factor at various lengths of the two-leg cylinder. We find that its magnitude scales linearly with the system size, signifying long-range order and thus the existence of persisting currents in the ground states. In Supplementary Figure 3(b), we confirm the exponential decay of the spin current correlations at long distances, which rules out the possibility of spin-valley locked polarization.

D. Fermi density for four-leg case
To verify that the evidence of valley polarization still shows up for more-leg case, we plot the Fermi density for fourleg case in Supplementary Figure 4. We see that the Fermi surface size is still twice what would be for non-interacting system, and the Fermi density inside the pocket is slightly less than one.

E. Truncation error extrapolation
In Supplementary Figure 5, we show the truncation error extrapolation data for some main results shown in the main text for the two-leg case. Each data point is extrapolated to zero truncation error by a quadratic fit (marked as = 0), utilizing data collected with six truncation errors ranging from 3 × 10 −6 to 1 × 10 −7 . We have also checked that the extrapolated results do not differ significantly from the result obtained with a very small truncation error. For more leg cases, we did the same extrapolation for all results shown in the main text. For those "sanity check" simulations with very small truncation error, the maximum bond dimension is 3125 and the simulation took 126 cpu hours for 4-leg system N x = 100, n = 0.05 and truncation error 7 × 10 −11 ; the maximum bond dimension is 2602 and the simulation took 260 cpu hours for 6-leg system N x = 60, n = 0.05 and truncation error 7 × 10 −10 . We remark that, although the bond dimensions look relatively small (due to the smallness of Hilbert space in dilute limit), all the simulations are still quite time-consuming, due to the complicated form of the Hamiltonian.