Research Article

Stripe order in the underdoped region of the two-dimensional Hubbard model

See allHide authors and affiliations

Science  01 Dec 2017:
Vol. 358, Issue 6367, pp. 1155-1160
DOI: 10.1126/science.aam7127

Numerics converging on stripes

The Hubbard model (HM) describes the behavior of interacting particles on a lattice where the particles can hop from one lattice site to the next. Although it appears simple, solving the HM when the interactions are repulsive, the particles are fermions, and the temperature is low—all of which applies in the case of correlated electron systems—is computationally challenging. Two groups have tackled this important problem. Huang et al. studied a three-band version of the HM at finite temperature, whereas Zheng et al. used five complementary numerical methods that kept each other in check to discern the ground state of the HM. Both groups found evidence for stripes, or one-dimensional charge and/or spin density modulations.

Science, this issue p. 1161, p. 1155


Competing inhomogeneous orders are a central feature of correlated electron materials, including the high-temperature superconductors. The two-dimensional Hubbard model serves as the canonical microscopic physical model for such systems. Multiple orders have been proposed in the underdoped part of the phase diagram, which corresponds to a regime of maximum numerical difficulty. By combining the latest numerical methods in exhaustive simulations, we uncover the ordering in the underdoped ground state. We find a stripe order that has a highly compressible wavelength on an energy scale of a few kelvin, with wavelength fluctuations coupled to pairing order. The favored filled stripe order is different from that seen in real materials. Our results demonstrate the power of modern numerical methods to solve microscopic models, even in challenging settings.

Competing inhomogeneous orders are a common feature in many strongly correlated materials (1). A famous example is found in the underdoped region of the phase diagram of the high-temperature cuprate superconductors (HTSCs). Here, multiple probes—including neutron scattering, scanning tunneling microscopy, resonant x-ray scattering, and nuclear magnetic resonance spectroscopy—all lend support to various proposed inhomogeneous orders, such as charge, spin, and pair density waves, with suggested patterns ranging from unidirectional stripes to checkerboards (2, 3). Recent experiments on cuprates indicate that the observed inhomogeneous orders are distinct from, and compete with, pseudogap physics (4, 5).

Much theoretical effort has been directed toward explaining the origin of the inhomogeneities (6). Numerical calculations on microscopic lattice models have provided illuminating examples of the possible orders. The prototypical lattice model used to understand HTSCs is the two-dimensional (2D) Hubbard model on a square lattice, with the HamiltonianEmbedded Imagewhere a (a) denotes the usual fermion creation (annihilation) operators; n is the number operator; t and U are the kinetic and repulsion energies, respectively; and i and j are lattice site indices, where 〈ij〉 indicates that the summation is over nearest neighbors. A large number of numerical techniques have been applied to compute the low-temperature and ground-state phase diagram of this model. Early evidence for unidirectional stripe ordering in the Hubbard model came from Hartree-Fock calculations (710), whereas the nonconvex energy versus filling curves in exact diagonalization of small clusters of the t-J model (derived from the Hubbard model at large U where double occupancy is eliminated) were interpreted as signs of macroscopic phase separation (11, 12). Since then, inhomogeneous orders have been obtained, both in the Hubbard and t-J models, from calculations using the density matrix renormalization group (DMRG) (1315), variational quantum Monte Carlo (16) and constrained-path (CP) auxiliary-field quantum Monte Carlo (AFQMC) (17), infinite projected entangled pair states (iPEPS) (18), density matrix embedding theory (DMET) (19), and functional renormalization group (20), among others, although the type of inhomogeneity can vary depending on the model and numerical method. However, there are other sophisticated simulations—for example, with variational and projector quantum Monte Carlo (21, 22) and cluster dynamical mean-field theory—which do not see, or are unable to resolve, the inhomogeneous order (23, 24). The most recent studies with iPEPS (18) and DMET (19), as well as some earlier variational calculations (16, 2527), further show that both homogeneous and inhomogeneous states can be stabilized within the same numerical methodology, with a small energy difference between homogeneous and inhomogeneous states on the order of ~0.01t per site.

