Variational Monte Carlo study of stripes as a function of doping in the t − t ′ Hubbard model

. We perform variational Monte Carlo simulations of the single-band Hubbard model on the square lattice with both nearest ( t ) and next-nearest ( t ′ ) neighbor hoppings. Our work investigates the consequences of increasing hole doping on the instauration of stripes and the behavior of the superconducting order parameter, with a discussion on how the two phenomena affect each other. We consider two different values of the next-nearest neighbor hopping parameter, that are appropriate for describing cuprate superconductors. We observe that stripes are the optimal state in a wide doping range; the stripe wavelength reduces at increasing doping, until stripes melt into a uniform state for large values of doping. Superconducting pair-pair correlations, indicating the presence of superconductivity, are always suppressed in the presence of stripes. Our results suggest that the phase diagram for the single-band Hubbard model is dominated by stripes, with superconductivity being possible only in a narrow doping range between striped states and a nonsuperconducting metal.


Introduction
The concept of a stripe phase is one of the unconventional features that emerged over the years when interpreting a broad range of experimental results on copper oxide superconductors.Once we allow holes to wander in an antiferromagnetic background, the creation of striped inhomogeneities can be a consequence.The reason for the self-organization of local inhomogeneities can be found in the competition between the tendency of the electrons to cluster in antiferromagnetic regions, hence producing a shortrange tendency to phase separation, and the longrange Coulomb interaction that instead frustrates it [1,2,3,4].Striped states indeed constitute the best compromise between these competing phenomena and allow the doping holes to be delocalized along linear stripes.
From an experimental point of view, the first evidence of stripes comes from a neutron scattering study on a single crystal of La 1.48 Nd 0.4 Sr 0.12 CuO 4 [5].Since then, a variety of experimental probes, based on neutron scattering, x-ray scattering and scanning tunneling microscopy, pointed to the presence of spin and charge orders [6,7,8,9,10,11].
Nuclear magnetic resonance (NMR) is also a valuable tool to study charge and spin modulations in cuprates.In particular, NMR measures indicated that cuprates that are not La-based may exhibit charge order without spin order [12,13].Moreover, NMR has been employed to investigate the nature of the pseudogap critical point when superconductivity is suppressed [14].
The simplest model that has been considered to reproduce the essential features of the cuprates' phase diagram is the single-band Hubbard model, where only the d x 2 −y 2 orbital of Cu atoms is retained and the impact of oxygen atoms is neglected.However, despite its simplicity, obtaining accurate approximations for the ground state and for low-energy excitations is far from being trivial and several states, very close in energy, have been proposed, obtaining different conclusions from different numerical and analytical methods [15].The model is reported here: where c † R,σ (c R,σ ) creates (destroys) an electron with spin σ on site R and n R,σ = c † R,σ c R,σ is the electron density per spin σ on site R.In the following, we indicate the coordinates of the sites with R = (x, y).It is important that the Hubbard model includes not only the nearest neighbor hopping t and the on-site electron-electron repulsion U , but also the next-nearest neighbor hopping t ′ that has been shown to be a relevant feature in all cuprates, as it constitutes an essential material-dependent parameter, with t ′ /t < 0 [16].The electron density is given by n = N/L, where N is the number of electrons and L is the total number of sites.The hole doping is defined as δ = 1−n.
In the t ′ = 0 case, the presence of stripes in the Hubbard model originates from density-matrix renormalization group (DMRG) studies on 6-leg ladders [17,18] and from further works supporting the idea that charge and spin inhomogeneities may pervade the phase diagram of the Hubbard model [19,20].Charge modulations have been also proposed to be present and to possibly enhance superconductivity by the Dynamic Cluster Approximation and by determinant Quantum Nonte Carlo [21,22].Later, a work which combined a variety of numerical techniques [23], focused on the representative doping δ = 1/8 at U/t = 8, settling that the lowestenergy stripe is a bond-centered one with periodicity λ = 8 (named the stripe wavelength) in the charge sector and 2λ = 16 in the spin sector.
As a consequence, the enlarged unit cell of length λ contains, on average, one hole, as obtained by previous Hartree-Fock calculations [24,25,26,27].Electron pairing was not found in this case and also further studies highlighted the absence of superconductivity at doping δ = 1/8, the system being possibly an insulator [28,29,30].Away from doping δ = 1/8, stripes have been proposed to be stable by different numerical methods, possibly coexisting with superconductivity [28,31,32,33].Recently, an accurate variational auxiliary field quantum Monte Carlo study proposed a phase diagram where stripe phases with no superconductivity are present close to half filling while a superconductive region emerges around δ ∼ 0.2 [29].
In the t ′ ̸ = 0 case, stripes emerge already in the Hartree-Fock approximation, for t ′ /t < 0 [34].Then, their presence is confirmed by different numerical methods, which agree in observing a reduction of the stripe wavelength when t ′ /t < 0, while the stripe wavelength increases when t ′ /t > 0 [35,36,37,38].There is consensus on the absence of superconductivity at δ = 1/8, while different outcomes on the existence of superconductivity are reported when other dopings are considered.The debate is still open since, recently, a combined study based on DMRG and AFQMC indicates that partially filled stripes coexist with superconductivity in a large doping range of the t − t ′ Hubbard model [39], while a DMRG study on 6-leg cylinders suggests that superconducting correlations decay exponentially for t ′ /t < 0 [40].
Alternatively, it has been suggested that superconductivity can be enhanced when models less simplified than the Hubbard one are taken into account.For instance, superconductivity can be recovered, without clear long-range stripe order, with an ab-initio approach that highlights the presence of the realistic offsite interactions [41].Moreover, the three-band Hubbard (or Emery) model has been proposed as a way to enhance superconductivity [42,43], while hopping modulation in a stripe-like manner has been suggested to enhance superconductivity even in the pure Hubbard model [44].Fluctuating stripes have been also proposed to coexist with superconductivity, at difference with the static ones, in the attractive Hubbard model [45].
In this paper, we employ the variational Monte Carlo method with backflow correlations to investigate the effect of doping on stripes and superconductivity in the t − t ′ Hubbard model, for t ′ /t = −0.25 and t ′ /t = −0.40 at U/t = 8.Our simulations are performed on 6-leg ladders, with L x rungs, the total number of sites being L = L x × 6.This geometry has been employed in DMRG calculations and is expected to capture the properties of truly two-dimensional clusters [23], while it allows us to accommodate long stripes along the rungs.First of all, we show that bond-centered and site-centered stripes have similar energies, the only relevant quantity being the stripe wavelength λ.Then, we observe that stripes are the optimal state in a wide doping range: The stripe wavelength reduces at increasing doping, until stripes melt into a uniform state for large values of the doping.We show that superconducting correlations are always suppressed in the presence of stripes, regardless of their insulating or metallic character.Instead, when the stripes melt into a uniform state, a narrow region of superconductivity is observed, around δ ∼ 0.30 for t ′ /t = −0.25 and around δ ∼ 0.21 for t ′ /t = −0.40.We also report on the effect of the next-nearest neighbor hopping, showing that a larger value of |t ′ /t| induces a faster disruption of both stripes and superconductivity as a function of doping.The results discussed here are based on the Master's thesis of the first author of the paper.[46].

