Charge stripe and superconductivity tuned by interlayer interaction in a sign-problem-free bilayer extended Hubbard model

Abstract

Competing orders represent a central challenge in understanding strongly correlated systems. In this work, we employ projector quantum Monte Carlo simulations to study a sign-problem-free bilayer extended Hubbard model. In this model, a charge stripe phase, characterized by a peak at (0\pm 2\text{π}\delta ,0) is induced by highly anisotropic interlayer antiferromagnetic spin-exchange coupling J_z, and strongly suppressed by the decreasing spin-exchange anisotropy (i.e., introducing the spin-flip term J_{\mathrm{\perp }}); in contrast, the introduction of J_{\perp } favors the emergence of interlayer pairing superconductivity. We further demonstrate that the anisotropy of the interlayer spin-exchange directly governs the competition between these two phases, while the on-site interaction U plays a complex role in tuning both the charge stripe and superconductivity. Furthermore, against the background of charge stripes we discover the spin stripes, characterized by twice the period of charge stripes, indicating the profound internal connection between them. Our work identifies the key factors driving charge stripe formation, highlights the sensitivity of both the charge stripe and superconducting phases to interaction parameters, and thereby provides valuable insights into competing orders in strongly correlated systems.

Graphical abstract

Keywords

extended Hubbard model / stripe order / superconductivity / anisotropic spin-exchange

1 Introduction

One of the interesting phenomena in the Hubbard model is the competition between different phases [1–3]. This also involves the crucial open problem of whether the Hubbard model can describe high-T_c superconductors or not [4–6]. Despite the simplicity of the pure Hubbard model, there are controversies even in the phase diagram since the solutions of the Hubbard model in more than one dimension in the strong coupling region are hard to obtain [7]. Some experimental findings regarding high-T_c superconductivity in cuprate-based materials may not be totally reproduced by the pure Hubbard model, especially the details of charge ordering and superconductivity [8–11]. One comprehensive study using the constrained-path quantum Monte Carlo and density matrix renormalization group algorithm reveals that in the pure Hubbard model under strong coupling and at a special electron filling, the stripe phase is dominant and superconductivity is suppressed [12]. Furthermore, the periodicity of the charge stripe observed in experiments is not consistent with the numerical results from the pure Hubbard model [13]. This may imply that some extra terms should be introduced into the pure Hubbard model to stabilize superconductivity and reproduce the correct charge ordering.

For example, long-range hopping can impair the stripe phase in some parameter regions and yield the correct periodicity of charge ordering [13–17]. Long-range hopping can also capture the doping asymmetry in experimental data. The nearest neighbor interaction is another factor that is being investigated intensively. Proposed by the investigations of one-dimensional cuprate chains [18], the nearest-neighbor attractive interaction has been found to enhance the superconducting correlation and suppress the charge ordering [19–22]. Besides, the electron−phonon coupling is also considered for stripes [23]. Other studies have introduced stripe order explicitly, which can be achieved by adding a spatially modulated potential. The impacts of these modulations are complicated and it has been found that pairing strength and pairing symmetry can be affected [24–26]. Furthermore, the spin-order correlation undergoes a \text{π}-phase shift when crossing a stripe [27].

Above all, the charge stripe phase is considered to be a very important feature that affects superconductivity in the Hubbard model. A lot of work focuses on identifying factors that can suppress the stripe phase and thereby enhance superconductivity. Understanding the charge stripe is of key importance in explaining unconventional superconductivity. Various numerical methods have been developed to study Hubbard models, but most of them are only efficient under certain conditions [28–34]. For example, the determinant quantum Monte Carlo (DQMC) algorithm [35] and its projector variant, the projector quantum Monte Carlo (PQMC), struggle with sign problem in many cases of interest.

Sign problems occur in the quantum Monte Carlo algorithms when the Boltzmann factor is negative. This can happen after decoupling the interaction terms with the Hubbard−Stratonovich (HS) transformation and tracing out the fermion degrees of freedom. The sign problem is the main obstacle when using the DQMC and PQMC algorithms to study the rich physics in the Hubbard model. One approach to overcome the sign problem is by utilizing sign-problem-free models. Among these models, the bilayer model, which belongs to the Kramers-class of sign-free models, shows potential for uncovering many-body effects in strongly correlated systems. Previous work has found a transition from a Mott insulator to superconductivity, a transition analogous to that in cuprate superconductors [36]. Some other studies using a similar model have found a stable charge stripe phase with highly anisotropic interlayer ferromagnetic interactions [37, 38]. This demonstrates the value of this model in the investigation of stripe order phenomena, and also motivates us to undertake a more refined investigation of how interlayer interactions and on-site interactions influence stripe order.

In this work, we use the PQMC algorithm to compute the ground-state properties of a sign-problem-free bilayer Hubbard-like model, featured by the interlayer antiferromagnetic (AFM) spin-exchange anisotropy. We discover the direct evidence of competition between stripes and superconductivity, which is sensitive to the anisotropic. We find that the charge stripe phase is present when the interlayer spin-exchange is highly anisotropic. Moreover, we find that this stripe phase is unstable upon introducing the spin-flip spin-exchange term J_{\mathrm{\perp }}, and that an interlayer pairing superconductivity enhances when the stripe phase is suppressed. Due to the flexibility of our model, we are able to tune the value of U\in [\frac{J_z}{4}-\frac{J_{\mathrm{\perp }}}{2},\frac{J_z}{4}+\frac{J_{\mathrm{\perp }}}{2}], without encountering sign problems. At small doping, the on-site interaction U is detrimental to both the stripe phase and superconductivity. While the superconductivity is enhanced by U at large doping, where charge stripe is absent. When charge stripes exist, we also discover the spin stripes, characterized by twice the period of charge stripes, indicating the profound internal connection between them. We discuss the possible origin of these phenomena and their connections to the standard Hubbard model.