The small energy differences between orders mean that very small biases in ground-state simulations, such as those from an incomplete treatment of fluctuations, using insufficiently accurate constraints to control the sign problem, or from finite-size effects, can easily stabilize one order over the other. Similarly, the low temperatures needed to resolve between orders is a challenge for finite temperature methods (28, 29). Settling the resulting debate between candidate states has thus far been beyond reach. In this work, we demonstrate that, with the latest numerical techniques, obtaining a definitive characterization of the ground-state order in the underdoped region of the 2D Hubbard model is now an achievable goal. As a representative point in the phase diagram, we chose the well-known 1/8 doping point at strong coupling (U/t = 8). Experimentally, this doping corresponds to a region of maximal inhomogeneity in many HTSCs, and, in the strong coupling regime, it is recognized as a point of maximum numerical difficulty and uncertainty in simulations (24). Using state-of-the-art computations with detailed cross-checks and validation—including newer methodologies such as iPEPS and DMET, as well as recent developments in established methodologies such as CP-AFQMC and DMRG—and with exhaustive accounting for finite size effects combined with calculations directly in the thermodynamic limit, we are able to finally answer the question, What is the order and physics found in the underdoped ground state of the 2D Hubbard model?

Computational strategy

An important strategy we use to address this part of the Hubbard-model phase diagram is to combine the insights of multiple numerical tools with complementary strengths and weaknesses. This approach, pioneered in (24), greatly increases the confidence of the numerical characterization. To understand what each method contributes, we briefly summarize the theoretical background and corresponding sources of error. Further details are provided in (30).

Auxiliary-field quantum Monte Carlo

AFQMC expresses the ground state of a finite system through imaginary time evolution, Embedded Image, where Embedded Image is an initial state. The projection is Trotterized, and the evolution reduces to a stochastic single-particle evolution in the presence of auxiliary fields generated by the Hubbard-Stratonovich decoupling of the Hubbard repulsion. Away from half-filling, this decoupling has a sign problem. We use the CP approximation to eliminate the sign problem at the cost of a bias dependent on the quality of the trial state (31, 32). In this work, the Trotter error is well converged, and we report the statistical error bar. To minimize the CP bias, we use several different trial states, including a self-consistently optimized trial state (33). The calculations are carried out on finite cylinders with open, periodic, and twist-averaged boundary conditions, with widths of up to 12 sites and lengths of up to 72 sites. This method can reach large sizes, and finite-size effects are minimized. The uncontrolled error is from the CP approximation.

Density matrix renormalization group

DMRG is a variational wave function approximation using matrix product states (MPS), which are low-entanglement states with a 1D entanglement structure. The quality of the approximation is determined by the bond dimension (matrix dimension) of the MPS. The calculations are carried out on finite cylinders with widths of up to 7 sites and lengths of up to 64 sites, with periodic boundary conditions in the short direction and open boundaries in the long direction. Two different DMRG algorithms were used: one formulated in a pure (real-space) lattice basis and the other in a mixed momentum and lattice (hybrid) basis, with the momentum representation along the short periodic direction (34). We remove the bond-dimension error and finite size error in the long direction by well-known extrapolation procedures and report the associated error bar (35). Consistency between the lattice and hybrid DMRG algorithms provides a strong validation of this error bar. The remaining uncontrolled error is the finite-width error in the periodic direction.

Density matrix embedding

DMET is a quantum embedding method that works directly at the thermodynamic limit, although interactions are only accurately treated within an impurity cluster (36). To solve the impurity problem, which consists of a supercell of the original lattice coupled to a set of auxiliary bath sites, we use a DMRG solver. We treat supercells with up to 18 sites. The error bar reported in DMET corresponds to the estimated error from incomplete self-consistency of the impurity problem. The remaining uncontrolled error is the finite impurity-size error.

Infinite projected entangled pair states

iPEPS is a variational approach that uses a low-entanglement tensor network ansatz natural to 2D systems (3739). The calculations are carried out directly in the thermodynamic limit, where different supercell sizes, including up to 16 sites, are used to stabilize different low-energy states (with different orders commensurate with the supercell). As in DMRG, the accuracy of the ansatz is systematically controlled by the bond dimension D of the tensors. Estimates of quantities in the exact D limit are obtained by using an empirical extrapolation technique, which is a potential source of uncontrolled error.

Cross-checks: Systematic errors and finite-size biases

The use of multiple techniques allows us to estimate the uncontrolled errors from one technique using information from another. For example, by carrying out simulations on the same finite clusters in the AFQMC and DMRG calculations, we can estimate the CP bias in AFQMC. Similarly, in the AFQMC calculations, we can treat larger-width cylinders than is possible in the DMRG simulations; thus, we can estimate the finite-width error in DMRG.

