Introduction

The high-Tc superconductivity1,2 is probably one of the most exciting and also challenging open problems in condensed matter physics. The strong coupling between the spin and charge degrees of freedom leads to various competing orders at low temperature, resulting in rich phase diagrams3. Specifically, the hole doping near \({\bar{n}}_{{\rm{h}}}=1/8\) provides an ideal system to experimentally study the ground state of the pseudogap4 which is one of the most salient phenomena in high-Tc superconductivity5. Near this doping level, the charge and spin stripe orders are observed in some cuprate compounds, e.g., La1.875Ba0.125CuO4, by various experimental techniques, including angle-resolved photoemission and scanning tunneling microscopy4, neutron and X-ray scattering6, etc.

It is widely believed that the physics of superconductivity could be understood as doped Mott insulators5, which could be described by the two-dimensional Hubbard model7,8 and the t–J model9, the strong coupling limit of Hubbard model. However, the theoretical results about the ground state near hole doping \({\bar{n}}_{{\rm{h}}}=1/8\) in the t–J model are still highly controversial10,11,12. The question about whether the ground state has the stable stripe order, and the relation between the superconductivity and the stripe order are under intensive debates10,11,12,13,14,15,16,17,18. Early works on this issue have been reviewed in refs. 19,20. Very recently, variational quantum Monte Carlo (vQMC) simulations combined with few Lanczos steps16, suggest that the ground state at \({\bar{n}}_{{\rm{h}}}=1/8\) is homogeneous without stripes order. These results are contradictory to the results of the early density matrix renormalization group (DMRG) calculations14,21,22,23. More recently, infinite projected entangled pair states (iPEPS) with full update calculations17,18 suggest that the ground state has stable stripes. Nevertheless, the calculations17,18 suggest that the uniform phase is energetically very close to the stripe phase, and the energy difference becomes even smaller with increasing bond dimension. Therefore, it is hard to determine what the true ground state is unless fully converged calculations are performed.

The projected entangled pair states method (PEPS)24,25,26,27,28,29, and its generalization to fermionic systems (fPEPS)30,31,32,33 provide systematically improvable variational wave functions for many-body problems. In recent works, we developed a gradient method combined with Monte Carlo sampling techniques to optimize the (f)PEPS wave functions with controlled accuracy34,35,36. This method significantly reduces the scaling with respect to the bond dimension D, thereby allowing a much larger bond dimension to be used, resulting in highly accurate and converged results for large finite systems. In this work, we apply this recently developed and highly accurate fPEPS method to explore the true ground state of the t–J model at the doping level \({\bar{n}}_{{\rm{h}}} = 1/8\). The computational results allow us to shed some new light on this long-standing open problem. From our computation, we obtained the hole energy Eh = −1.621 for J/t = 0.4 in the thermodynamic limit. Remarkably, we find that the ground state of the t–J model at 1/8 hole doping has stable stripes that are along the diagonal directions instead of the vertical direction suggested by the previous works10,12,14,17,18,37 with stripe hole filling ρl = 0.5. We further show that the long-range superconductivity order at this point is suppressed.

Results and discussion

The model

The t–J model is defined on a two-dimensional square lattice as