2 Model and method

Our model was first proposed in Ref. [39], and subsequent works have found exotic phases and superconductivity within it [36, 40–43]. This model originated from the \mathrm{S}\mathrm{O}(5) symmetric model [44], and can be generalized to a bilayer model that contains both on-site interaction and interlayer spin-exchange. The Hamiltonian can be written as follows:

(1)\begin{array}{rl}H=& -t\sum _{⟨i,j⟩l\sigma }\left(c_{il\sigma }^{†}{c}_{jl\sigma }+\mathrm{H}.\mathrm{c}.\right)+U\sum _{il}{n}_{il\uparrow }{n}_{il\downarrow }\\ & +J_z\sum _i{S}_{i1}^z{S}_{i2}^z+\frac{J_{\mathrm{\perp }}}{2}\sum _i\left(S_{i1}^{+}{S}_{i2}^{-}+\mathrm{H}.\mathrm{c}.\right),\end{array}

where c_{il\sigma }(c_{il\sigma }^{†}) denotes the annihilation (creation) operator at site i, layer l, and spin \sigma. n_{il\sigma }=c_{il\sigma }^{†}{c}_{il\sigma } is density operator, S_{il}^z=\frac{n_{il\uparrow }-n_{il\downarrow }}{2} is z component spin operator, and S_{il}^{+}=c_{il\uparrow }^{†}{c}_{il\downarrow } (its Hermitian conjugate S_{il}^{-}) is spin annihilation (creation) operator, contributing to the xy component spin. The non-interaction term only includes the intralayer nearest-neighbor hopping t, and the interaction terms involve the intralayer on-site repulsion U, the z-component term J_z and the spin-flip (xy-component) term J_{\perp } of interlayer AFM spin-exchange.

The absence of the sign problem in this model originates from time-reversal symmetry. To illustrate this, we define {\psi }_i={\left(c_{i1\uparrow },c_{i1\downarrow },c_{i2\uparrow },c_{i2\downarrow }\right)}^{\mathrm{T}}, and rewrite the Hamiltonian in the following form:

(2)\begin{array}{rl}H=& -t\sum _{⟨i,j⟩}\left({\psi }_i^{†}{\psi }_j+H.c.\right)-\sum _{\alpha =1}^5\frac{g_{\alpha }}{2}{\left({\psi }_i^{†}{\mathrm{\Gamma }}^{\alpha }{\psi }_i\right)}^2\\ & -\frac{g_0}{2}\sum _i{\left({\psi }_i^{†}{\psi }_j-2\right)}^2,\end{array}

where we define

(3)\begin{array}{rl}& {\mathrm{\Gamma }}^{1-3}=\left(\begin{array}{cc}\sigma & 0\\ 0& -\sigma \end{array}\right),{\mathrm{\Gamma }}^4=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right),\\ & {\mathrm{\Gamma }}^5=\left(\begin{array}{cc}0& \mathrm{i}I\\ -\mathrm{i}I& 0\end{array}\right),\end{array}

and \sigma represents the Pauli matrices. We focus on those five {\mathrm{\Gamma }}^{\alpha } matrices which are proven to be invariant under the time-reversal transformation \mathcal{T}=R\mathbb{C}, where \mathbb{C} denotes complex conjugation and R is given by

(4)R=\left(\begin{array}{cc}0& -\mathrm{i}{\sigma }_2\\ -\mathrm{i}{\sigma }_2& 0\end{array}\right).

This ensures that the fermion determinant is positive. Furthermore, the mapping from Eq. (1) to Eq. (2) is performed by the parameter relation:

(5)\begin{array}{rl}& g_{1,2}=\frac{1}{4}U+\frac{1}{8}{J}_{\mathrm{\perp }}-\frac{1}{16}{J}_z,\\ & g_3=\frac{1}{4}U-\frac{1}{8}{J}_{\mathrm{\perp }}+\frac{3}{16}{J}_z,\\ & g_{0,4,5}=-\frac{1}{4}U+\frac{1}{8}{J}_{\mathrm{\perp }}+\frac{1}{16}{J}_z,\end{array}

where the sign problems are avoided by g_{1-5}\ge 0. For these details, see Ref. [39] and Ref. [41].

In this study, we chose a small enough time-step \mathrm{\Delta }\tau =0.1 so that the trotter errors [O(\mathrm{\Delta }{\tau }^2)] are smaller than those associated with the statistical sampling. Furthermore, we chose a large enough projection time parameter \mathrm{\Theta }=24-28 for different sizes to converge the system to the ground state. Specifically, with L_y=8 fixed, we chose \mathrm{\Theta }=24 for L_x\le 16, \mathrm{\Theta }=25 for L_x=24 and \mathrm{\Theta }=28 for L_x=32. These parameters are sufficient to converge the system, and an evidence of convergence is provided in Appendix A.

3 Results and discussion

When the interlayer interactions are purely anisotropic, previous works have shown that this model exhibits both charge and spin stripe phases [37, 38]. Moreover, this state remains robust even when next-nearest-neighbor hopping is incorporated. Those works only consider the z-component ferromagnetic spin-exchange interaction between layers due to the limitation of sign-problem, leaving the effect of AFM interlayer spin-exchange interaction and its anisotropy unclear. To reveal these effects, we define the charge correlation function as follows:

(6)S_c\left(\mathit{k}\right)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i,j}{\mathrm{e}}^{-\mathrm{i}\mathit{k}\cdot \left({\mathit{R}}_i-{\mathit{R}}_j\right)}\left[⟨n_i{n}_j⟩-⟨n_i⟩⟨n_j⟩\right],}

where N is the number of unit cells and n_i={\sum }_{l\sigma }{n}_{il\sigma } for convenience of expression. In the following, we present the charge correlation function at k_y=0, where the peak position along k_x reflects the periodicity of charge modulations. Unless otherwise specified, we use the system size L_y=8 and L_x=2/\delta with periodic boundary conditions, where \delta denotes the hole doping density, and we set t=1 as the energy scale.

As shown in Fig. 1(a), when the system contains only the U and J_z terms, the peak in S_c(k_x,0) is at k_x=\pm 2\text{π}\delta. The peak intensity of S_c is maximal at 1/8 doping, a feature consistent with previous work on cuprate superconductors [17, 45], and then suppressed by the further increasing doping. This phenomenon is interesting because the sign-problem-free model we employ shows stripe behavior similar to that in cuprate superconductors, despite the different geometries of the two systems. This suggests that the underlying mechanisms might be similar. Previous work (Ref. [38]) also reveals that this model shares another remarkable similarity with the single-layer Hubbard model, and the introduction of a next-nearest-neighbor hopping term shortens the stripe wavelength. Focusing on the topic around effects of interlayer AFM spin-exchange and its anisotropy, we summarize our conclusions as a rough strip-shaped phase diagram at \delta =1/8 or 1/6 in Fig. 1(b). The stripe phase is robust at the high anisotropy of spin-exchange region, while gradually vanishes and transforms into the superconducting phase as the anisotropy decreases (i.e., as J_{\perp }/J_z increases). It reflects a key competition between stripe and superconductivity. Besides, with J_{\perp }/J_z increasing, the charge stripe may become the finite-size artifact or directly vanish. These results will be detailedly discussed in the following.

We now take a deeper look at the stripe phase. Since the interlayer Ising term J_z favors antiparallel spin pairs on a rung, holes introduced into the system do not spread homogeneously because doing so would break these spin pairs. Holes instead concentrate on specific sites to minimize the system’s energy. However, this order is readily disrupted by introducing fluctuations to the interlayer spin-exchange; namely, adding the J_{\mathrm{\perp }} term suppresses the stripe phase. To maintain this model as sign-problem-free, the relation J_z+2J_{\mathrm{\perp }}=4U must be satisfied. With U=2.0 fixed, we tune the ratio of J_z to J_{\mathrm{\perp }} and observe that this peak disappears rapidly, as shown in Fig. 2.

When we introduce J_{\mathrm{\perp }}=0.6 at \delta =1/8, we observe that the peak position shifts, as shown in Fig. 2. The peak moves to \text{π}/8, which corresponds to a charge period that matches the system size L_x given the condition L_x=2/\delta. Whether this period shift is real or a finite-size effect in this system is a critical question, which means that the larger L_x should be considered.

To verify the authenticity of the stripes, we perform the finite-size analysis. With U=2.0 fixed, we further introduce J_{\perp } to change the anisotropy of this interlayer spin-exchange interaction, as shown in Fig. 3. In the pure J_z case, the charge correlation is strengthened as L_x increases at the constant peak momentum k_x=\pm 2\text{π}\delta for both (a) \delta =1/8 and (d) \delta =1/6, which indicates the real and robust (or long-range) charge stripes. When increasing J_{\perp } to 0.4 in Figs. 3(b) and (e), the position of the peak momentum still persists. However, the peak for \delta =1/8 declines from L_x=16 to L_x=32, indicating real but weak (or short-range) charge stripes. Further increasing J_{\perp } to 0.6, the charge stripes are absent, which have different performances for \delta =1/8 and \delta =1/6, as shown in Figs. 3(c) and (f). For \delta =1/8, the position of momentum peak shifts to minimum non-zero momentum as mentioned in Fig. 2(a) and it changes with L_x different, demonstrating a finite-size artifact. We emphasize that this peak momentum shift with a fixed size is a strong signal of finite-size effect and the stripe has vanised. While for \delta =1/6, the peak directly vanishes. The finite-size analysis confirms the real charge stripes in our system, and further elucidates that charge stripes are suppressed and destroyed by decreasing anisotropy of interlayer AFM spin-exchange and how it evolves.

In the single-layer Hubbard model, a comparable spin-exchange interaction arises from a perturbative treatment of electron hopping, leading to an effective t-J model [46, 47]. The stripe state in the single-layer model maintains the energetically favored AFM order between stripes. In this sense, the mechanism underlying stripe formation in the two models is similar. We have confirmed this mechanism in our system (See Appendix B for details). Spin-exchange in the single-layer lattice also involves attraction between antiparallel spins. The presence of spin-exchange between four nearest-neighboring sites can induce a larger energy difference. Consequently, the stripe phase in the single-layer lattice is expected to be more stable than that in the bilayer lattice, where only interlayer spin-exchange exists. Furthermore, Ref. [48] studied nonlocal interactions in the extended Hubbard model and found that nearest-neighbor attraction enhances AFM correlations. To conclude this section, the stripe phase in the two-dimensional Hubbard model is likely more stable, which reflects the difference between our sign-free model and the pure Hubbard model.

