Stable diagonal stripes in the t-J model at $\bar{n}_h$=1/8 doping from fPEPS calculations

We investigate the 2D t-J model at a hole doping of $\bar{n}_h$=1/8 using recently developed high accuracy fermionic projected entangled pair states(fPEPS) method. By applying stochastic gradient descent method combined with Monte Carlo sampling technique, we obtain the ground state hole energy $E_{\rm hole}$=-1.6186 for $J/t$=0.4. We show that the ground state has stable diagonal stripes instead of vertical stripes with width of 4 unit cells, and stripe filling $\rho_l$=0.5. We further show that the long range superconductivity order is suppressed at this point.

The high-Tc superconductivity [1,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 diagrams [3]. Specifically, the hole doping nearn h =1/8 provides an ideal system to experimentally study the ground state of the pseudogap [4] which is one of the most salient phenomena in high-Tc superconductivity [5]. Near this doping level, the charge and spin stripe orders are observed in some cuprate compounds, e.g., La 1.875 Ba 0.125 CuO 4 , by various experimental techniques, including angle-resolved photoemission and scanning tunneling microscopy [4], neutron and xray scattering [6], etc.
It is widely believed that the physics of superconductivity could be understood as doped Mott insulators [5], which could be described by the 2-dimensional Hubbard model [7,8] and the t-J model [9], the strong coupling limit of Hubbard model. However, the theoretical results about the ground state near hole dopingn h =1/8 in the t-J model are still highly controversial [10][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 debates. [10][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 few Lanczos steps [16], suggest that the ground state atn h =1/8 is homogeneous without stripes order. These results are contradictory to the results of the early DMRG calculations [14,[21][22][23]. More recently, iPEPS with full update calculations [17,18] suggest that the ground state has stable stripes. Nevertheless, the calculations [17,18] suggest that the uniform phase is energetically very close to the stripe phase, and the energy difference become 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 accuracy [34][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 doping level n 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 E h = -1.6186 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 previous works [10,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.
The t-J model is defined on a two-dimensional square lattice as, where i, j are the nearest-neighbor sites. c i,σ (c † i,σ ) is the electron annihilation (creation) operator of spin σ (σ =↑, ↓) on site i, whereas n i = σ c † i,σ c i,σ and S i are the electron number and the spin-1/2 operators respectively. Double occupations are not allowed. We solve the model by using recently developed fPEPS [17, 30-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 techniques. [34,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-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 techniques [39,40] 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 D c to converge the results. [36] In this work, open boundary conditions (BC) are used. We focus on the parameters of t=1 and J/t = 0.4 with hole doping ofn h = 1 8 . Figure 1 depicts the ground state energies of t-J model at hole dopingn h =1/8, for different system sizes L 1 ×L 2 , ranging from 4×4 to 12×12. In all calculations, the virtual bond dimension D is fixed to 12, whereas the truncated dimension is set to D c =48 to ensure that the energies are well converged at the given D. As shown in our previous work, [36] for t-J model atn h =1/8, usually the choice of D c =3D=36 can ensure convergence, which is only weakly size dependent. For the 4×4 system, the ground state energy obtained by our calculation is E fPEPS =-0.56428, compared with the exact energy E ED =-0.56436, where the energy difference is about 1×10 −4 .
We extrapolate the ground state energies to the thermodynamic limit using a second-order polynomial fitting on √ L 1 L 2 . The extrapolated ground state energy in the thermodynamic limit is E ∞ =-0.6701. The corresponding energy per hole is defined as E hole = [E(n h )−E 0 ]/n h , where E 0 =-0.467775 is the energy at zero doping, taken from Ref. 41. We compare the ground state hole energies atn h =1/8 obtained by various methods in Table I. The hole energy we obtained is E fPEPS hole =-1.619, which is slightly lower than the hole energy obtained from DMRG calculation, [17] E DMRG hole =-1.612 with χ extrapolated to ∞. The result is also lower than the one from iPEPS SU calculations E iPEPS hole =-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 iPEPS hole =-1.578 forn h =0.12, [18] which is expected to has lower (more negative) hole energy than that of n h =1/8. As a comparison, the recent variational QMC simulation gives E QMC hole =-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 Energy VQMC Ref. [16] iPEPS Ref. [17] DMRG Ref. [17] (L L )  close to E QMC hole . 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 method [13] in Fig. S1 of the Supplementary materials (SM [42]). They are in remarkably good agreement. The calculations suggest that ground state of the t-J model with hole filling ofn h = 1 8 on the 4×4 lattice is in an 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 become non-uniform when increase the size of the system. For the 4×8 system, the ground state of the system is not uniform any more. The holes clusters are more localized in the center of the lattice without local magnetic order [see Fig. S2(a) in the SM [42]]. When the lattice size  is further increased to 6×8 and larger, the holes form stripes along the diagonal direction [see Fig. S2(b,c) in the SM [42]]. 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 diagonal direction on the antiferromagnetic background, with a π phase-shifted magnetic order across the domain wall. [14] To investigate the structure of the diagonal stripe states, we plot the hole density n h i,j =1-n i,j , where n i,j is the average number of electrons on site (i, j), and staggered magnetization (−1) i−1 S z i,j 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 sitecentered nature of the hole stripes is evident from the hole density n h i,j (green line). These hole and spin pattern with 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 has a hole filling ρ l =W ·n h =0.5, where W is the width of the domain wall, i.e., half filling. [14] To further test the robustness of the diagonal stripes against the size of the system, we simulated on lattices of different sizes (see Table. S1 in the SM [42]), 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.
The t-J model atn h =1/8 has been intensively investigated by various methods, and the results are still highly controversial. [10][11][12] On the one hand, the recent variational quantum Monte Carlo (QMC) simulations combined with few Lanczos steps [16] suggest that the ground state atn h =1/8 is homogeneous without stripes order. Its hole energies are very close to those obtained from DMRG [14,[21][22][23] and recent iPEPS calculations [18]. On the other hand, DMRG calculations show that the ground state has stable stripes. [14,18,[21][22][23] Our results support that the stripe phase is stable, which is in agreement with the DMRG calculations. [14,18,[21][22][23] In DMRG calculations, the stripes are further characterized as the site-centered vertical stripes, and the width of the stripes are 4 atn h ∼ 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 BCs used in the calculations. In the DMRG calculations, a periodic BC in y-axis, and an open BC in x-axis are used, which favors the vertical stripes along the y direction. [14] 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 calculations [17,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, andn h ∼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 stripes. [17,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 L 2 ) tensors were independent. With these constrains, 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 L 1 × L 2 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 require 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 D, [18] 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 method [43].
As a comparison, we also calculate the anisotropic t-J model with t x /t y =0.85, and J x /J y =0.85 2 following Ref. 18. We show the result of J x = 0.4t x , t y = 0.85t x and J y = 0.289t x for different sizes in Fig. S3 of the SM [42]. Compared with the isotropic case, the stripe orients along the bonds with stronger couplings, which is in agreement with previous results, [18,44,45]. The results can be understood as the kinetic energy can be effectively lowered by hopping along this directions. The anisotropic interaction converts the site-centered diagonal stripes into bond-centered vertical ones in these simulations.
To further investigate the relationship between the stripe order and the superconductivity, we calculate the hole pair correlation functions, which are defined as, where s, d denote the s-or d-wave paring. The superconductivity order parameter ∆ s,d (r i ) is defined as, with r i being the coordinate at site i, with "+" for s-wave and "−" for d-wave paring. In Fig. 4, we show both the sand the d-wave pair correlation functions P s,d (i, j) with r i fixed at (6,2) and r j changed from (6,2) to (6,12) in the 12× 12 lattice. To obtain highly accurate results, D c =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 atn h =1/8 decay rather quickly with distance. Even though the pair correlations can be fitted roughly by power law decay functions, i.e., P s,d (r) ∼ r −α , with α ∼ 4.9 and 4.4 for s-wave and dwave pairing respectively, since α ≫1, the long range order of superconductivity is suppressed atn h =1/8.
To summarize, we investigate the ground state of t-J model at hole dopingn 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 above calculations the stripes are vertical. We further show that the long range order of superconductivity is suppressed in this phase.