In all of the methods, there is a bias toward orders commensurate with the shape of the simulation cell, be it the finite lattice and boundary conditions in AFQMC and DMRG, the impurity cluster in DMET, or the supercell in iPEPS. Using this bias, together with different boundary conditions and pinning fields, we can stabilize different metastable orders. For example, by setting up clusters commensurate with multiple inhomogeneous orders and observing the order that survives, we can determine the relative energetics of the candidate states. We can fit the orders along the short or the long axis of the cluster to obtain two independent estimates of the energy. We have carried out exhaustive studies of about 100 different combinations of clusters, cells, and boundary conditions to fully investigate the low-energy landscape of states. These detailed results are presented in (30). To characterize the orders, we use the local hole density Embedded Image, magnetic moment Embedded Image, and pairing order Embedded Image (i adjacent to j).

Characterizing the ground state at 1/8 doping

Using the above methods, we carried out calculations for the ground state of the 2D Hubbard model at 1/8 doping at U/t = 8. The first check of reliability is the independent convergence of the methods for the energy per site. Although the quality of the ground-state energy may be a poor proxy for the quality of the corresponding state when the overall accuracy is low (as there are always many degenerate states far above the ground state), calculations with well-converged energies tightly constrain the ground-state order, as any degeneracies must be below the energy-convergence threshold. Figure 1 shows the best energy estimate for the ground state from the different methods (30). The two different DMRG formulations (in real-space and hybrid basis) are in good agreement, providing a strong independent check of the calculations; in subsequent figures, we report only the single consistent result. Note that the error bars for AFQMC, DMRG, and DMET do not reflect the uncontrolled systematic errors in the methods. However, as described above, the systematic errors can be estimated by cross-checks between the methods. For example, DMRG and AFQMC calculations on finite clusters with identical boundary conditions provide an estimate of the small CP bias [see (30,33)] consistent with the difference in the DMRG and AFQMC energies in Fig. 1; similarly, AFQMC extrapolations to the thermodynamic limit indicate that the DMRG energies are essentially converged with respect to cylinder width.

Fig. 1 Ground-state energies.

Best estimates of ground-state energy for the 1/8-doped 2D Hubbard model at U/t = 8 from DMET, AFQMC, iPEPS, and DMRG in units of t. Inset shows best estimates of ground-state energy for the half-filled 2D Hubbard model at U/t = 8. Here and elsewhere, error bars indicate only the estimable numerical errors of each method; uncontrolled systematic errors are not included. For details, see (30).

There is good agreement between all the methods, and all energies lie in the range –0.767 ± 0.004t. If, for a typical HTSC material, we estimate t ~ 3000 K, then this corresponds to a range of about ±10K per site, or ±100K per hole. For a numerical comparison, this is also more than an order of magnitude lower than the temperatures accessible in simulations using finite-temperature methods in the thermodynamic limit in this part of the phase diagram, indicating that we are potentially accessing different physics (24, 29). Shown in the inset are the corresponding best estimates at half-filling from the same methods, where the spread in energies is less than 0.001t. This illustrates the substantially greater numerical challenge encountered in the underdoped region. Nonetheless, the accuracy and agreement reached here represent a 10-fold improvement over recent comparisons of numerical methods at this point in the phase diagram (24).

Ground-state stripe order

For all the methods used, the lowest energies shown in Fig. 1 correspond to a vertically striped state. This is a codirectional charge and spin-density wave state, with the region of maximum hole density coinciding with a domain wall in the antiferromagnetism. As mentioned, unidirectional stripes of various kinds are a long-standing candidate order in the doped Hubbard and related models. Hartree-Fock calculations give filled stripes (i.e., one hole per unit cell of the charge order) in both vertical and diagonal orientations, whereas one of the first applications of the DMRG to 2D systems found strong evidence for half-filled stripes in the t-J model (13). Finally, one of the earliest examples of inhomogeneity in doped HTSCs was the static half-filled stripes in LaSrCuO at 1/8 doping (40).

The convergence to the same inhomogeneous order in the ground state in the current study, from multiple methods with very different approximations, strongly suggests that stripes indeed represent the true ground-state order of the Hubbard model in the underdoped regime and further highlights the accuracy that we achieve with different techniques. However, the stripe order that we find has some unusual characteristics. We return to the details of the stripe order, its associated physics, and its relationship with experimentally observed stripes further below. First, however, we examine the possibility of other competing metastable states.