Next, we gradually suppress the U term in the Hamiltonian under the condition of J_{\mathrm{\perp }}\ne 0. To avoid the sign problem, the value of U must satisfy \frac{J_z}{4}-\frac{J_{\mathrm{\perp }}}{2}\le U\le \frac{J_z}{4}+\frac{J_{\mathrm{\perp }}}{2}. Given that the on-site interaction U prohibits double occupation, it tends to favor the formation of antiparallel spins. However, U also favors other ordered phases, including superconductivity and in-plane magnetism. The involvement of possible intralayer order constitutes another difference between this model and the single-layer Hubbard model. To reveal the effect of U in our system, we measure the charge correlation with sufficiently large J_{\perp } to avoid the sign problem, which may result in the finite-effect artifact of stripe as mentioned above. Therefore we exhibit the U-effect at both real stripe and its finite-size artifact regions as shown in Figs. 4(a) and (b) respectively. In Fig. 4(a) with J_z=7.2 and J_{\perp }=0.4, increasing U suppresses the stripe phase for both \delta =1/8 and \delta =1/6, which suggests that U may be more inclined towards intralayer order. However in Fig. 4(b), the position of peak for \delta =1/8 under J_z=11.0 and J_{\perp }=0.5 shifts to the minimum momentum, indicating a finite-size artifact, which is detailedly discussed in Appendix C. Although the stripe becomes the finite-size artifact at \delta =1/8 or vanishes at \delta =1/6, the suppression of U still persists.

The superconductivity in this system also shows a complex dependence on doping and interaction strength. We focus on the interlayer pairing and its correlation functions which can be defined as

(7)S_{sc}(\mathit{k})={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i,j}{\mathrm{e}}^{-\mathrm{i}\mathit{k}\cdot \left({\mathit{R}}_i-{\mathit{R}}_j\right)}⟨{\mathrm{\Delta }}_i^{†}{\mathrm{\Delta }}_j⟩,}

where {\mathrm{\Delta }}_i=c_{i1\uparrow }{c}_{i2\downarrow }-c_{i1\downarrow }{c}_{i2\uparrow }, and the uniform pairing correlation function corresponds to S_{sc}(\mathit{k}=0), abbreviated as S_{sc}.

Now we present the pairing correlation function as a function of electron filling ⟨n⟩=1-\delta in Fig. 5. This figure clearly reveals the evidence of competition between charge stripe and superconductivity. On the one hand, the isotropy of interlayer AFM spin-exchange (i.e., introducing J_{\mathrm{\perp }}) enhances superconductivity by disrupting the charge stripe. On the other hand, for the pure J_z case, superconductivity decreases sharply around 1/4 doping where charge stripe emerges or strengthens for both U = 2 and U = 3. This abrupt change is particularly pronounced for U=3. Thus, for the pure J_z case, the pairing correlation function exhibits a flat profile below 1/4 doping, attributed to suppression by charge stripes. In contrast, in the presence of J_{\mathrm{\perp }}, the pairing correlation function shows a smooth increase as the doping level rises to 1/4, since charge stripes are weakened and gradually disappear in this region. This provides clear evidence for the competition between charge stripes and superconductivity. In addition, a finite-size analysis of pairing correlation function for U=2 case demonstrates that long-range pairing correlation is absent when the charge stripe exists and it emerges when the stripe have vanished, further strongly supporting this competition, and see the detail in Appendix D.

However, the relationship between superconductivity and charge stripes in this model is complex. We gradually decrease U when J_{\mathrm{\perp }} is present, as shown in Fig. 6. Due to the constraints on the model parameters, a smaller U can be achievable when J_{\mathrm{\perp }} is large. We previously found that the charge stripes are suppressed as U increases, and we now observe that superconductivity is also suppressed at low doping levels with increasing U, as shown in Fig. 6(b). This provides another piece of evidence for the influence of intralayer ordered phases. At lower electron densities, where charge stripes do not form, increasing U enhances superconductivity. This reflects the multiple roles played by parameters U and J, while the electron density also plays a significant role in regulating the system’s properties.

The charge stripes are often accompanied by some form of spin stripes. In the single layer Hubbard model [13, 14], the period of the spin stripe is always twice that of the charge stripe. There is an internal connection between charge and spin stripes, for example, the AFM order is favored between charge stripes. Furthermore, we also explore the spin stripe in our system. We define the spin correlation function as

(8)S_s\left(\mathit{k}\right)={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i,j}{\mathrm{e}}^{-\mathrm{i}\mathit{k}\cdot \left({\mathit{R}}_i-{\mathit{R}}_j\right)}⟨O_i^s{O}_j^s⟩,}

with O_i^s=\left(c_{i1\uparrow }^{†}{c}_{i1\uparrow }-c_{i1\downarrow }^{†}{c}_{i1\downarrow }-c_{i2\uparrow }^{†}{c}_{i2\uparrow }+c_{i2\downarrow }^{†}{c}_{i2\downarrow }\right) characterized by the interlayer AFM alignment.

The spin correlation function is calculated in the same parameters, as shown in Fig. 7. We clearly observe the spin stripe within the background of AFM order at the momentum (\text{π}\pm \text{π}\delta ,\text{π}), indicating the twice period of charge stripe, and further finite-size analysis (i.e., the peak momentum remains as L_x increases) demonstrates its authenticity. In the pure J_z case, there are maximal peaks for both (a) \delta =1/8 and (d) \delta =1/6, and the peak grows as L_x increases, implying the presence of robust spin stripes. There are still similar behaviors after introducing J_{\perp }=0.4 in Figs. 7(b) and (e), but the peaks are suppressed by J_{\perp }. While for J_{\perp }=0.6, the position of peak shifts to the momentum closest to (\text{π},\text{π}) as L_x changes as shown in Figs. 7(c) and (f), indicating the collapse of spin stripes in the form of finite-size artifacts. It is very consistent and synchronous with the disappearance of the charge stripes, implying the important internal connection between spin and charge stripes. The detailed connection can be explored in the real space, and see the Appendix B. Besides, in all cases, the spins along the L_y direction always maintain AFM arrangement, characterized by the maximum at (\text{π},\text{π}), as shown in the insets of Fig. 7.

