Stripe order enhanced superconductivity in the Hubbard model.

Unidirectional ("stripe") charge density wave order has now been established as a ubiquitous feature in the phase diagram of the cuprate high-temperature superconductors, where it generally competes with superconductivity. Nonetheless, on theoretical grounds it has been conjectured that stripe order (or other forms of "optimal" inhomogeneity) may play an essential positive role in the mechanism of high-temperature superconductivity. Here, we report density matrix renormalization group studies of the Hubbard model on long four- and six-leg cylinders, where the hopping matrix elements transverse to the long direction are periodically modulated-mimicking the effect of putative period 2 stripe order. We find that even modest amplitude modulations can enhance the long-distance superconducting correlations by many orders of magnitude and drive the system into a phase with a substantial spin gap and superconducting quasi-long-range order with a Luttinger exponent, [Formula: see text].

A complex relation between multiple ordering tendencies appears to be a universal feature of highly correlated electronic systems (1). For example, charge density wave (CDW), spin density wave (SDW), and d-wave superconducting (SC) orders all arise in significantly overlapping regimes of the phase diagram of the cuprate high-temperature superconductors. Moreover, a similarly delicate balance between these same ordering tendencies appears unavoidable in studies (2, 3) of the Hubbard model with a repulsive U of order the bandwidth, .

There are clear senses in which these orders “compete.” This can be seen phenomenologically in the cuprates where suppressing SC order with a magnetic field enhances the strength of the observed CDW and where the most robust SC often appears in regions of the phase diagram where the CDW order is relatively weaker (4). A similar feature is vividly apparent in density matrix renormalization group (DMRG) studies of the Hubbard model on long but relatively narrow cylinders and ladders (5–17). Here, the closest possible approximation of an SC state is a Luther–Emery liquid (18), in which the SC and CDW susceptibilities are determined by quantum mechanically dual variables. Thus, any change in the parameters (e.g., details of the band structure or the strength of the interactions) that enhances the long-distance correlations of one necessarily decreases the other. It has even been suggested that this competition is so ferocious that the Hubbard model with may never be SC in the two-dimensional limit (14).

However, the fact that high-temperature superconductivity and CDW (not to mention SDW) orders all seem to appear together suggests that they may be linked in a more multifacetted manner than the word “competing” suggests (19). Indeed, two distinct theoretical proposals carry the implication that CDW order can enhance SC. 1) It was proposed in refs. 20 and 21 that CDW fluctuations—associated with proximity to a putative CDW quantum critical point—could serve as an effective pairing “glue” and thereby, enhance SC even under conditions in which fully developed CDW order might depress SC by opening gaps on portions of the Fermi surface. 2) It was proposed in ref. 22 and further developed in a variety of subsequent papers (23–28) that static or slowly fluctuating CDW order could produce a form of “optimally inhomogeneous” electronic structure that could enhance SC.

In the present paper, we use DMRG studies of the square lattice Hubbard model on four- and six-leg cylinders with length L = 32 and 48 to explore the second of these propositions. We consider the model with only nearest-neighbor (NN) interactions t, with , and for electron density per site with and 1/12. Moreover, we assume an ordered period 2 explicit CDW with ordering vector perpendicular to the long axis of the cylinder, so that the hopping matrix elements in this direction are alternately enhanced or depressed, as shown in Fig. 1.

Fig. 1.

The Hubbard model on the square cylinder. Periodic and open boundary conditions are imposed, respectively, along the directions specified by the lattice basis vectors and . t = t and () are hopping integrals between NN sites in the and directions. U is the on-site Coulomb repulsion, and L and L are the numbers of sites.

For dt = 0, this is the uniform Hubbard model, which in this range of parameters appears (9–11, 14, 15) to favor an insulating phase with spontaneous translation symmetry breaking corresponding to an array of “full stripes” (i.e., the CDW period along the cylinder is ) (29–31). As might be expected, this state has exponentially falling SC correlations at long distances. For dt = t, this system consists of decoupled two-leg ladders. While the behavior of the two-leg ladder depends on the ratio of , as long as this ratio does not exceed a critical value (16), the two-leg ladder is known (16, 23, 32–34) to support a Luther–Emery liquid phase with power-law SC correlations that fall with distance r as with K between one and two.* Here, we explore the effect of relatively weak modulations, .