Competing states: Uniform d-wave state

Recent work using iPEPS and DMET on the t-J and Hubbard models suggested close competition between a uniform d-wave superconducting ground state and a striped order (18, 19). Uniform states did not spontaneously appear in any of our calculations, which indicates that they lie higher in energy than do striped orders. We found that we could stabilize a uniform d-wave state in the DMET calculations by constraining the impurity cluster to a 2 site–by–2 site or Embedded Image site–by–Embedded Image site geometry and, in the iPEPS calculations, by using a 2 site–by–2 site unit cell. DMET calculations on similarly shaped larger clusters (such as a 4 site–by–4 site cluster) spontaneously broke symmetry to create a nonuniform state. From these calculations, we estimate that the uniform state lies ~0.01t above the lowest-energy state and is not competitive at the energy resolution we can now achieve (30).

Competing states: Other short-range orders

Although other types of orders have been proposed in the underdoped region, such as spiral magnetic phases (20, 41) and checkerboard order (42), we find no evidence for other kinds of short-range orders at this point in the phase diagram. The lack of checkerboard order, which would easily fit within the large clusters in our simulations (e.g., up to 64 site–by–6 site in the DMRG calculations), appears to rule them out as low-energy states, in agreement with earlier DMRG simulations on the t-J model (43). Though we cannot rule out incommensurate orders, we have found that the variation of energy with unit-cell wavelength (see below) is quite smooth, and thus we do not expect a dramatic energy gain from incommensurability. We note that studies that have found incommensurate magnetic orders have focused on smaller values of U (20).

Diagonal versus vertical stripes

We find the ground-state order to be a vertical stripe–type order, but other studies of stripes indicate that different orientations can form (44). On short-length scales, the relevant question is whether diagonal stripes [with a (π, π) wave vector] are competitive with vertical stripes [with a (0, π) wave vector]. With the boundary conditions used in this work, diagonal stripes would be frustrated in the DMRG and AFQMC calculations and thus did not spontaneously appear. To stabilize diagonal stripes in the DMET and iPEPS calculations, we used tilted Embedded Image site–by–Embedded Image site impurity clusters (n = 2 or n = 5) for DMET, and a 16 site–by–16 site simulation cell with 16 independent tensors in iPEPS. The 16 site–by–16 site iPEPS cell gave a diagonal stripe (Fig. 2) that was substantially higher in energy than the vertical stripe, by 0.009t. The DMET cluster gave rise to a frustrated diagonal order that we also estimate to be higher in energy by ~0.005t (30). Although it is likely that the orientation of the stripe will depend on doping and coupling, vertical stripes appear to be clearly preferred at this point in the phase diagram.

Fig. 2 Competing states.

Shown are the energies of important competing states relative to the striped ground state from DMET and iPEPS and the sketches of the corresponding orders. (A) Relative energy of competing states in units of t compared to the vertically striped state. Charge, spin, and pairing orders of the uniform d-wave state from (B) DMET [blue circle in (A)] and (C) iPEPS [green squares in (A)]. (D) Charge and spin orders of the diagonally striped state from iPEPS. Note that the spins are flipped in the neighboring supercells. [For (B), (C), and (D), circle radius is proportional to hole density, arrow height is proportional to spin density (red and blue arrows indicate spin direction), and bond width (blue and tan lines) is proportional to pairing density]. For more details, see (30).

Ground-state stripes: Detailed analysis

We now return to a more detailed discussion of the vertical-stripe order found in the ground state. Within the family of vertical stripes, we can consider questions of wavelength (charge and spin periodicity), type and strength of charge and spin modulation (e.g., bond- versus site-centered), and coexistence with pairing order.