$$H=-t{\sum \limits_{\langle i,j\rangle ,\sigma }}({c}_{i,\sigma }^{\dagger }{c}_{j,\sigma }+H.c.)+J{\sum \limits_{\langle i,j\rangle}}({{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}-\frac{1}{4}{n}_{i}{n}_{j})$$
(1)

where 〈i, j〉 are the nearest-neighbor sites. ci,σ (\({c}_{i,\sigma }^{\dagger }\)) is the electron annihilation (creation) operator of spin σ (σ = ↑,↓) on site i, whereas \({n}_{i}={\sum }_{\sigma }{c}_{i,\sigma }^{\dagger }{c}_{i,\sigma }\) and Si are the electron number and the spin-1/2 operators, respectively. Double occupations are not allowed. We solve the model by using recently developed fPEPS17,30,31,32,33,36,38 method, which is described in the Methods section.

Ground state energies

For the 4 × 4 system, the ground state energy obtained by our calculation is EfPEPS = −0.56428, compared with the exact energy EED = −0.56436, where the energy difference is about 1 × 10−4.

Fig. 1: The extrapolation of the ground state energies.
figure 1

The ground state energies of t–J model with t = 1, and J/t = 0.4 at hole doping \({\bar{n}}_{{\rm{h}}}=\frac{1}{8}\). The red squares represent the energies on different lattice sizes. The energies are extrapolated to the thermodynamic limit via a second-order polynomial function of \(\sqrt{{L}_{1}{L}_{2}}\). The green, blue, and black lines are the ground state energies obtained by QMC, iPEPS simple update, and DMRG methods, respectively.

Full size image

Figure 1 depicts the ground state energies of t–J model at hole doping \({\bar{n}}_{{\rm{h}}}=1/8\), for different system sizes L1 × L2, ranging from 6 × 8 to 12 × 12, and these energies are given in Table S239. To reduce the boundary effects in the extrapolation, we exclude the data of two smallest sizes: the 4 × 4 and 4 × 8 lattices. We extrapolate the ground state energies to the thermodynamic limit via a second-order polynomial fitting on \(\sqrt{{L}_{1}{L}_{2}}\). The extrapolated ground state energy in the thermodynamic limit is E∞ = −0.6704 (E∞ = −0.6711 if the two smallest lattices are included). The corresponding energy per hole is defined as \({E}_{{\rm{hole}}}=[E({\bar{n}}_{{\rm{h}}})-{E}_{0}]/{\bar{n}}_{{\rm{h}}}\), where E0 = −0.467775 is the energy at zero doping, taken from ref. 40. We compare the ground state hole energies at \({\bar{n}}_{{\rm{h}}}=1/8\) obtained by various methods in Table 1. The hole energy we obtained is \({E}_{{\rm{hole}}}^{{\rm{fPEPS}}}=-1.621\) (\({E}_{{\rm{hole}}}^{{\rm{fPEPS}}}=-1.626\) if the two smallest lattices are included), which is lower than the hole energy obtained from DMRG calculation17, \({E}_{{\rm{hole}}}^{{\rm{DMRG}}}=-1.612\) with χ extrapolated to ∞. The result is also lower than the one from iPEPS SU calculations \({E}_{{\rm{hole}}}^{{\rm{iPEPS}}}=-1.593\)17, which was obtained by extrapolating bond dimension D to ∞. The recent iPEPS full update calculations with D = 14 give the hole energy \({E}_{{\rm{hole}}}^{{\rm{iPEPS}}}=-1.578\) for \({\bar{n}}_{{\rm{h}}}=0.12\)18, which is expected to have lower (more negative) hole energy than that of \({\bar{n}}_{{\rm{h}}}=1/8\). As a comparison, the recent variational QMC simulation gives \({E}_{{\rm{hole}}}^{{\rm{QMC}}}=-1.546\)16. We note that the ground state energy obtained in this work is significantly lower than the ground state energies obtained by DMRG and iPEPS calculation before extrapolation, which are close to \({E}_{{\rm{hole}}}^{{\rm{QMC}}}\).

Table 1 The ground state energies.
Full size table

Diagonal stripe orders

We now take a closer look at the t–J model. The hole density and magnetization of the 4 × 4 lattice obtained from fPEPS are compared with those obtained by diagonalization method13 in Fig. S139. They are in remarkably good agreement. The calculations suggest that the ground state of the t–J model with hole filling of \({\bar{n}}_{{\rm{h}}}=\frac{1}{8}\) on the 4 × 4 lattice is in a uniform phase with virtually no local magnetization (the local magnetization is less than 10−4), which is in good agreement with the conclusions of ref. 13.

Would the uniform state be stable when the size of the system is increased? Unfortunately, the hole distribution becomes non-uniform when increasing the size of the system. For the 4 × 8 system, the ground state of the system is not uniform anymore. The holes clusters are more localized in the center of the lattice without local magnetic order (see Fig. S2a39). When the lattice size is further increased to 6 × 8 and larger, the holes form stripes along the diagonal direction (see Fig. S2b, c39).

Fig. 2: The ground state hole density and spin moment.
figure 2

The ground state hole density and spin moment on the 12 × 12 lattice. The diameter of the circles represents the magnitude of holes density and the length of the arrow represents the local magnetic moments.

Full size image

Figure 2 depicts the ground state hole distribution and local magnetization of the 12 × 12 lattice. The sizes of the circles and arrows represent the magnitude of the hole density and local magnetic moments. The systems show clear stripe order along the diagonal direction on the antiferromagnetic background, with a π phase-shifted magnetic order across the domain wall14.

To investigate the structure of the diagonal stripe states, we plot the hole density \(\langle {n}_{i,j}^{{\rm{h}}}\rangle\) = 1 − 〈ni,j〉, where 〈ni,j〉 is the average number of electrons on site (i, j), and staggered magnetization \({(-1)}^{i-1}{S}_{i,j}^{z}\) along the diagonal direction perpendicular to the stripes in Fig. 3. The staggered magnetization (red line) shows a period of 8, whereas the hole density shows a period of 4. The site-centered nature of the hole stripes is evident from the hole density \(\langle {n}_{i,j}^{{\rm{h}}}\rangle\) (green line). These hole and spin patterns with a period of 4 and 8, which are robust for different size of systems, can be used to explain why the stripe order cannot be stabled in the small 4 × 4 and 4 × 8 systems, which are too small to accommodate such stripes. The stripes have a hole filling \({\rho }_{{\rm{l}}}=W\cdot {\bar{n}}_{{\rm{h}}}=0.5\), where W is the width of the domain wall, i.e., half filling14.

Fig. 3: The average hole density along the diagonal direction.
figure 3

The average hole density 〈nh〉 (green squares) and spin structure function (red circles) along the diagonal direction on the 12 × 12 lattice.

Full size image

To further test the robustness of the diagonal stripes against the size of the system, we simulated on lattices of different sizes (see Table S139), and aspect ratios. Except in some extreme cases, e.g., the width of the lattice is less than 4, we always obtain the diagonal stripes with similar spin and hole distribution structures. These results clearly demonstrate the diagonal stripes ground states are robust against the size and shape of the lattices. We also tried to induce the vertical stripes by applying zig-zag magnetic field on the boundary14. However, the stripes slightly distort after the magnetic field is switched off, and the energy of the state is about 0.0015 higher than the ground state with diagonal stripe orders in the 12 × 12 system.

The t–J model at \({\bar{n}}_{{\rm{h}}}=1/8\) has been intensively investigated by various methods, and the results are still highly controversial10,11,12. On the one hand, the recent variational quantum Monte Carlo (QMC) simulations combined with few Lanczos steps16 suggest that the ground state at \({\bar{n}}_{{\rm{h}}}=1/8\) is homogeneous without stripes order. Its hole energies are very close to those obtained from DMRG14,21,22,23 and recent iPEPS calculations18. On the other hand, DMRG calculations show that the ground state has stable stripes14,18,21,22,23.

Our results support that the stripe phase is stable, which is in agreement with the DMRG calculations14,18,21,22,23. In DMRG calculations, the stripes are further characterized as the site-centered vertical stripes, and the width of the stripes is 4 at \({\bar{n}}_{{\rm{h}}} \sim 1/8\)14, which are also in good agreement with our results. However, in our calculations, the stripes are along the diagonal direction in contrast to the vertical stripes obtained from DMRG calculations. This discrepancy may come from different boundary conditions (BCs) used in the calculations. In the DMRG calculations, a periodic BC in the y-axis, and an open BC in the x-axis are used, which favors the vertical stripes along the y direction14. One possible way to clarify this problem would be to perform DMRG calculations with periodic BC along the diagonal direction, to enforce a ground state with diagonal stripes, and compare the energy with that of vertical stripes.

Very recently, the iPEPS calculations17,18 also suggest that the ground state is a vertical stripe phase, where the width of stripes and stripe hole filling depend on the exchange parameter J. At J/t = 0.4, and \({\bar{n}}_{{\rm{h}}} \sim 1/8\), they obtain stripe filling ρl ~ 0.5, which is in agreement with DMRG calculations. They also compare the energies of diagonal stripes and vertical stripes, and it has been found that the diagonal stripes have somewhat higher energies than the vertical stripes17,18. However, the diagonal stripes obtained by the iPEPS calculations are very different from our cases. In ref. 18, the L × L (for L up to 11) supercells were used to obtain the diagonal stripe phase, but only L (instead of L2) tensors were independent. With these constraints, the resulting diagonal stripes are insulating with filling ρl = 1 holes per unit length, compared to ρl = 0.5 holes per unit length in our calculations. We remark that our calculations are unbiased in the sense that we do not have any constrains on the tensors. All L1 × L2 tensors are independent and free to change during the optimization. We always obtain the same stripe ground states for randomly chosen initial states. Another important difference between our method and the iPEPS method is that the iPEPS method directly works in the thermodynamic limit, which requires extremely large D to converge the results. In practice, such large D is infeasible, and therefore, the final results rely heavily on the extrapolation on the bond dimension D. In fact, it has been found that the energy differences between the uniform state and the stripes states become smaller and smaller with increasing D18, Therefore, no definite conclusion can be made based on their current numerical results. On the contrary, we work on finite systems, where the results can be fully converged with the given D. We then extrapolate the results to the thermodynamic limit by the well established finite-scaling method41.

As a comparison, we also calculate the anisotropic t–J model with tx/ty = 0.85, and Jx/Jy = 0.852 following ref. 18. We show the result of Jx = 0.4tx, ty = 0.85tx, and Jy = 0.289tx for different sizes in Fig. S339. Compared with the isotropic case, the stripe orients along the bonds with stronger couplings, which is in agreement with the previous results18,42,43. The results can be understood as the kinetic energy can be effectively lowered by hopping along these directions. The anisotropic interaction converts the site-centered diagonal stripes into bond-centered vertical ones in these simulations.

Superconductivity

To further investigate the relationship between the stripe order and the superconductivity, we calculate the hole pair correlation functions, which are defined as

$${P}_{{\rm{s,d}}}(i,j)=\left\langle {\Delta }_{{\rm{s,d}}}^{\dagger }\left({{\boldsymbol{r}}}_{i}\right){\Delta }_{{\rm{s,d}}}\left({{\boldsymbol{r}}}_{j}\right)+{\Delta }_{{\rm{s,d}}}\left({{\boldsymbol{r}}}_{i}\right){\Delta }_{{\rm{s,d}}}^{\dagger }\left({{\boldsymbol{r}}}_{j}\right)\right\rangle$$
(2)

where s, d denote the s- or d-wave paring. The superconductivity order parameter Δs,d(ri) is defined as

$$\begin{array}{rcl}{\Delta }_{{\rm{s,d}}}\left({{\boldsymbol{r}}}_{i}\right)={\sum \limits_\pm }\frac{1}{2}\left\{\left({c}_{i\uparrow }{c}_{i\pm \hat{x}\downarrow }-{c}_{i\downarrow }{c}_{i\pm \hat{x}\uparrow }\right)\pm \left({c}_{i\uparrow }{c}_{i\pm \hat{y}\downarrow }\right.\right. \left.\left.-{c}_{i\downarrow }{c}_{i\pm \hat{y}\uparrow }\right)\right\}\end{array}$$
(3)

with ri being the coordinate at site i, with “+” for s-wave and “−” for d-wave paring. In Fig. 4, we show both the s- and the d-wave pair correlation functions Ps,d(i, j) with ri fixed at site (6,2) and rj changed from (6,2) to (6,12) in the 12 × 12 lattice. To obtain highly accurate results, Dc = 6D is used to calculate the correlation functions. The dips in the correlation functions around r ~ 3–5 are presumably related to the hole stripe structures. As shown in the figure, the superconductivity order for both s-wave and d-wave pairing in the t–J model at \({\bar{n}}_{{\rm{h}}}\) = 1/8 decay rather quickly with distance. Even though the pair correlations can be fitted roughly by power law decay functions, i.e., Ps,d(r) ~ r−α, with α ~ 4.9 and 4.4 for s-wave and d-wave pairing respectively, since α ≫ 1, the long-range order of superconductivity is suppressed at \({\bar{n}}_{{\rm{h}}}=1/8\).

Fig. 4: The superconductivity pair correlation functions.
figure 4

The pair correlation functions for the d-wave and the s-wave superconductivity on 12 × 12 lattice. We fix ri = (6, 2) and scan rj = (6, 2) to rj = (6, 12), where r is defined as ∣rj − ri∣.

Full size image

To summarize, we investigate the ground state of t–J model at hole doping \({\bar{n}}_{{\rm{h}}}=1/8\), using the recently developed highly accurate fPEPS method. We obtain the most competitive ground state hole energy. We find that the ground state has stable stripes along the diagonal direction, with stripe hole filling ρl = 0.5. These results partially agree with recent DMRG and iPEPS calculations, in the sense that the ground state has stable stripes, except that in the above calculations the stripes are vertical. We further show that the long-range order of superconductivity is suppressed in this phase. The work provides a new scenario of the ground state of the long-standing open problem.

Methods

We solve the t–J model by using recently developed fPEPS17,30,31,32,33,36,38 method. The fPEPS wave functions are first optimized via an imaginary time evolution with simple update (SU)26 scheme, followed by gradient optimization combined with Monte Carlo sampling techniques34,36. U(1) symmetry is enforced during the calculations to conserve the number of electrons in the system. More details of the methods are discussed in refs. 34,35,36. It is well known that the environment effects are oversimplified in the SU method. Therefore, the use of SU may introduce large errors. However, the results can be used as a good starting point for the subsequent gradient optimization. The gradient optimization method treats the environment effects exactly with controllable errors, and therefore can obtain much more accurate results. The Monte Carlo sampling techniques44,45 are used to calculate the energies and their gradients, which may greatly reduce the complexity of the calculations, and allow us to use large virtual bond dimension D and truncation parameters Dc to converge the results36. In this work, open boundary conditions (BC) are used. Using a periodic BC may have less boundary effects. However, due to the extremely high computational costs [O(D8) for periodic BC vs O(D6) for open BC], it is still infeasible to treat periodic systems via fPEPS. We focus on the parameters of t = 1 and J/t = 0.4 with hole doping of \({\bar{n}}_{{\rm{h}}}=\frac{1}{8}\). In all calculations, the virtual bond dimension D is fixed to 12, and the truncated dimension is set to Dc = 56 (~4.7D) to ensure that the energies are well converged at the given D, where the errors due to Dc are less than 10−5.