In all cases, we find that the modulation enhances the SC correlations at long distances relative to the uniform cylinder (dt = 0) by many orders of magnitude. Indeed, the modulated cylinder seemingly forms a Luther–Emery liquid. The spin–spin correlator and the single-particle Green function fall exponentially with distance with a correlation length of order a lattice constant, indicating the existence of a spin gap. Moreover, there are clear CDW correlations with wavelength for the four-leg cylinder and for the six-leg cylinder. However, while it is plausible that they also have power-law correlations characterized by Luttinger exponent K, the expected duality relation is only barely consistent with the DMRG results for the four-leg cylinder and clearly inconsistent with them for the six-leg cylinder. Thus, unambiguous identification of the conformal field theory that characterizes the long-distance properties of the six-leg cylinder is still a work in progress ().

The Model

We employ DMRG (35) to study the ground-state properties of the Hubbard model on the square lattice, which is defined by the Hamiltonian

Here, () is the electron creation (annihilation) operator on site with spin polarization σ, and is the electron number operator. We take the lattice geometry to be cylindrical with periodic (open) boundary condition in the () direction, as shown in Fig. 1. denotes NN sites. t = t, , and are the electron hopping integrals between NN sites in the and directions, respectively. Here, we focus on cylinders with width L and length L, where L and L are the numbers of sites along the and directions, respectively. The total number of sites is , the number of electrons is N, and the doping level of the system is defined as , where is the number of doped holes relative to the half-filled insulator that arises when N = N.

In the present study, we chose units of energy such that t = 1 and consider . We consider U = 12 at and doping levels and focus on L = 4- and 6-leg cylinders of length up to L = 48. We perform around 60 sweeps and keep up to states for L = 4 cylinders with a typical truncation error and up to states for L = 6 cylinders with a typical truncation error .

The results of our calculations (as explained below) are summarized for in the remaining figures and quantified in Table 1. More details, including further analysis of truncation error and results for , are provided in .

Table 1.

Summary of extracted parameters

Ly	dt	Ksc	Δd	Δs	Δ_π	Kcdw	ξs	ξG
4	0.0	1.38 (3)	0.0	0.0	0.066	1.27 (1)	8.6 (4)	3.9 (2)
4	0.1	1.22 (3)	0.019	–0.011	0.074	1.35 (1)	7.1 (2)	3.6 (2)
4	0.2	1.08 (2)	0.032	–0.016	0.082	1.46 (1)	4.7 (2)	3.0 (1)
4	0.3	1.02 (2)	0.042	–0.021	0.091	1.48 (1)	2.9 (1)	2.5 (1)
6	0.0	∞	0.0	0.0	0.0	0.3 (3)	3.9 (4)	2.4 (3)
6	0.3	1.04 (9)	0.070	0.004	0.038	3.5 (2)	1.7 (1)	1.8 (1)
6	0.4	1.03 (8)	0.062	–0.011	0.065	3.3 (2)	1.3 (1)	2.2 (1)

The parameters are obtained by fitting the DMRG results to theoretically expected asymptotic forms of various correlation functions for and the given values of L and dt. Exponentially falling correlations are represented by a Luttinger exponent of . Precise levels of uncertainty due to finite size effects—especially with regard to the Luttinger exponents—are difficult to estimate.

SC Pair-Field Correlations

We have calculated the equal-time spin-singlet SC pair-field correlation function

Here, is the spin-singlet pair creation operator on the NN bond from site (x, y) oriented in the or direction. We are interested in the decay of this quantity at large distances along the cylinder, r, as a function of both the relative orientation of the two bonds, α and β, and their relative displacement transverse to the cylinder, y. We take (x0, y0) to be the “origin,” chosen to be a site near the center of the system with and . At long distances (), exhibits power-law decay (i.e., quasi–long-range order [QLRO]) characterized by the Luttinger exponent K:

The nature of the pairing is encoded in the behavior of the amplitudes, . Specifically, were there true long-range order (i.e., in the limit ), we could classify SC states (e.g., d wave vs. s wave) by the behavior under symmetry transformations of these amplitudes. Thus, to develop some intuition concerning the meaning of these amplitudes, we analyze what they would mean in this limit. The spatial symmetries of the striped model are such that there are two inequivalent y-directed bonds and a unique x-directed bond. In a state with SC long-range order and if we assume that the translation symmetry of the model is not spontaneously broken, then the most general singlet order parameter on NN bonds can be parameterized as