We first discuss the wavelength λ. At 1/8 doping, the filling of the stripe is related to the wavelength by λ/8. As described, we can access different wavelength metastable stripes and their relative energetics by carefully choosing different total cluster dimensions and boundary conditions (in the DMRG and AFQMC calculations) or unit cell (iPEPS) and impurity (DMET) sizes (30). Figure 3 shows the energy per site of the stripe versus its wavelength λ for the multiple methods. Earlier DMRG calculations on the Hubbard model had focused on λ = 4 (half-filled stripes), which are seen in HTSCs (13, 14), but we now observe that these are relatively high in energy. A striking feature is that for λ = 5 to 8 the energies are nearly degenerate. This is clearly seen in the DMET data, where stripes of all wavelengths can be stabilized, as well as from the averaged energy of the methods between λ = 5 to 8 (stars in Fig. 3). The energy difference between the λ = 5 and λ = 8 stripe in the different methods is estimated to be between 0.0005t (DMRG) and 0.0041t (iPEPS). This suggests that the magnetic domain walls can fluctuate freely, consistent with proposals for fluctuating stripes. In particular, the stripes may be distorted at a small cost over long length scales.

Fig. 3 Wavelength of the vertical-stripe order.

Energies of stripes with different wavelengths relative to that of the λ = 8 stripe from DMET, AFQMC, iPEPS, and DMRG in units of t. To aid readability, the data points are shifted horizontally. Inset shows relative energies of stripes with different wavelengths from UHF, with an effective coupling U/t = 2.7. For details of calculations, see (30).

Although the different wavelengths are nearly degenerate, there appears to be a slight minimum near wavelength λ = 8 (a filled stripe) in all the methods. Very recently, similar filled stripes have been reported as the ground state in part of the frustrated t-J model phase diagram (45). The wavelength λ = 9 appears much higher in energy in both DMET and DMRG. In the DMRG calculations, the λ = 9 state was not even metastable, as boundary conditions and initial states were varied, so the high-energy state shown was forced with a static potential. The AFQMC results show a much weaker dependence on wavelength for longer wavelengths; for example, the λ = 8 and λ = 10 stripe energies per site appear to be within 0.0015t. However, when a mixture of the λ = 8 and λ = 10 stripe states is set up on a 40 site–length cluster that is commensurate with both, the state that survives is the λ = 8 stripe, suggesting a preference for this wavelength. The increase in energy at wavelengths λ > 8 coincides with unfavorable double occupancy of the stripe. This simple interpretation is supported by a mean-field [unrestricted Hartree-Fock (UHF)] calculation with an effective interaction U/t = 2.7 chosen within the self-consistent AFQMC procedure (Fig. 3, inset). The mean-field result shows a clear minimum at a wavelength λ = 8 vertical stripe. [This requires the use of an effective U/t; at the bare U/t = 8, mean-field theory would produce a diagonal stripe (46).] The correspondence between the energies and densities in the effective mean-field and correlated calculations suggests that mean-field theory with a renormalized interaction may be surprisingly good at describing the energetics of stripes. However, mean-field theory appears to somewhat underestimate the degeneracy of the stripes as a function of wavelength, particularly at shorter wavelengths.

The vertical-stripe order for the λ = 8 stripe from the different methods is depicted in Fig. 4. We show the full period (16) for the spin modulation. The stripe is a bond-centered stripe in the AFQMC, DMRG, and DMET calculations. In the iPEPS calculation, the stripe is nominally site centered. In all the calculations, the width of the hole domain-wall spans several sites, blurring the distinction between bond- and site-centered stripes, and we conclude that the energy difference between the two is very small. There is substantial agreement in the order observed by the different numerical techniques, with only some differences in the modulation of the hole and spin densities.

Fig. 4 Charge and spin orders.

Shown are sketches of the charge and spin orders in the λ = 8 stripes from (A) DMET, (B) AFQMC, (C) iPEPS, and (D) DMRG. The local magnetic moments and hole densities are shown above and below the order plots, respectively. (Circle radius is proportional to hole density, and arrow height is proportional to spin density.) The gray dashed lines represent the positions of maximum hole density and the magnetic domain wall. For more details, see (30).

For even wavelength stripes, the spin wavelength must be twice that of the charge modulation to accommodate the stripe as well as the antiferromagnetic order. At odd wavelengths, site-centered stripes appear in all the calculations, and here, charge and spin order can have the same wavelength. [This odd-even alternation does not affect the peaks of the structure factor near (π, π); see (30).]

Pairing order, fluctuations, and superconductivity