Method
Our numerical results are obtained with the variational Monte Carlo method (VMC), which is based on the definition of correlated variational wave functions, whose physical properties can be evaluated within a Monte Carlo scheme [47].In particular, electronelectron correlation is inserted by means of a densitydensity Jastrow factor [48,49] on top of a Slater determinant or a Bardeen-Cooper-Schrieffer (BCS) state.In addition, backflow correlations, as described in [50,51], are implemented; the latter ingredient is important to get accurate variational states.
The wave function is defined by: where J d is the density-density Jastrow factor and |Φ 0 ⟩ is a state that is constructed from the ground state of an auxiliary noninteracting Hamiltonian by applying backflow correlations [50,51].The variational wave functions described below are similar to the ones we built to study the parameter regimes discussed in Refs.[28,38].
The Jastrow factor is given by where n R = σ n R,σ is the electron density on site R and v R,R ′ are pseudopotentials that are optimized for every independent distance |R − R ′ | of the lattice.
In the case of stripe states, the auxiliary noninteracting Hamiltonian is defined as: The first one defines the kinetic energy of the electrons: where the value of the nearest neighbor hopping parameter t is fixed to be equal to the one in the Hubbard Hamiltonian of Eq. ( 1), in order to set the energy scale.Then, the second and third terms describe linear stripes along the x direction that can be either bond-centered or site-centered: and If x 0 = 1/2 the stripes are symmetric with respect to the bond halfway in between two neighboring lattice sites, hence they are called bond centered.Conversely, for x 0 = 0, the stripes are called site centered, as the symmetry axis lies exactly on a lattice site.The periodicity of the charge modulation in both cases is given by λ = 2π/Q.On the other hand, the spin modulation has a π-phase shift across the sites with maximal hole density, resulting in a spin modulation of 2λ = 4π/Q when λ is even and 2π/Q when λ is odd.The spin modulation along the y direction is assumed to have Néel order in all cases.For clarity, we report in Fig. 1 a sketch of the spin modulation along the x direction, for bond-and site-centered stripes, in the case of even and odd wavelengths λ; the length of the arrows is proportional to the size of the local magnetic moments.The effect of the π-shift is clearly visible in the figures.The last term in Eq. ( 4) introduces BCS electron pairing: where the pairing amplitude is modulated in real space: This modulation has been named "in phase" in [52].A chemical potential µ is also included in H BCS .
In the case of uniform states, the auxiliary Hamiltonian is defined as: The kinetic term H 0 is defined as for striped states.
The BCS electron pairing is now defined without modulation in real space, i.e. ∆ R,R+x = ∆ x and ∆ R,R+y = −∆ y , for each site R.In addition, a standard Néel order with pitch vector K = (π, π) can be considered in the uniform state: The auxiliary Hamiltonians of Eqs. ( 4) and ( 10) can be diagonalized by standard methods.Its ground state is then constructed.On top of it, backflow correlations are inserted to define |Φ 0 ⟩ of Eq. ( 2), following our previous works [50,51].
All the parameters in the trial wave function are optimized with the stochastic reconfiguration method [53], in order to minimize the variational energy.In particular, for striped states, we fix a given stripe wavelength λ and optimize ∆ x , ∆ y , µ, ∆ c , ∆ s , and t′ (as well as all the pseudopotentials in the Jastrow factor and the backflow parameters).For a uniform state, we do not have ∆ c and ∆ s as parameters and we optimize instead ∆ AF .Once the energy and all the parameters converge to stable values, the optimization run can be concluded.The values of the parameters are fixed to their averages and a run at fixed parameters is performed to compute the quantum averages needed for the correlation function or the superconducting order parameter.
As already discussed in [28,38], finite values of ∆ c and ∆ s lead effectively to charge and spin modulations, as signaled by a peak (diverging in the thermodynamic limit) at a given Q vector in the static spin and charge structure factors.In order to assess the metallic or the insulating nature of the ground state, we can investigate the small-q behavior of the static charge structure factor N (q), defined as: where ⟨. ..⟩ indicates the expectation value over the variational wave function.Indeed, charge excitations are gapless when N (q) ∝ |q| for |q| → 0, while a charge gap is present whenever N (q) ∝ |q| 2 for |q| → 0 [51,54].The possible existence of superconductivity is investigated by computing correlation functions between Cooper pairs at distance r along the x direction.In particular, we can consider pairs along the y direction, so that: where P R = c R+y,↓ c R,↑ − c R+y,↑ c R,↓ destroys two electrons at nearest-neighbor sites (along y).Then, superconductivity exists whenever D(r) does not decay to zero at large values of r.
Our simulations are performed on ladders with L = L x × 6 sites and periodic boundary conditions along both the x and the y directions.In order to fit charge and spin patterns in the cluster, we take L x = 2kλ (with k integer).Some analysis of the dependence of numerical results on L y and L x can be found in Refs.[28,38].