In the limit dt = 0, each of these parameters would be associated with a state with different symmetries—nonzero values of Δ or Δ would characterize an “extended s-wave” or “d-wave state,” while nonzero would correspond to a period 2 pair density wave (also known as a “π-pairing” state). Note that, by symmetry, the pair field vanishes on all x-directed bonds in the π-pairing state. However, for nonzero dt, the symmetry distinction between these states is removed, so some mixture of all three is expected. However, it is still reasonable (and conventional) to refer to the case in which is the largest component as “d wave–like” pairing.

For noninfinite L, the amplitudes in Eq. can be viewed as reflecting the local symmetry of the pairing and as indicators of the preferred form of pairing that should be expected in the limit. Importantly, for dt = 0, even for noninfinite L, there is a sharp distinction between π pairing (with and ) and d wave–like pairing (with and ). To date, there is no evidence of a tendency toward π pairing on anything other than the four-leg cylinder. However, since for L = 4, π pairing is equivalent to d-wave pairing on plaquettes oriented perpendicular to the long axis of the cylinder, such a state has been seen and has been referred to in this context as “true d-wave” (9) or “plaquette d-wave” (13) pairing. More generally, for , we can loosely identify distinct states by which component is largest (dominant). (These symmetry arguments are made more precise in .)

Fig. 2 shows (i.e., between t bonds) for L = 4 cylinders at . The exponent K, obtained by fitting the results using Eq. , is for the uniform case, dt = 0.0, while for . We have also computed other components of . is shown in Fig. 2, and and are shown in . For the isotropic case with dt = 0.0, and decay exponentially as with (8, 13) and (i.e., the amplitudes are consistent with π-pairing QLRO with and ). This is consistent with previous studies of the L = 4 Hubbard and t–J models with dt = 0 (8, 10, 11, 13). The key observation is that and are significantly enhanced for dt > 0, so that they decay as a power law with a similar K as . In particular, not only is K decreased from its dt = 0 value, increases rapidly as well. For example, for dt = 0.3, , and . (More complete results are presented in Table 1.)

Fig. 2.

SC pair-field correlations. (A) and (B) on cylinders at with different dt and (C) and (D) on cylinders at with different dt on double-logarithmic scales. (C, Inset and D, Inset) and in double-logarithmic scales with dt = 0.4 on both and cylinders. r is the distance between two Cooper pairs in the direction. Note that only the central half region with is shown and used in the fitting, whereas the remaining data points from each end are removed to minimize boundary effects. The dashed lines denote power-law fitting to .

The results are still more dramatic for L = 6. Consistent with previous studies on the isotropic Hubbard model, on L = 6 cylinders with dt = 0, we find that the SC correlations are relatively weak and appear to decay exponentially with distance as shown, for , in Fig. 2 . However, as was the case for L = 4 cylinders, we find that the SC pair-field correlations are dramatically enhanced by a finite dt > 0, where we find that with . Moreover, the SC pairing symmetry is d-wave like with . For example, for dt = 0.3, , and . As summarized in , the results we have obtained for are qualitatively similar to those with . For instance, for dt = 0.3 at , and .

CDW Correlations

To measure the charge order, we define the rung density operator and its expectation value . Fig. 3 shows the charge density distribution n(x) for L = 4 cylinders, which is consistent with “half-filled charge stripes” with wavelength . This corresponds to an ordering wave vector (i.e., viewing the cylinder as a one-dimensional [1D] system, two holes per 1D unit cell). The charge density profile n(x) for L = 6 cylinders is shown in Fig. 3 , which has wavelength , consistent with “two third–filled” charge stripes. This corresponds to an ordering wave vector (i.e., four holes per 1D unit cell).

Fig. 3.

Charge density profiles. Charge density distribution n(x) at doping level on cylinders with (A) dt = 0.0 and (B) dt = 0.3 and on cylinders with (C) dt = 0.0 and (D) dt = 0.4. The exponent K is extracted using Eq. , where the red lines are fitting curves. A few data points in light gray are neglected to minimize boundary effects.

At long distance, the spatial decay of the CDW correlation is dominated by a power law with the Luttinger exponent K. The exponent K can be obtained by fitting the charge density oscillations induced by the boundaries of the cylinder (17, 33)

Here, A is an amplitude, is a phase shift, is the mean density, and for L = 4 cylinders and for L = 6 cylinders. Note that to improve the fitting quality, a few data points (corresponding to the light gray points in Fig. 3) are excluded to minimize the boundary effect. Values of K are summarized in Table 1. The fact that K > K for all cases in which dt > 0 suggests that CDW order is secondary compared with SC. The one exception is L = 6 and dt = 0, where the CDW correlations are at best slowly decaying and are clearly stronger than the SC. Our results are consistent with CDW QLRO with a value of , consistent with previous results for the t–J model (14). Note that similar values of K can also be obtained from the asymptotic falloff of the density–density correlation function, as shown in .