A key question is whether pairing order coexists with stripe order. Previous work on the t-J model with iPEPS found coexisting d-wave order for partially filled (λ < 8) stripes. We did not find d-wave order in the Hubbard λ = 8 stripe with any technique. However, d-wave order can be found at other wavelengths. For example, for λ = 5 and λ = 7 stripes, iPEPS produces d-wave order along the bonds (see Fig. 5), with a maximum d-wave expectation value of 0.026 and 0.021, respectively. DMRG calculations with pinning pairing fields on the boundary for a 32 site–by–4 site cylinder also find d-wave order, with a maximum d-wave order of 0.025, consistent with the iPEPS results. In the DMET calculations, the lowest energy λ = 5 stripe has no d-wave order; however, at slightly higher energy (~0.003t), a λ = 5 state similar to the iPEPS stripe can be found with coexisting d-wave order, but with a substantially smaller maximum order parameter of 0.01. Overall, our results support the coexistence of modulated d-wave order with the striped state, although the strength of pairing is dependent on the stripe wavelength and filling. The pairing modulation that we find (Fig. 5) is in-phase between cells. Other kinds of pairing inhomogeneities, such as pair density waves, have also been discussed in the literature (6).

Fig. 5 d-wave pairing.

Shown are metastable stripe states with d-wave pairing from iPEPS, DMET, and DMRG. (A) and (B) iPEPS stripes with λ = 5 and λ = 7. (C) DMET metastable λ = 5 stripe with pairing. (Circle radius is proportional to hole density, arrow height is proportional to spin density, and bond width is proportional to pairing density). (D) DMRG pairing order parameters on a 32 site–by–4 site cylinder. The positive values are from the vertical bonds, and the negative values are from the horizontal bonds. The x axis is the site number along the long axis of the cylinder. For details, see (30).

It has long been argued that fluctuating stripes could promote superconductivity (4749). Our work provides some support for this conjecture, as there is a low-energy scale associated with the deformation of stripe wavelength as well as evidence for coupling between the wavelength and the pairing channel. We can imagine fluctuations in wavelength both at low temperatures as well as in the ground state. In the latter case, this could lead to a stripe-liquid ground state rather than a stripe crystal. Our calculations are consistent with both possibilities.

Varying the coupling

We may also ask whether the U/t = 8, 1/8 doping point is an anomalous point in the Hubbard phase diagram and if, for example, moving away from this point would cause the unusual stripe compressibility (with respect to wavelength at fixed doping) to be lost. In Fig. 6, we show the energies of various striped states and the uniform state at U/t = 6 and U/t = 12, 1/8 doping, computed using AFQMC, DMET, and DMRG. At both couplings, the stripes around wavelength λ = 8 are nearly degenerate, with the degeneracy increasing as the coupling increases. At U/t = 6, we find that the ground state is a filled stripe state with wavelength λ = 8, with a larger energy stabilization than at U/t = 8. The trend is consistent with the state observed at U/t = 4, with a more sinusoidal spin-density wave, more delocalized holes, and a more pronounced minimum wavelength (17). At U/t = 12, we find a filled stripe with AFQMC and DMRG (width of six sites), but DMET and DMRG on a narrower cylinder (width of four sites) find λ = 5 and λ = 6. The similarity of the DMET and DMRG (width of four sites) data suggests that the shorter wavelength is associated with a finite-width effect. We note that two-thirds-filled stripes consistent with λ = 5 and λ = 6 were also seen in earlier DMRG studies on six-site-width cylinders (15), but a more detailed analysis shows that the filled stripe becomes favored when extrapolated to infinite bond dimension (30). Thus, we conclude that the ground state at U/t = 12 is also the λ = 8 stripe, although stripes of other wavelengths become even more competitive than at U/t = 8. Overall, the similarity in the physics over a wide range of U/t indicates that striped orders with low-energy fluctuations of domain walls remain a robust feature in the moderate to strongly coupled underdoped region.

Fig. 6 Varying the interaction strength.

Relative energies of stripe states (versus wavelength) and the uniform d-wave state at 1/8 doping for (A) weaker and (B) stronger couplings. For details, see (30).

Connection to stripe order in HTSCs