4 Conclusions

Using the PQMC algorithm, we study a sign-problem-free bilayer Hubbard-like model to investigate the competitive interplay between charge stripes and interlayer superconductivity, which is controlled by interlayer spin-exchange anisotropy (J_z and J_{\mathrm{\perp }}). Due to parameter constraint J_z+2J_{\perp }=4U which avoids the sign problem, the accessible parameter space is so limited that we can not study ideal parameters such as fixing J_z and introducing J_{\perp }. However we mainly focus on the study on the interlayer spin-exchange anisotropy. On the one hand, we gradually increase J_{\perp } and decrease J_z by small intervals (like \mathrm{\Delta }{J}_{\perp }=0.2), which promotes the reliability of our conclusions. On the other hand, the anisotropy can be effectively quantified as the ratio of J_{\perp } and J_z (J_{\perp }/J_z), which further weakens the dependence on fixing J_z or J_{\perp }. Hence our following conclusions are still reliable.

A charge stripe phase, analogous to that in the pure two-dimensional Hubbard model, emerges in the bilayer model with anisotropic interlayer AFM spin-exchange interaction. We find that the interlayer spin-exchange anisotropy plays a crucial role in the formation of both charge stripes and interlayer superconductivity. The charge stripe phase, which is characterized by a peak at (0\pm 2\text{π}\delta ,0), is stable for highly anisotropic spin-exchange (especially for pure J_z case) but strongly suppressed by decreasing the anisotropy (i.e., introducing the spin-flip term J_{\perp }). Concurrently, superconducting correlations are enhanced and become robust as charge stripes diminish, indicating a direct competition between these two phases. We also discuss the possible relationship between this model and the pure two-dimensional Hubbard model, in which such competition also exists. Besides, the on-site interaction U plays a complex role: As U increases, it suppresses both charge stripes and superconductivity under low doping, suggesting the emergence of other intralayer orders. However, with further enhancing doping, increasing U enhances superconductivity in the absence of charge stripes. Furthermore, we also discover the synchronous spin stripes with twice the period of the charge stripes, indicating a key internal connection between the two. This work identifies the key factors contributing to the formation of charge stripes, and highlights the sensitivity of stripe and superconducting phases to interaction parameters (especially interlayer AFM spin-exchange anisotropy), thus offering insights into competing orders in strongly correlated systems.

5 Appendix A PQMC algorithm and convergence

In this work, we use the PQMC method, through which the observables are calculated by the projection as follows:

(A.1)\begin{array}{c}⟨O⟩={\displaystyle \frac{⟨{\mathrm{\Phi }}_T|{\mathrm{e}}^{-\mathrm{\Theta }H}O{\mathrm{e}}^{-\mathrm{\Theta }H}|{\mathrm{\Phi }}_T⟩}{⟨{\mathrm{\Phi }}_T|{\mathrm{e}}^{-2\mathrm{\Theta }H}|{\mathrm{\Phi }}_T⟩}}.\end{array}

To treat the quantum system, we first utilize the Trotter decomposition, {\mathrm{e}}^{-\mathrm{\Theta }H}={\left({\mathrm{e}}^{-\mathrm{\Delta }\tau H}\right)}^{N_t}=\prod _{\tau }^{N_t}{\mathrm{e}}^{-\mathrm{\Delta }\tau V}{\mathrm{e}}^{-\mathrm{\Delta }\tau T}, where N_t\mathrm{\Delta }\tau =\mathrm{\Theta }, T is the noninteracting operator [i.e., hopping term -t{\sum }_{⟨i,j⟩}({\psi }_i^{†}{\psi }_j+H.c.)], and V represents the interacting operator (i.e., the remaining interacting terms in Hamiltonian Eq. (1)). In our Hamiltonian, each term in V is a perfect square of bilinear operator. We decouple these terms using the standard HS transformation:

(A.2)\begin{array}{c}\begin{array}{rl}{\mathrm{e}}^{-\mathrm{\Delta }\tau V}=& {\mathrm{e}}^{\sum _{\alpha =0}^5{\sum }_i\frac{\mathrm{\Delta }\tau g_{\alpha }}{2}{\left({\psi }_i^{†}{\mathrm{\Gamma }}^{\alpha }{\psi }_i-z_{\alpha }\right)}^2}\\ =& \prod _{\alpha ,i}\sum _{k=\pm 1,\pm 2}{\gamma }_{k,i}{\mathrm{e}}^{\sqrt{\mathrm{\Delta }\tau g_{\alpha }}{\eta }_{k,i}\left({\psi }_i^{†}{\mathrm{\Gamma }}^{\alpha }{\psi }_i-z_{\alpha }\right)}\\ =& \sum _k\prod _i{\gamma }_{k,i}U(k,i),\end{array}\end{array}