Results
In this section, we study the instauration of superconductivity and stripes of different wavelength λ when changing the hole doping δ.We consider two typical values of the hopping parameter for cuprates (t ′ /t = −0.25 and t ′ /t = −0.4) in order to see how the value of t ′ /t affects the stripe order.The on-site Coulomb repulsion U/t = 8, kept fixed throughout the simulations, is chosen to ensure strong enough correlations.Indeed, in [38], it is shown that for smaller values of U/t, such as U/t ≲ 4, the striped wave functions are not stable and converge to the uniform state with vanishing parameters ∆ c and ∆ s .
We consider both commensurate and incommensurate doping values.By "commensurate" doping we refer to the introduction of an integer number of holes every 1/δ lattice sites; conversely, this number is noninteger for "incommensurate" doping values.
The optimal variational state is found by comparing, for each considered value of t ′ /t and δ, the variational energies of the striped states for various λ and of the uniform state.The optimal state is then the one with the lowest energy.We start by considering the case t ′ /t = −0.25.The energy per site, in units of t, as a function of δ is reported in Table 1.Here, we compare the energy for the best striped state E stripe /t with that of the uniform state E uniform /t for a broad range of doping values.The striped state is almost always energetically favorable.As δ increases, the wavelength λ decreases more and more until, at doping δ = 1/3, the striped state and the uniform state become energetically indistinguishable.When discussing the behavior of the gap parameters, we will show that this effectively corresponds to a melting of the stripe.
In Sec. 2, we have introduced the gap parameters ∆ c , ∆ s and ∆ AF related to the "strength" of the charge, spin, and Néel order, respectively.In section we proceed by looking at their behavior as the hole-doping increases in the case t ′ /t = −0.25.Their values are plotted in Fig. 2. For small doping, the striped state is well established and indeed ∆ c Table 1.Energy per site (in units of t) for the best striped state E stripe and the uniform state E uniform , along with their relative difference ∆E = E stripe − E uniform , as a function of δ for t ′ /t = −0.25.Data are shown for Lx = 48 at dopings 1/12, 1/6, 1/4, 1/3, for Lx = 45 at doping 1/5, for Lx = 40 at doping 1/8, and for Lx = 70 at doping 1/10.These values of Lx are chosen to accommodate the selected dopings, fit charge and spin patterns of the stripes, and be large enough to have limited size effects.The error bar on the energy, of the order of 10 −4 t, takes into account the weak lattice size dependence of the energies.and ∆ s are finite.This corresponds to a well-defined order in both charge and spin.Also the uniform state, despite not being the optimal one, is able to develop Néel antiferromagnetism in the underdoped regime, as indicated by a finite ∆ AF .
As δ increases, we see that all these parameters decrease monotonically until, at large values of δ they become much smaller and eventually negligible.This corresponds, for the striped states, to the absence of any order: the stripe "melts" and effectively reduces to the uniform state.Hence the degeneracy in energy pointed out for δ = 1/3.Our results confirm the shrinking of the stripes at increasing doping, i.e., shorter modulations are favored when more holes are present in the system.
The metallic or insulating behavior of the optimal state can be assessed from the small-q behavior of the static structure factor N (q), see Eq. ( 12).In particular, we plot the quantity N (q x , 0)/q x at small q x , as shown in Fig. 3.
As a reference, we used a uniform state (squares), which is known to be metallic (except at half-filling, when each site is occupied one electron and the Coulomb repulsion prevents them from moving freely) even though it has a higher variational energy.We observe that, for the striped state at δ = 1/6 (diamonds), N (q x , 0)/q x clearly tends to zero, compatibly with an insulating behavior.On the other hand, for all the other striped states at δ = 1/5 and δ = 1/4 (circles), N (q x , 0)/q x tends to a finite value indicating that these states are metallic.
Finally, we address the coexistence of superconductivity and stripe order, by computing the superconducting order parameter of Eq. ( 13).In Fig. 4, we compare the uniform (but not optimal) state at δ = 1/6 (blue squares), taken as a reference, with some optimal striped states.All the striped states show strongly suppressed pair-pair correlations with respect to the uniform case.The stripes at δ = 1/5 and δ = 1/4(circles), despite having a metallic character, exhibit a suppression in D(r) similar to that of the insulating stripe at δ = 1/6 (diamonds).This supports the idea that the stripe order disrupts superconductivity, no matter their metallic or insulating character.
Since stripes are found to compete with superconductivity, we investigate then whether there is a region where hole-doping is strong enough to restore the uniform state but not too strong to suppress superconductivity.In order to answer this question, we look for the first value of δ at which the uniform state becomes energetically favorable and compute the pair-pair correlations.For t ′ /t = −0.25, as discussed in Table 1, the optimal state at δ = 1/4 is a stripe of wavelength λ = 3 while at δ = 1/3 we have already reached the uniform state.We then study incommensurate values of δ in the range 1  4 , 1 3 .Since the wavelengths of the stripes decrease at increasing doping, it is sufficient to compare the striped state with λ = 3 and the uniform state in this doping regime.Their variational energies are presented in Table 2.
We can identify as the critical doping, the value δ c = 0.29, where the optimal parameters of the stripe state ∆ c and ∆ s vanish.In Fig. 4, lower panel, the pair-pair correlations for this state (black squares) are plotted next to the uniform, but not optimal, superconducting state at δ = 1/6 (blue squares) and the nonsuperconducting striped state at δ = 1/6 with λ = 4 (diamonds), for comparison.In this "intermediate" state, superconductivity is suppressed with respect to the uniform state at smaller doping, due to an already strong hole-doping, but is still present,   The second set of simulations involved the same search for the optimal state but at a larger value of |t ′ /t|, namely t ′ /t = −0.4.The main effect of a larger value of |t ′ /t| is to suppress the stripe pattern and makes the optimal state converge to the uniform one faster.Indeed, while at δ = 1/5 the optimal state is a stripe with λ = 3, already at δ = 1/4 the striped state is no longer favorable and the uniform state is the optimal one.We can ascribe this behavior to the larger degree of frustration that is present for a larger value of |t ′ /t|.We investigate then whether the suppression of the stripes at a lower concentration of holes might be associated to the presence of stronger superconducting correlations around the critical doping δ c .Following the same reasoning as before, the doping at which the uniform state prevails is for δ in the range 1 5 , 1 4 .In particular, stripes are suppressed already at δ c = 0.21.In Fig. 5, we show the pair-pair correlations for the case t ′ /t = −0.4 for different values of doping.Again, we can see how correlations for the uniform state at the critical doping (full squares) are suppressed with respect to the uniform but not optimal, superconducting state at δ = 1/6 (empty squares), but still stronger than in the cases where stripe order is established.
To conclude the discussion, we compare the magnitude of the superconducting correlations for the two different values of t ′ /t.The curves are plotted in Fig. 6.We observe how, when |t ′ /t| is larger, superconductivity for the (nonoptimal) uniform state at δ = 1/6 is slightly suppressed.This is equally true for the uniform states at the critical dopings δ c where the uniform state is restored: The effect of a larger value of |t ′ /t| in suppressing superconductivity is clearly visible, since superconducting correlations at the critical dopings are smaller at t ′ /t = −0.4than at  t ′ /t = −0.25,even if δ c is definitely smaller in the first case.