In HTSCs, the accepted stripe wavelength at 1/8 doping (e.g., in LaSrCuO) is λ ≈ 4.3 (close to half-filled) (40). However, we find that the λ = 4 stripe is not favored in the 2D Hubbard model for the coupling range (U/t = 6 to 12) normally considered most relevant to cuprate physics. This implies that the detailed charge ordering of real materials arises from even stronger coupling or, more likely, quantitative corrections beyond the simple Hubbard model. With respect to the latter, one possibility is long-range hopping (such as a next-nearest-neighbor hopping), which has been seen to change the preferred stripe wavelength in the frustrated t-J model (45). Another possibility is the long-range Coulomb repulsion. Long-range repulsion can play a dual role, in both driving charge inhomogeneity as well as smoothing it out. In the Hubbard model, where stripes naturally form, the latter property can help drive the ground state toward shorter stripe wavelengths. We have estimated the effect of the long-range interactions on the stripe energetics by computing the Coulomb energy of the charge distributions in Fig. 4. We use a dielectric constant of 15.5 [in the range proposed for the cuprate plane (50)]. This gives a contribution favoring the shorter wavelength stripes that is on the order of ~0.01t per site for the λ = 4 versus λ = 8 stripe (30). Although this is only an order-of-magnitude estimate, it is on the same energy scale as the stripe energetics in Fig. 2 and thus provides a plausible competing mechanism for detailed stripe physics in real materials.


In this work, we have used state-of-the-art numerical methods to determine the ground state of the 1/8 doping point of the 2D Hubbard model at moderate to strong coupling. Through careful convergence of all the methods, and exhaustive cross-checks and validations, we are able to eliminate several of the competing orders that have been proposed for the underdoped region in favor of a vertically striped order with wavelength near λ ≈ 8. The striped order displays a remarkably low-energy scale associated with changing its wavelength, which implies strong fluctuations either at low temperature or in the ground state itself. This low-energy scale can roughly be accounted for at the mean-field level with a strongly renormalized U. We find coexisting pairing order with a strength dependent on the stripe wavelength, indicating a coupling of stripe fluctuations to superconductivity. The stripe degeneracy is robust, as the coupling strength is varied.

It has long been a goal of numerical simulations to provide definitive solutions of microscopic models. Our work demonstrates that even in one of the most difficult condensed matter models, such unambiguous simulations are now possible. In so far as the 2D Hubbard model is a realistic model of high-temperature superconductivity, the stripe physics observed here provides a firm basis for understanding the diversity of inhomogeneous orders seen in the materials, as well as a numerical foundation for the theory of fluctuations and its connections to superconductivity. However, our work also enables us to see the limitations of the Hubbard model in understanding real HTSCs. Unlike the stripes at this doping point in real materials, we find filled stripes rather than near-half-filled stripes. Given the very small energy scales involved, terms beyond the Hubbard model, such as long-range Coulomb interactions, will likely play a role in the detailed energetics of stripe fillings. The work we have presented provides an optimistic perspective that achieving a comprehensive numerical characterization of more-detailed models of the HTSCs will also be within reach.

Supplementary Materials

Materials and Methods

Figs. S1 to S41

Tables S1 to S10

References (5167)

References and Notes

  1. See supplementary materials.
  2. We construct the trial wave function with desired wavelength by solving the noninteracting Hamiltonian with pinning fields in the whole system with the same structure.
Acknowledgments: Work performed by B.-X.Z., C.-M.C., M.-P.Q., H.S., S.R.W., S.Z., and G.K.-L.C. was supported by the Simons Foundation through the Simons Collaboration on the Many Electron Problem. S.R.W. acknowledges support from the NSF (DMR-1505406), as do S.Z. and H.S. (DMR-1409510). M.-P.Q. was also supported by the U.S. Department of Energy (DOE) (DE-SC0008627). G.K.-L.C. acknowledges support from a Simons Investigators Award and the DOE (DE-SC0008624). DMET calculations were carried out at the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by DE-AC02-05CH11231. AFQMC calculations were carried out at the Extreme Science and Engineering Discovery Environment, supported by the NSF (ACI-1053575); at the Oak Ridge Leadership Computing Facility at Oak Ridge National Lab; and at the computational facilities at the College of William and Mary. P.C. was supported by the European Research Council under the European Union’s Horizon 2020 research and innovation program (grant no. 677061). G.E. and R.M.N. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) through grant no. NO 314/5-1 in Research Unit FOR 1807. Data used in this work are in the supplementary materials and online at The DMET calculations were performed with the DMET library, available online at The real-space DMRG calculations were performed with ITensor, available online at Access to the other computer codes can be arranged with the authors upon reasonable request: For the AFQMC code, please contact S.Z.; for the hybrid DMRG code, please contact R.M.N.; and for the iPEPS code, please contact P.C.
View Abstract

Navigate This Article