where z_0=2 and z_{1-5}=0. The coefficients \gamma and \eta follow the standard form for the four-component HS transformation, with details available in Ref. [49]. After the HS decomposition, the system transforms into an ensemble of free fermions. The Green’s functions for a given configuration can be calculated easily, thereby rendering the observables a weighted average over these free-fermion configurations. The weight of each configuration is defined as \omega =\prod {\gamma }_k^{\tau }det[P^{†}{B}_V({\mathit{k}}_{N_t})B_T\dotsB_V({\mathit{k}}_1)B_{T}P], where B_T denotes the matrix form of {\mathrm{e}}^{-\mathrm{\Delta }\tau T} in the single-particle basis, and B_V({\mathit{s}}_{\tau }) is the matrix form of \prod _{i}U(k,i) for a specific configuration {\overrightarrow{s}}_{\tau } at imaginary time \tau in the single-particle basis. P is a matrix that characterizes the Slater determinant of the trial wave function {\mathrm{\Phi }}_T. The weight \omega is crucial for PQMC simulations, and the sign problem emerges when its value is non-positive. Fortunately, if the matrix in the determinant remains invariant under the \mathcal{T} transformation, its determinant is consistently positive. The model employed in this work satisfies this condition, allowing us to perform quantum Monte Carlo simulations free of the sign problem. Moreover, all unique combinations of parameter values in the PQMC simulations are listed in Table A1 for the clarity.

In the PQMC simulating process, we also measure the convergence of projection time \mathrm{\Theta } to ensure that the system can be converged to the ground state. We provide an example for the charge correlation function at \delta =1/8 case with different interaction parameters and sizes, as shown in Fig. A1. Regardless of L_x=16 or L_x=32, or whether J_{\perp } is introduced or not, S_c deviates greatly when \mathrm{\Theta }=16, However when \mathrm{\Theta }\ge 20, all S_c remain basically consistent within the error range. Therefore \mathrm{\Theta } we considered in the text is sufficient to converge the system.

6 Appendix B The antiferromagnetic order between charge stripes

To verify the mechanism about AFM order between stripes as described in text, we define the spin correlation function in real space as

(B.1)\begin{array}{c}{C}_s(\mathit{R})={\displaystyle \frac{1}{N}}{\displaystyle \sum _{{\mathit{R}}_i-{\mathit{R}}_j=\mathit{R}}⟨O_i^s{O}_j^s⟩,}\end{array}

where O_i^s is the same definition as in Eq. (8). We plot the charge and spin correlation function in real space at J_z=8, J_{\perp }=0 and \delta =1/6 as an example, as shown in Fig. A2. The spin correlation is processed by the staggered factor [i.e., (-1{)}^{x+y}{C}_s(x,y)] to highlight the spin stripe and the boundary of \text{π}-phase shift. Although the pattern of a perfect period in real space is limited by the strict symmetry about x=0 and lattice sizes (instead its period is more accurate in momentum space), Fig. A2 is still impressive, especially in the corresponding relationship between charge and spin stripes.

The charge stripe in Fig. A2(a) is characterized by high electron (blue circle area) and high hole (red circle area) concentrations, while the separation of red and blue circles indicates the \pi-phase shift in the context of AFM order in Fig. A2(b). By comparison, the \pi-phase shift always emerges in high hole-concentration region, which is also the reason why there are always weak spin correlations (small circles) on the boundary of \pi-phase shift. It results in a favored AFM order (\pi-phase shift) between neighbor charge stripes. In turn, it requires that the period of the spin stripe is twice that of the charge stripe, and we have confirmed this viewpoint in the momentum space in text. These results are consistent with previous works [13, 14, 38], emphasizing the important internal connections between change and spin stripes.

7 Appendix C The J_z=12 case and its anisotropy

We also measure the larger interlayer spin-exchange interaction cases, namely J_z+2J_{\perp }=4U with fixed U=3. As shown in Fig. A3, there are similar behaviors with the U=2 case. Specifically, for both (a) \delta =1/8 and (b) \delta =1/6, the significant stripe peaks stand out with the consistent momentum K_x=2\pi \delta at pure J_z=12, while the introduction and increase of J_{\perp } shifts the peak momentum to a minimum non-zero value for \delta =1/8, indicating the emergence of finite-size artifacts just like Fig. 2 case. For \delta =1/6, the peak is suppressed by J_{\perp } with almost constant peak momentum. Regardless of the case of U=2 or U=3, we discover that the charge stripes for different doping almost simultaneously disappear, either becoming a finite-size product or being directly destroyed.

8 Appendix D The finite-effect analysis of pairing correlations

In the text, we discover the pairing correlation functions have the competitive behavior with stripes, as discussed in Fig. 5. However whether the superconducting correlation is long-ranged is unknown, which can be a more profound signal of this competition. To identify the existence of long-ranged interlayer pairing superconductivity in our system, we resort to the finite-size analysis. The finite-size results further confirm the competition between stripe and superconductivity, as shown in Fig. A4.

For the pure J_z cases, as can be seen in Fig. A4, the uniform pairing correlation at (0,0) is impaired as L increases at both (a) ⟨n⟩=0.875 (\delta =1/8) and (b) ⟨n⟩=0.833 (\delta =1/6), indicating the absence of long-range pairing correlations, where the stripes are instead prominent (or long-ranged). With doping increasing, the stripes are weakened as indicated in Fig. 1, and the superconducting pairing correlation becomes long-ranged at (c) ⟨n⟩=0.667 (\delta =1/3), characterized by an increasing trend with L. It is completely consistent with the sharp decline around 1/4 doping for the pure J_z cases of Fig. 5. After introducing J_{\perp }=1, the superconducting pairing correlations always are long-ranged at any electron filling in Figs. A4(d)−(f), while the stripes have been destroyed by the decreasing anisotropy of interlayer spin-exchange. These signals strongly support the competition between stripe and superconductivity.