Spin–Spin Correlations

To describe the magnetic properties of the ground state, we calculate the spin–spin correlation functions defined as

Here, is the spin operator on site , and is the reference site with . Fig. 4 shows F(r) for both L = 4 and L = 6 cylinders at with different dt. It is clear that F(r) decays exponentially as at long distances, with a finite correlation length ξ (i.e., there must be a finite gap in the spin sector). Moreover, ξ decreases with increasing dt on both L = 4 and L = 6 cylinders. In addition, we also observe for both L = 4 and L = 6 cylinders that the spin–spin correlation has spatial modulation with a wavelength λ that is twice that of the charge (i.e., ). Values of ξ for and various values of dt are given in Table 1.

Fig. 4.

Spin–spin correlations at . (A) F(r) on cylinders with different dt and (B) F(r) on cylinders with different dt in semilogarithmic scale. Dashed lines denote exponential fit , where r is the distance between two sites in the direction.

Single-Particle Green Function

We have also calculated the single-particle Green function, defined as

Fig. 5 shows G(r) for both L = 4 and L = 6 cylinders at with different dt. The long-distance behavior of G(r) is consistent with exponential decay . The extracted correlation lengths for both L = 4 and L = 6 cylinders are comparable with ξ, as also shown in Table 1.

Fig. 5.

Single-particle Green function at . (A) G(r) on cylinders with different dt and (B) G(r) on cylinders with different dt on the semilogarithmic scale. Dashed lines denote exponential fitting , where r is the distance between two sites in the direction.

Summary of Results

What we have generically found, both for L = 4 and L = 6, over the entire investigated range of stripe modulation strength, dt, and doped hole concentration, δ, is a form of SC QLRO with exponentially falling spin and single-particle correlations and with typically weaker but presumably also power law–correlated CDW QLRO. These results are summarized in Table 1 where the values of the Luttinger exponents K and K, the various superconducting amplitudes, Δ, Δ, and , and the correlation lengths ξ and ξ are given as a function of dt for both the 4 and 6 leg cylinders.

Conclusions

It is both conceptually and practically important to understand what aspects of electronic structure are optimal for superconductivity. Circumstantial evidence has been adduced in several ways that certain organized forms of spatially inhomogeneous structure can enhance superconductivity, but we feel that the present results constitute the clearest and most unambiguous evidence to date that this is a real and robust effect. They also are interesting in the context of the still more basic question of whether the two-dimensional repulsive Hubbard model can support high-temperature superconductivity; the present results offer encouraging evidence of an affirmative answer, as they constitute some of the strongest long-range SC correlations documented to date on systems wider than four legs. It is worth acknowledging that the present results on period 2 CDW order cannot be directly compared with the situation in the cuprates, where the CDW order typically has period closer to three (Y1Ba2Cu3O) or four (Bi2Sr2CaCuO and LaSrCuO4). Nonetheless, it suggests that a more nuanced approach to the intertwining of CDW and SC orders may be appropriate in the cuprate context.

Finally, there is the question of obtaining a conceptual understanding of the numerical results we have reported. This is an ongoing endeavor. However, it is worth mentioning a possible connection between the present results and recent DMRG results that exhibit enhanced superconductivity in a lightly doped quantum spin liquid (36). Indeed, in the discussion of the “spin-gap proximity effect” in ref. 22, an analogy was made between the effects of stripe order and a mechanism based on a doped spin liquid.

It is reasonable to conclude that the low-energy magnetic fluctuations associated with antiferromagnetic order or near order (i.e., with energies small compared with the SC gap) are detrimental to SC; they would generally be expected to be pair breaking (a clear discussion is in ref. 37). However, higher-energy, short-range correlated antiferromagnetic fluctuations can produce precisely the sort of momentum-dependent interactions that are most conducive to d-wave SC. In this sense, a fully gapped spin liquid would seem to have just the right spectrum of magnetic fluctuations to be an optimal parent to a high-temperature superconductor. Indeed, it is possible to view the gap in such a state as the pairing gap of a superconductor that is waiting to be liberated. In a similar sense, the undoped (δ = 0) two-leg Hubbard ladder has a spin gap and can be viewed as a Mott insulator of preexisting Cooper pairs (rung singlets). In this sense, doping into a modulated array of effective two-leg ladders may be analogous to doping a fully gapped quantum spin liquid.