Conclusions
We have explored the consequences of increasing hole doping on the instauration of stripe order, superconductivity and their reciprocal interplay, for two prototypical values of t ′ /t.All the main results of the present work are collected and summarized by the final phase diagram reported in Fig. 7. Superconductivity is considered to be present when the average of D(r) over the last 10 distances is above a threshold value of 3 × 10 −4 §.
By looking for the optimal state for different values of the hole-doping δ, we found that stripes are present over a broad range of doping values, as they are energetically favorable in comparison to the uniform state.Remarkably, site and bond-centered stripes have been found to be essentially degenerate in energy, suggesting that there is no relevant difference between the two configurations.Upon increasing δ, the wavelength of the stripes shrinks until eventually the uniform state is restored.A larger |t ′ /t| is associated to a shrinking of the wavelength λ and leads to a faster dissolution of the stripes, with the uniform state being the optimal one at a smaller value of δ.
The coexistence of superconductivity and stripe order is addressed by looking at the pair-pair superconducting correlations D(r).For both values of t ′ /t, superconductivity is found to be suppressed whenever stripes (no matter their metallic or insulating nature) are present, suggesting that the two phenomena interfere with each other.There is then a small interval in δ among which the hole-doping is strong enough to restore the uniform state but not too strong to completely suppress superconductivity.Furthermore, our results show that, at t ′ /t = −0.4all superconducting correlations are weaker than at t ′ /t = −0.25,even if stripes melt at a smaller value of δ.
In conclusion, our results confirm that the phase diagram of the Hubbard model is dominated by stripe states, possibly overestimating this phase with respect to superconductivity, when connected to the cuprate physics.