References

[1] D. J. Scalapino , A common thread: The pairing interaction for unconventional superconductors, Rev. Mod. Phys. 84(4), 1383 (2012)

[2] S. Raghu , S. A. Kivelson , and D. J. Scalapino , Superconductivity in the repulsive Hubbard model: An asymptotically exact weak-coupling solution, Phys. Rev. B 81(22), 224505 (2010)

[3] E. Dagotto , Correlated electrons in high-temperature superconductors, Rev. Mod. Phys. 66(3), 763 (1994)

[4] J. E. Hirsch and H. Q. Lin , Pairing in the two-dimensional Hubbard model: A Monte Carlo study, Phys. Rev. B 37(10), 5070 (1988)

[5] H. Q. Lin , J. E. Hirsch , and D. J. Scalapino , Pairing in the two-dimensional Hubbard model: An exact diagonalization study, Phys. Rev. B 37(13), 7359 (1988)

[6] T. Maier , M. Jarrell , T. Pruschke , and J. Keller , d-wave superconductivity in the Hubbard model, Phys. Rev. Lett. 85(7), 1524 (2000)

[7] D. P. Arovas , E. Berg , S. A. Kivelson , and S. Raghu , The Hubbard model, Annu. Rev. Condens. Matter Phys. 13(1), 239 (2022)

[8] E. H. da Silva Neto , P. Aynajian , A. Frano , R. Comin , E. Schierle , E. Weschke , A. Gyenis , J. Wen , J. Schneeloch , Z. Xu , S. Ono , G. Gu , M. Le Tacon , and A. Yazdani , Ubiquitous interplay between charge ordering and high-temperature superconductivity in cuprates, Science 343(6169), 393 (2014)

[9] J. M. Tranquada , B. J. Sternlieb , J. D. Axe , Y. Nakamura , and S. Uchida , Evidence for stripe correlations of spins and holes in copper oxide superconductors, Nature 375(6532), 561 (1995)

[10] J. Chang , E. Blackburn , A. T. Holmes , N. B. Christensen , J. Larsen , J. Mesot , R. Liang , D. A. Bonn , W. N. Hardy , A. Watenphul , M. Zimmermann , E. M. Forgan , and S. M. Hayden , Direct observation of competition between superconductivity and charge density wave order Âain YBa₂Cu₃O_6.67, Nat. Phys. 8(12), 871 (2012)

[11] R. Comin and A. Damascelli , Resonant X-ray scattering studies of charge order in cuprates, Annu. Rev. Condens. Matter Phys. 7(1), 369 (2016)

[12] M. Qin , C. M. Chung , H. Shi , E. Vitali , C. Hubig , U. Schollwöck , S. R. White , and S. Zhang (Simons Collaboration on the Many-Electron Problem) , Absence of superconductivity in the pure two-dimensional Hubbard model, Phys. Rev. X 10(3), 031016 (2020)

[13] B. X. Zheng , C. M. Chung , P. Corboz , G. Ehlers , M. P. Qin , R. M. Noack , H. Shi , S. R. White , S. Zhang , and G. K. L. Chan , Stripe order in the underdoped region of the two-dimensional Hubbard model, Science 358(6367), 1155 (2017)

[14] H. Xu , C. M. Chung , M. Qin , U. Schollwöck , S. R. White , and S. Zhang , Coexistence of superconductivity with partially filled stripes in the Hubbard model, Science 384(6696), eadh7691 (2024)

[15] C. Zhang,J.W. Li,J. von Delft, Frustration-induced superconductivity in the t-t' Hubbard model, arXiv: 2023)

[16] H. C. Jiang and T. P. Devereaux , Superconductivity in the doped Hubbard model and its interplay with next-nearest hopping t', Science 365(6460), 1424 (2019)

[17] E. W. Huang , C. B. Mendl , H.-C. Jiang , B. Moritz , and T. P. Devereaux , Stripe order from the perspective of the Hubbard model, npj Quantum Mater. 3, 22 (2018)

[18] Z. Chen , Y. Wang , S. N. Rebec , T. Jia , M. Hashimoto , D. Lu , B. Moritz , R. G. Moore , T. P. Devereaux , and Z. X. Shen , Anomalously strong near-neighbor attraction in doped 1d cuprate chains, Science 373(6560), 1235 (2021)

[19] Z. Zhou , W. Ye , H. G. Luo , J. Zhao , and J. Chang , Robust superconducting correlation against intersite interactions in the extended two-leg Hubbard ladder, Phys. Rev. B 108(19), 195136 (2023)

[20] H. C. Jiang , Pair density wave in the doped three-band Hubbard model on two-leg square cylinders, Phys. Rev. B 107(21), 214504 (2023)

[21] C. Peng , Y. Wang , J. Wen , Y. S. Lee , T. P. Devereaux , and H. C. Jiang , Enhanced superconductivity by near-neighbor attraction in the doped extended Hubbard model, Phys. Rev. B 107(20), L201102 (2023)

[22] L. Zhang , T. Guo , Y. Mou , Q. Chen , and T. Ma , Enhancement of d-wave pairing in the striped phase with nearest neighbor attraction, Phys. Rev. B 105(15), 155154 (2022)

[23] S. Karakuzu , A. Tanjaroon Ly , P. Mai , J. Neuhaus , T. A. Maier , and S. Johnston , Stripe correlations in the two-dimensional Hubbard−Holstein model, Commun. Phys. 5(1), 311 (2022)