Figure 1 .
Figure1.Top left panel: spin modulation of a bond-centered stripe with even wavelength.The minimum and the maximum of the local magnetic momenta are located on bonds, as highlighted by the red rectangles.Top right panel: spin modulation of a site-centered stripe with even wavelength.The minimum and the maximum of the local magnetic momenta are located on sites, as highlighted by the blue rectangles.Bottom panels: spin modulations of a stripe with odd wavelength; in the left panel the modulation is bond-centered where the local magnetic moment is minimal (red rectangle) and site-centered where the local magnetic moment is maximal (blue rectangle); in the right panel the situation is reversed.

Figure 2 .
Figure 2. Behavior of ∆c (full squares) and ∆s (empty squares) for the best striped state and ∆ AF (circles) for the uniform state, as a function of δ for t ′ /t = −0.25.The best striped state for each doping is shown in figure.The error bars are smaller than the symbol size.

Figure 4 .
Figure 4. Pair-pair correlations D(r) as a function of r on a log-log scale, at t ′ /t = −0.25.Upper panel: Data are reported for the optimal stripe states at different hole-dopings δ (circles) and for the (nonoptimal) uniform state at doping δ = 1/6 (blue squares).Lower panel: Data are reported for the uniform state at the critical doping δc = 0.29 (black squares) along with the (nonoptimal) uniform superconducting state at δ = 1/6 (blue squares) and the nonsuperconducting striped state with λ = 4 at δ = 1/6 (diamonds).

Figure 5 .
Figure 5. Pair-pair correlations D(r) as a function of r on a loglog scale.Data are reported at t ′ /t = −0.4 for the (nonoptimal) uniform state at doping δ = 1/6 (empty squares) and for the optimal striped states at dopings δ = 1/6 (diamonds), 1/5 (circles), as well as for the uniform state at the critical doping δc = 0.21 (full squares).

Figure 6 .
Figure 6.Comparison of the pair-pair correlations D(r) for t ′ /t = −0.25 and t ′ /t = −0.4,as a function of r on a loglog scale.Data are reported for the uniform but not optimal superconducting states at δ = 1/6 (empty red squares at t ′ /t = −0.25 and empty blue squares at t ′ /t = −0.4) and for the two uniform states at the critical dopings δc = 0.29 for t ′ /t = −0.25 (full red squares) and δc = 0.21 for t ′ /t = −0.4(full blue squares).

Table 2 .
Energy per site (in units of t) for the best striped state