[24] C. Chen , P. Zhong , X. Sui , R. Ma , Y. Liang , S. Hu , T. Ma , H. Q. Lin , and B. Huang , Charge stripe manipulation of superconducting pairing symmetry transition, Nat. Commun. 15(1), 9502 (2024)

[25] C. Chen , Z. Fan , R. Ma , Y. Pan , Y. Liang , B. Huang , and T. Ma , Stripe order manipulated dominant pairing symmetry in the Hubbard model, Phys. Rev. B 109(4), 045101 (2024)

[26] C. Cheng , R. Mondaini , and M. Rigol , Singlet pairing and superconductivity in t−j ladders with mott insulating stripes, Phys. Rev. B 98(12), 121112 (2018)

[27] R. Mondaini , T. Ying , T. Paiva , and R. T. Scalettar , Determinant quantum Monte Carlo study of the enhancement of d-wave pairing by charge inhomogeneity, Phys. Rev. B 86(18), 184506 (2012)

[28] J. P. F. LeBlanc , A. E. Antipov , F. Becca , I. W. Bulik , G. K. L. Chan , C. M. Chung , Y. Deng , M. Ferrero , T. M. Henderson , C. A. Jiménez-Hoyos , E. Kozik , X. W. Liu , A. J. Millis , N. V. Prokofev , M. Qin , G. E. Scuseria , H. Shi , B. V. Svistunov , L. F. Tocchio , I. S. Tupitsyn , S. R. White , S. Zhang , B. X. Zheng , Z. Zhu , and E. Gull (Simons Collaboration on the Many-Electron Problem) , Solutions of the two-dimensional Hubbard model: Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5(4), 041041 (2015)

[29] G. Knizia and G. K. L. Chan , Density matrix embedding: A simple alternative to dynamical mean-field theory, Phys. Rev. Lett. 109(18), 186404 (2012)

[30] N. V. Prokofév and B. V. Svistunov , Polaron problem by diagrammatic quantum Monte Carlo, Phys. Rev. Lett. 81(12), 2514 (1998)

[31] S. R. White , Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69(19), 2863 (1992)

[32] G. Kotliar , S. Y. Savrasov , K. Haule , V. S. Oudovenko , O. Parcollet , and C. A. Marianetti , Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78(3), 865 (2006)

[33] T. Maier , M. Jarrell , T. Pruschke , and M. H. Hettler , Quantum cluster theories, Rev. Mod. Phys. 77(3), 1027 (2005)

[34] T. Guo , L. Zhang , T. Ma , and H. Q. Lin , Dual role of the stripe phase in superconducting correlations in a bilayer square lattice, Phys. Rev. B 111(19), 195127 (2025)

[35] R. Blankenbecler , D. J. Scalapino , and R. L. Sugar , Monte Carlo calculations of coupled boson-fermion systems. i, Phys. Rev. D 24(8), 2278 (1981)

[36] T. Ma , D. Wang , and C. Wu , Doping-driven antiferromagnetic insulator−superconductor transition: A quantum Monte Carlo study, Phys. Rev. B 106(5), 054510 (2022)

[37] F. F. Assaad , V. Rousseau , F. Hebert , M. Feldbacher , and G. G. Batrouni , Spin and charge dynamics of stripes in doped mott insulators, Europhys. Lett. 63(4), 569 (2003)

[38] Y. Caplan and D. Orgad , Quantum Monte Carlo study of a bilayer u(2)×u(2)-symmetric Hubbard model, Phys. Rev. B 108(16), 165131 (2023)

[39] C. Wu , J. P. Hu , and S. C. Zhang , Exact SO(5) symmetry in the spin-3/2 fermionic system, Phys. Rev. Lett. 91(18), 186402 (2003)

[40] S. Capponi , C. Wu , and S. C. Zhang , Current carrying ground state in a bilayer model of strongly correlated systems, Phys. Rev. B 70(22), 220505 (2004)

[41] C. Wu and S. C. Zhang , Sufficient condition for absence of the sign problem in the fermionic quantum Monte Carlo algorithm, Phys. Rev. B 71(15), 155115 (2005)

[42] R. Ma and T. Ma , Competition between antiferromagnetism and superconductivity in a doped Hubbard model with anisotropic interaction, Phys. Rev. B 107(21), 214509 (2023)

[43] R. Ma , Z. Fan , T. Ma , and C. Wu , Parameters dependent superconducting transition temperature in a sign-problem-free bilayer model, Chin. Phys. Lett. 42(11), 110705 (2025)

[44] D. Scalapino , S. C. Zhang , and W. Hanke , SO(5) symmetric ladder, Phys. Rev. B 58(1), 443 (1998)

[45] Y. F. Jiang , J. Zaanen , T. P. Devereaux , and H. C. Jiang , Ground state phase diagram of the doped Hubbard model on the four-leg cylinder, Phys. Rev. Res. 2(3), 033073 (2020)

[46] K. A. Chao , J. Spałek , and A. M. Oles , Canonical perturbation expansion of the Hubbard model, Phys. Rev. B 18(7), 3453 (1978)

[47] Y. Tian, Mapping the Hubbard model to the t−j model using ground state unitary transformations, arXiv: 2019)

[48] M. Yao , D. Wang , and Q. H. Wang , Determinant quantum Monte Carlo for the half-filled Hubbard model with nonlocal density−density interactions, Phys. Rev. B 106(19), 195121 (2022)

[49] F. Goth, Higher order auxiliary field quantum Monte Carlo methods, arXiv: 2020)

RIGHTS & PERMISSIONS

Higher Education Press