Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

Variational Monte Carlo study of stripes as a function of doping in the tt𝑡superscript𝑡t-t^{\prime}italic_t - italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT Hubbard model

Antonio Lechiara, Vito Marino, Luca F. Tocchio111Author to whom any correspondence should be addressed Institute for Condensed Matter Physics and Complex Systems, DISAT, Politecnico di Torino, I-10129 Torino, Italy luca.tocchio@polito.it
(July 8, 2024)
Abstract

We perform variational Monte Carlo simulations of the single-band Hubbard model on the square lattice with both nearest (t𝑡titalic_t) and next-nearest (tsuperscript𝑡t^{\prime}italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT) 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.

\ioptwocol

1 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 short-range tendency to phase separation, and the long-range 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 La1.48Nd0.4Sr0.12CuO4 [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 dx2y2subscript𝑑superscript𝑥2superscript𝑦2d_{x^{2}-y^{2}}italic_d start_POSTSUBSCRIPT italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_y start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_POSTSUBSCRIPT 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:

=tR,R,σcR,σcR,σ𝑡subscript𝑅superscript𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absentsuperscript𝑅𝜎\displaystyle{\cal H}=-t\sum_{\langle R,R^{\prime}\rangle,\sigma}c^{\dagger}_{% R,\sigma}c^{\phantom{\dagger}}_{R^{\prime},\sigma}caligraphic_H = - italic_t ∑ start_POSTSUBSCRIPT ⟨ italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ⟩ , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_σ end_POSTSUBSCRIPT
tR,R,σcR,σcR,σ+H.c.+URnR,nR,,superscript𝑡subscriptdelimited-⟨⟩𝑅superscript𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absentsuperscript𝑅𝜎H.c.𝑈subscript𝑅subscript𝑛𝑅subscript𝑛𝑅\displaystyle-t^{\prime}\sum_{\langle\langle R,R^{\prime}\rangle\rangle,\sigma% }c^{\dagger}_{R,\sigma}c^{\phantom{\dagger}}_{R^{\prime},\sigma}+\textrm{H.c.}% +U\sum_{R}n_{R,\uparrow}n_{R,\downarrow},- italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT ⟨ ⟨ italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ⟩ ⟩ , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_σ end_POSTSUBSCRIPT + H.c. + italic_U ∑ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT , (1)

where cR,σsubscriptsuperscript𝑐𝑅𝜎c^{\dagger}_{R,\sigma}italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT (cR,σsubscriptsuperscript𝑐absent𝑅𝜎c^{\phantom{\dagger}}_{R,\sigma}italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT) creates (destroys) an electron with spin σ𝜎\sigmaitalic_σ on site R𝑅Ritalic_R and nR,σ=cR,σcR,σsubscript𝑛𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absent𝑅𝜎n_{R,\sigma}=c^{\dagger}_{R,\sigma}c^{\phantom{\dagger}}_{R,\sigma}italic_n start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT = italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT is the electron density per spin σ𝜎\sigmaitalic_σ on site R𝑅Ritalic_R. In the following, we indicate the coordinates of the sites with 𝐑=(x,y)𝐑𝑥𝑦{\bf R}=(x,y)bold_R = ( italic_x , italic_y ). It is important that the Hubbard model includes not only the nearest neighbor hopping t𝑡titalic_t and the on-site electron-electron repulsion U𝑈Uitalic_U, but also the next-nearest neighbor hopping tsuperscript𝑡t^{\prime}italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT that has been shown to be a relevant feature in all cuprates, as it constitutes an essential material-dependent parameter, with t/t<0superscript𝑡𝑡0t^{\prime}/t<0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t < 0 [16]. The electron density is given by n=N/L𝑛𝑁𝐿n=N/Litalic_n = italic_N / italic_L, where N𝑁Nitalic_N is the number of electrons and L𝐿Litalic_L is the total number of sites. The hole doping is defined as δ=1n𝛿1𝑛\delta=1-nitalic_δ = 1 - italic_n.

In the t=0superscript𝑡0t^{\prime}=0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT = 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𝛿18\delta=1/8italic_δ = 1 / 8 at U/t=8𝑈𝑡8U/t=8italic_U / italic_t = 8, settling that the lowest-energy stripe is a bond-centered one with periodicity λ=8𝜆8\lambda=8italic_λ = 8 (named the stripe wavelength) in the charge sector and 2λ=162𝜆162\lambda=162 italic_λ = 16 in the spin sector. As a consequence, the enlarged unit cell of length λ𝜆\lambdaitalic_λ 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𝛿18\delta=1/8italic_δ = 1 / 8, the system being possibly an insulator [28, 29, 30]. Away from doping δ=1/8𝛿18\delta=1/8italic_δ = 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.2similar-to𝛿0.2\delta\sim 0.2italic_δ ∼ 0.2 [29].

In the t0superscript𝑡0t^{\prime}\neq 0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ≠ 0 case, stripes emerge already in the Hartree-Fock approximation, for t/t<0superscript𝑡𝑡0t^{\prime}/t<0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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<0superscript𝑡𝑡0t^{\prime}/t<0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t < 0, while the stripe wavelength increases when t/t>0superscript𝑡𝑡0t^{\prime}/t>0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t > 0 [35, 36, 37, 38]. There is consensus on the absence of superconductivity at δ=1/8𝛿18\delta=1/8italic_δ = 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 tt𝑡superscript𝑡t-t^{\prime}italic_t - italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT Hubbard model [39], while a DMRG study on 6-leg cylinders suggests that superconducting correlations decay exponentially for t/t<0superscript𝑡𝑡0t^{\prime}/t<0italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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 off-site 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 tt𝑡superscript𝑡t-t^{\prime}italic_t - italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT Hubbard model, for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 and t/t=0.40superscript𝑡𝑡0.40t^{\prime}/t=-0.40italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.40 at U/t=8𝑈𝑡8U/t=8italic_U / italic_t = 8. Our simulations are performed on 6-leg ladders, with Lxsubscript𝐿𝑥L_{x}italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT rungs, the total number of sites being L=Lx×6𝐿subscript𝐿𝑥6L=L_{x}\times 6italic_L = italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT × 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 λ𝜆\lambdaitalic_λ. 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.30similar-to𝛿0.30\delta\sim 0.30italic_δ ∼ 0.30 for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 and around δ0.21similar-to𝛿0.21\delta\sim 0.21italic_δ ∼ 0.21 for t/t=0.40superscript𝑡𝑡0.40t^{\prime}/t=-0.40italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.40. We also report on the effect of the next-nearest neighbor hopping, showing that a larger value of |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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].

2 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, electron-electron correlation is inserted by means of a density-density 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:

|Ψ=𝒥d|Φ0,ketΨsubscript𝒥𝑑ketsubscriptΦ0|\Psi\rangle={\cal J}_{d}|\Phi_{0}\rangle\,,| roman_Ψ ⟩ = caligraphic_J start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT | roman_Φ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ , (2)

where 𝒥dsubscript𝒥𝑑{\cal J}_{d}caligraphic_J start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT is the density-density Jastrow factor and |Φ0ketsubscriptΦ0|\Phi_{0}\rangle| roman_Φ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ 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

𝒥d=exp(12R,RvR,RnRnR),subscript𝒥𝑑12subscript𝑅superscript𝑅subscript𝑣𝑅superscript𝑅subscript𝑛𝑅subscript𝑛superscript𝑅{\cal J}_{d}=\exp\left(-\frac{1}{2}\sum_{R,R^{\prime}}v_{R,R^{\prime}}n_{R}n_{% R^{\prime}}\right)\,,caligraphic_J start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT = roman_exp ( - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∑ start_POSTSUBSCRIPT italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ) , (3)

where nR=σnR,σsubscript𝑛𝑅subscript𝜎subscript𝑛𝑅𝜎n_{R}=\sum_{\sigma}n_{R,\sigma}italic_n start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT italic_σ end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT is the electron density on site R𝑅Ritalic_R and vR,Rsubscript𝑣𝑅superscript𝑅v_{R,R^{\prime}}italic_v start_POSTSUBSCRIPT italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT are pseudopotentials that are optimized for every independent distance |𝐑𝐑|𝐑superscript𝐑|{\bf R}-{\bf R}^{\prime}|| bold_R - bold_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT | of the lattice.

In the case of stripe states, the auxiliary noninteracting Hamiltonian is defined as:

aux=0+charge+spin+BCS.subscriptauxsubscript0subscriptchargesubscriptspinsubscriptBCS{\cal H}_{\rm{aux}}={\cal H}_{0}+{\cal H}_{\rm{charge}}+{\cal H}_{\rm{spin}}+{% \cal H}_{\rm{BCS}}.caligraphic_H start_POSTSUBSCRIPT roman_aux end_POSTSUBSCRIPT = caligraphic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + caligraphic_H start_POSTSUBSCRIPT roman_charge end_POSTSUBSCRIPT + caligraphic_H start_POSTSUBSCRIPT roman_spin end_POSTSUBSCRIPT + caligraphic_H start_POSTSUBSCRIPT roman_BCS end_POSTSUBSCRIPT . (4)

The first one defines the kinetic energy of the electrons:

0=tR,R,σcR,σcR,σt~R,R,σcR,σcR,σ+H.c.,subscript0𝑡subscript𝑅superscript𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absentsuperscript𝑅𝜎superscript~𝑡subscriptdelimited-⟨⟩𝑅superscript𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absentsuperscript𝑅𝜎H.c.{\cal H}_{0}=-t\sum_{\langle R,R^{\prime}\rangle,\sigma}c^{\dagger}_{R,\sigma}% c^{\phantom{\dagger}}_{R^{\prime},\sigma}-\tilde{t}^{\prime}\sum_{\langle% \langle R,R^{\prime}\rangle\rangle,\sigma}c^{\dagger}_{R,\sigma}c^{\phantom{% \dagger}}_{R^{\prime},\sigma}+\textrm{H.c.},caligraphic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = - italic_t ∑ start_POSTSUBSCRIPT ⟨ italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ⟩ , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_σ end_POSTSUBSCRIPT - over~ start_ARG italic_t end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT ⟨ ⟨ italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ⟩ ⟩ , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT , italic_σ end_POSTSUBSCRIPT + H.c. , (5)

where the value of the nearest neighbor hopping parameter t𝑡titalic_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𝑥xitalic_x direction that can be either bond-centered or site-centered:

charge=ΔcRcos[Q(xx0)](cR,cR,+cR,cR,),subscriptchargesubscriptΔcsubscript𝑅𝑄𝑥subscript𝑥0subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅{\cal H}_{\rm{charge}}=\Delta_{\rm{c}}\sum_{R}\cos\left[Q\left(x-x_{0}\right)% \right]\left(c^{\dagger}_{R,\uparrow}c^{\phantom{\dagger}}_{R,\uparrow}+c^{% \dagger}_{R,\downarrow}c^{\phantom{\dagger}}_{R,\downarrow}\right),caligraphic_H start_POSTSUBSCRIPT roman_charge end_POSTSUBSCRIPT = roman_Δ start_POSTSUBSCRIPT roman_c end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT roman_cos [ italic_Q ( italic_x - italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] ( italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT + italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT ) , (6)

and

spin=ΔsR(1)x+ysin[Q2(xx0)]×\displaystyle{\cal H}_{\rm{spin}}=\Delta_{\rm{s}}\sum_{R}(-1)^{x+y}\sin\left[% \frac{Q}{2}\left(x-x_{0}\right)\right]\timescaligraphic_H start_POSTSUBSCRIPT roman_spin end_POSTSUBSCRIPT = roman_Δ start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ( - 1 ) start_POSTSUPERSCRIPT italic_x + italic_y end_POSTSUPERSCRIPT roman_sin [ divide start_ARG italic_Q end_ARG start_ARG 2 end_ARG ( italic_x - italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] ×
×(cR,cR,cR,cR,).absentsubscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅\displaystyle\times\left(c^{\dagger}_{R,\uparrow}c^{\phantom{\dagger}}_{R,% \uparrow}-c^{\dagger}_{R,\downarrow}c^{\phantom{\dagger}}_{R,\downarrow}\right).× ( italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT - italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT ) . (7)

If x0=1/2subscript𝑥012x_{0}=1/2italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 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 x0=0subscript𝑥00x_{0}=0italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 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𝜆2𝜋𝑄\lambda=2\pi/Qitalic_λ = 2 italic_π / italic_Q. On the other hand, the spin modulation has a π𝜋\piitalic_π-phase shift across the sites with maximal hole density, resulting in a spin modulation of 2λ=4π/Q2𝜆4𝜋𝑄2\lambda=4\pi/Q2 italic_λ = 4 italic_π / italic_Q when λ𝜆\lambdaitalic_λ is even and 2π/Q2𝜋𝑄2\pi/Q2 italic_π / italic_Q when λ𝜆\lambdaitalic_λ is odd. The spin modulation along the y𝑦yitalic_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𝑥xitalic_x direction, for bond- and site-centered stripes, in the case of even and odd wavelengths λ𝜆\lambdaitalic_λ; the length of the arrows is proportional to the size of the local magnetic moments. The effect of the π𝜋\piitalic_π-shift is clearly visible in the figures. The last term in Eq. (4) introduces BCS electron pairing:

BCS=R,η=x,yΔR,R+η(cR,cR+η,cR,cR+η,)+H.c.subscriptBCSsubscriptformulae-sequence𝑅𝜂𝑥𝑦subscriptΔ𝑅𝑅𝜂subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐𝑅𝜂subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐𝑅𝜂H.c.\displaystyle{\cal H}_{\rm{BCS}}=\sum_{R,\eta=x,y}\Delta_{R,R+\eta}(c^{\dagger% }_{R,\uparrow}c^{\dagger}_{R+\eta,\downarrow}-c^{\dagger}_{R,\downarrow}c^{% \dagger}_{R+\eta,\uparrow})+\textrm{H.c.}caligraphic_H start_POSTSUBSCRIPT roman_BCS end_POSTSUBSCRIPT = ∑ start_POSTSUBSCRIPT italic_R , italic_η = italic_x , italic_y end_POSTSUBSCRIPT roman_Δ start_POSTSUBSCRIPT italic_R , italic_R + italic_η end_POSTSUBSCRIPT ( italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R + italic_η , ↓ end_POSTSUBSCRIPT - italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R + italic_η , ↑ end_POSTSUBSCRIPT ) + H.c.
μR,σcR,σcR,σ,𝜇subscript𝑅𝜎subscriptsuperscript𝑐𝑅𝜎subscriptsuperscript𝑐absent𝑅𝜎\displaystyle-\mu\sum_{R,\sigma}c^{\dagger}_{R,\sigma}c^{\phantom{\dagger}}_{R% ,\sigma},- italic_μ ∑ start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , italic_σ end_POSTSUBSCRIPT , (8)

where the pairing amplitude is modulated in real space:

ΔR,R+x=Δx|cos[Q2(x+12x0)]|subscriptΔ𝑅𝑅𝑥subscriptΔ𝑥𝑄2𝑥12subscript𝑥0\displaystyle\Delta_{R,R+x}=\Delta_{x}\left|\cos\left[\frac{Q}{2}(x+\frac{1}{2% }-x_{0})\right]\right|roman_Δ start_POSTSUBSCRIPT italic_R , italic_R + italic_x end_POSTSUBSCRIPT = roman_Δ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT | roman_cos [ divide start_ARG italic_Q end_ARG start_ARG 2 end_ARG ( italic_x + divide start_ARG 1 end_ARG start_ARG 2 end_ARG - italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] |
ΔR,R+y=Δy|cos[Q2(xx0)]|.subscriptΔ𝑅𝑅𝑦subscriptΔ𝑦𝑄2𝑥subscript𝑥0\displaystyle\Delta_{R,R+y}=-\Delta_{y}\left|\cos\left[\frac{Q}{2}\left(x-x_{0% }\right)\right]\right|.roman_Δ start_POSTSUBSCRIPT italic_R , italic_R + italic_y end_POSTSUBSCRIPT = - roman_Δ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT | roman_cos [ divide start_ARG italic_Q end_ARG start_ARG 2 end_ARG ( italic_x - italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) ] | . (9)

This modulation has been named “in phase” in [52]. A chemical potential μ𝜇\muitalic_μ is also included in BCSsubscriptBCS{\cal H}_{\rm{BCS}}caligraphic_H start_POSTSUBSCRIPT roman_BCS end_POSTSUBSCRIPT.

Refer to caption
Figure 1: 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.

In the case of uniform states, the auxiliary Hamiltonian is defined as:

aux=0+BCS+AF.subscriptauxsubscript0subscriptBCSsubscriptAF{\cal H}_{\rm{aux}}={\cal H}_{0}+{\cal H}_{\rm{BCS}}+{\cal H}_{\rm{AF}}.caligraphic_H start_POSTSUBSCRIPT roman_aux end_POSTSUBSCRIPT = caligraphic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + caligraphic_H start_POSTSUBSCRIPT roman_BCS end_POSTSUBSCRIPT + caligraphic_H start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT . (10)

The kinetic term 0subscript0{\cal H}_{0}caligraphic_H start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is defined as for striped states. The BCS electron pairing is now defined without modulation in real space, i.e. ΔR,R+x=ΔxsubscriptΔ𝑅𝑅𝑥subscriptΔ𝑥\Delta_{R,R+x}=\Delta_{x}roman_Δ start_POSTSUBSCRIPT italic_R , italic_R + italic_x end_POSTSUBSCRIPT = roman_Δ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT and ΔR,R+y=ΔysubscriptΔ𝑅𝑅𝑦subscriptΔ𝑦\Delta_{R,R+y}=-\Delta_{y}roman_Δ start_POSTSUBSCRIPT italic_R , italic_R + italic_y end_POSTSUBSCRIPT = - roman_Δ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT, for each site R𝑅Ritalic_R. In addition, a standard Néel order with pitch vector 𝐊=(π,π)𝐊𝜋𝜋{\bf K}=(\pi,\pi)bold_K = ( italic_π , italic_π ) can be considered in the uniform state:

AF=ΔAFR(1)x+y(cR,cR,cR,cR,).subscriptAFsubscriptΔAFsubscript𝑅superscript1𝑥𝑦subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅subscriptsuperscript𝑐𝑅subscriptsuperscript𝑐absent𝑅{\cal H}_{\rm{AF}}=\Delta_{\rm{AF}}\sum_{R}(-1)^{x+y}\left(c^{\dagger}_{R,% \uparrow}c^{\phantom{\dagger}}_{R,\uparrow}-c^{\dagger}_{R,\downarrow}c^{% \phantom{\dagger}}_{R,\downarrow}\right).caligraphic_H start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT = roman_Δ start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ( - 1 ) start_POSTSUPERSCRIPT italic_x + italic_y end_POSTSUPERSCRIPT ( italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT - italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT ) . (11)

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 |Φ0ketsubscriptΦ0|\Phi_{0}\rangle| roman_Φ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⟩ 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 λ𝜆\lambdaitalic_λ and optimize ΔxsubscriptΔ𝑥\Delta_{x}roman_Δ start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT, ΔysubscriptΔ𝑦\Delta_{y}roman_Δ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT, μ𝜇\muitalic_μ, ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT, and t~superscript~𝑡\tilde{t}^{\prime}over~ start_ARG italic_t end_ARG start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT (as well as all the pseudopotentials in the Jastrow factor and the backflow parameters). For a uniform state, we do not have ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT as parameters and we optimize instead ΔAFsubscriptΔAF\Delta_{\rm{AF}}roman_Δ start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT. 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 ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT lead effectively to charge and spin modulations, as signaled by a peak (diverging in the thermodynamic limit) at a given 𝐐𝐐{\bf Q}bold_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𝑞qitalic_q behavior of the static charge structure factor N(𝐪)𝑁𝐪N({\bf q})italic_N ( bold_q ), defined as:

N(𝐪)=1LR,RnRnRei𝐪(𝐑𝐑),𝑁𝐪1𝐿subscript𝑅superscript𝑅delimited-⟨⟩subscript𝑛𝑅subscript𝑛superscript𝑅superscript𝑒𝑖𝐪𝐑superscript𝐑N({\bf q})=\frac{1}{L}\sum_{R,R^{\prime}}\langle n_{R}n_{R^{\prime}}\rangle e^% {i{\bf q}\cdot({\bf R}-{\bf R^{\prime}})},italic_N ( bold_q ) = divide start_ARG 1 end_ARG start_ARG italic_L end_ARG ∑ start_POSTSUBSCRIPT italic_R , italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ⟨ italic_n start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT italic_n start_POSTSUBSCRIPT italic_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT end_POSTSUBSCRIPT ⟩ italic_e start_POSTSUPERSCRIPT italic_i bold_q ⋅ ( bold_R - bold_R start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ) end_POSTSUPERSCRIPT , (12)

where delimited-⟨⟩\langle\dots\rangle⟨ … ⟩ indicates the expectation value over the variational wave function. Indeed, charge excitations are gapless when N(𝐪)|𝐪|proportional-to𝑁𝐪𝐪N({\bf q})\propto|{\bf q}|italic_N ( bold_q ) ∝ | bold_q | for |𝐪|0𝐪0|{\bf q}|\to 0| bold_q | → 0, while a charge gap is present whenever N(𝐪)|𝐪|2proportional-to𝑁𝐪superscript𝐪2N({\bf q})\propto|{\bf q}|^{2}italic_N ( bold_q ) ∝ | bold_q | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT for |𝐪|0𝐪0|{\bf q}|\to 0| bold_q | → 0 [51, 54].

The possible existence of superconductivity is investigated by computing correlation functions between Cooper pairs at distance r𝑟ritalic_r along the x𝑥xitalic_x direction. In particular, we can consider pairs along the y𝑦yitalic_y direction, so that:

D(r)=1LRPRPR+rx,𝐷𝑟1𝐿subscript𝑅delimited-⟨⟩subscriptsuperscript𝑃absent𝑅subscriptsuperscript𝑃𝑅𝑟𝑥D(r)=\frac{1}{L}\sum_{R}\langle P^{\phantom{\dagger}}_{R}P^{\dagger}_{R+rx}\rangle,italic_D ( italic_r ) = divide start_ARG 1 end_ARG start_ARG italic_L end_ARG ∑ start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT ⟨ italic_P start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT italic_P start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R + italic_r italic_x end_POSTSUBSCRIPT ⟩ , (13)

where PR=cR+y,cR,cR+y,cR,subscriptsuperscript𝑃absent𝑅subscriptsuperscript𝑐absent𝑅𝑦subscriptsuperscript𝑐absent𝑅subscriptsuperscript𝑐absent𝑅𝑦subscriptsuperscript𝑐absent𝑅P^{\phantom{\dagger}}_{R}=c^{\phantom{\dagger}}_{R+y,\downarrow}c^{\phantom{% \dagger}}_{R,\uparrow}-c^{\phantom{\dagger}}_{R+y,\uparrow}c^{\phantom{\dagger% }}_{R,\downarrow}italic_P start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT = italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R + italic_y , ↓ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↑ end_POSTSUBSCRIPT - italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R + italic_y , ↑ end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_R , ↓ end_POSTSUBSCRIPT destroys two electrons at nearest-neighbor sites (along y𝑦yitalic_y). Then, superconductivity exists whenever D(r)𝐷𝑟D(r)italic_D ( italic_r ) does not decay to zero at large values of r𝑟ritalic_r.

Our simulations are performed on ladders with L=Lx×6𝐿subscript𝐿𝑥6L=L_{x}\times 6italic_L = italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT × 6 sites and periodic boundary conditions along both the x𝑥xitalic_x and the y𝑦yitalic_y directions. In order to fit charge and spin patterns in the cluster, we take Lx=2kλsubscript𝐿𝑥2𝑘𝜆L_{x}=2k\lambdaitalic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 2 italic_k italic_λ (with k𝑘kitalic_k integer). Some analysis of the dependence of numerical results on Lysubscript𝐿𝑦L_{y}italic_L start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT and Lxsubscript𝐿𝑥L_{x}italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT can be found in Refs. [28, 38].

3 Results

In this section, we study the instauration of superconductivity and stripes of different wavelength λ𝜆\lambdaitalic_λ when changing the hole doping δ𝛿\deltaitalic_δ. We consider two typical values of the hopping parameter for cuprates (t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 and t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4) in order to see how the value of t/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t affects the stripe order. The on-site Coulomb repulsion U/t=8𝑈𝑡8U/t=8italic_U / italic_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𝑈𝑡U/titalic_U / italic_t, such as U/t4less-than-or-similar-to𝑈𝑡4U/t\lesssim 4italic_U / italic_t ≲ 4, the striped wave functions are not stable and converge to the uniform state with vanishing parameters ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT.

We consider both commensurate and incommensurate doping values. By ”commensurate” doping we refer to the introduction of an integer number of holes every 1/δ1𝛿1/\delta1 / italic_δ 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/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t and δ𝛿\deltaitalic_δ, the variational energies of the striped states for various λ𝜆\lambdaitalic_λ and of the uniform state. The optimal state is then the one with the lowest energy.

3.1 t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25

We start by considering the case t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. The energy per site, in units of t𝑡titalic_t, as a function of δ𝛿\deltaitalic_δ is reported in Table 1. Here, we compare the energy for the best striped state Estripe/tsubscript𝐸stripe𝑡E_{\rm{stripe}}/titalic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT / italic_t with that of the uniform state Euniform/tsubscript𝐸uniform𝑡E_{\rm{uniform}}/titalic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT / italic_t for a broad range of doping values. The striped state is almost always energetically favorable. As δ𝛿\deltaitalic_δ increases, the wavelength λ𝜆\lambdaitalic_λ decreases more and more until, at doping δ=1/3𝛿13\delta=1/3italic_δ = 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.

Table 1: Energy per site (in units of t𝑡titalic_t) for the best striped state Estripesubscript𝐸stripeE_{\rm{stripe}}italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT and the uniform state Euniformsubscript𝐸uniformE_{\rm{uniform}}italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT, along with their relative difference ΔE=EstripeEuniformΔ𝐸subscript𝐸stripesubscript𝐸uniform\Delta E=E_{\rm{stripe}}-E_{\rm{uniform}}roman_Δ italic_E = italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT - italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT, as a function of δ𝛿\deltaitalic_δ for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. Data are shown for Lx=48subscript𝐿𝑥48L_{x}=48italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 48 at dopings 1/12,1/6,1/4,1/31121614131/12,1/6,1/4,1/31 / 12 , 1 / 6 , 1 / 4 , 1 / 3, for Lx=45subscript𝐿𝑥45L_{x}=45italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 45 at doping 1/5, for Lx=40subscript𝐿𝑥40L_{x}=40italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 40 at doping 1/8, and for Lx=70subscript𝐿𝑥70L_{x}=70italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 70 at doping 1/10. These values of Lxsubscript𝐿𝑥L_{x}italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT 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 104tsuperscript104𝑡10^{-4}t10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT italic_t, takes into account the weak lattice size dependence of the energies.
δ𝛿\deltaitalic_δ Estripesubscript𝐸stripeE_{\rm{stripe}}italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT Euniformsubscript𝐸uniformE_{\rm{uniform}}italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT ΔEΔ𝐸\Delta Eroman_Δ italic_E
1/12 -0.6646 (λ=8𝜆8\lambda=8italic_λ = 8) -0.6572 -0.0074
1/10 -0.6920 (λ=7𝜆7\lambda=7italic_λ = 7) -0.6842 -0.0078
1/8 -0.7322 (λ=5𝜆5\lambda=5italic_λ = 5) -0.7238 -0.0084
1/6 -0.7936 (λ=4𝜆4\lambda=4italic_λ = 4) -0.7847 -0.0088
1/5 -0.8281 (λ=3𝜆3\lambda=3italic_λ = 3) -0.8258 -0.0023
1/4 -0.8749 (λ=3𝜆3\lambda=3italic_λ = 3) -0.8727 -0.0021
1/3 -0.9197 (λ=3𝜆3\lambda=3italic_λ = 3) -0.9197 0

In Sec. 2, we have introduced the gap parameters ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT and ΔAFsubscriptΔAF\Delta_{\rm{AF}}roman_Δ start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT related to the ”strength” of the charge, spin, and Néel order, respectively. In this section we proceed by looking at their behavior as the hole-doping increases in the case t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. Their values are plotted in Fig. 2. For small doping, the striped state is well established and indeed ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT 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 ΔAFsubscriptΔAF\Delta_{\rm{AF}}roman_Δ start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT.

As δ𝛿\deltaitalic_δ increases, we see that all these parameters decrease monotonically until, at large values of δ𝛿\deltaitalic_δ 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𝛿13\delta=1/3italic_δ = 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.

Refer to caption
Figure 2: Behavior of ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT (full squares) and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT (empty squares) for the best striped state and ΔAFsubscriptΔAF\Delta_{\rm{AF}}roman_Δ start_POSTSUBSCRIPT roman_AF end_POSTSUBSCRIPT (circles) for the uniform state, as a function of δ𝛿\deltaitalic_δ for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. The best striped state for each doping is shown in figure. The error bars are smaller than the symbol size.

The metallic or insulating behavior of the optimal state can be assessed from the small-q𝑞qitalic_q behavior of the static structure factor N(q)𝑁qN(\textbf{q})italic_N ( q ), see Eq. (12). In particular, we plot the quantity N(qx,0)/qx𝑁subscript𝑞𝑥0subscript𝑞𝑥N(q_{x},0)/q_{x}italic_N ( italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , 0 ) / italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT at small qxsubscript𝑞𝑥q_{x}italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT, as shown in Fig. 3.

Refer to caption
Figure 3: Static structure factor (divided by qxsubscript𝑞𝑥q_{x}italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT) N(q)/qx𝑁qsubscript𝑞𝑥N(\textbf{q})/q_{x}italic_N ( q ) / italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT as a function of qxsubscript𝑞𝑥q_{x}italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT with qy=0subscript𝑞𝑦0q_{y}=0italic_q start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT = 0. Data are reported at t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25, for the (nonoptimal) uniform state at doping δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 for Lx=48subscript𝐿𝑥48L_{x}=48italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 48 (squares), and for the optimal striped states at dopings δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (diamonds), 1/5 (red circles), 1/4 (black circles) (for Lx=48,45,72subscript𝐿𝑥484572L_{x}=48,45,72italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 48 , 45 , 72, respectively). The error bars are smaller than the symbol size.

As a reference, we used a uniform state (squares), which is known to be metallic (except at half-filling, when each site is occupied by 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𝛿16\delta=1/6italic_δ = 1 / 6 (diamonds), N(qx,0)/qx𝑁subscript𝑞𝑥0subscript𝑞𝑥N(q_{x},0)/q_{x}italic_N ( italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , 0 ) / italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT clearly tends to zero, compatibly with an insulating behavior. On the other hand, for all the other striped states at δ=1/5𝛿15\delta=1/5italic_δ = 1 / 5 and δ=1/4𝛿14\delta=1/4italic_δ = 1 / 4 (circles), N(qx,0)/qx𝑁subscript𝑞𝑥0subscript𝑞𝑥N(q_{x},0)/q_{x}italic_N ( italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , 0 ) / italic_q start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT 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𝛿16\delta=1/6italic_δ = 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𝛿15\delta=1/5italic_δ = 1 / 5 and δ=1/4𝛿14\delta=1/4italic_δ = 1 / 4(circles), despite having a metallic character, exhibit a suppression in D(r)𝐷𝑟D(r)italic_D ( italic_r ) similar to that of the insulating stripe at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (diamonds). This supports the idea that the stripe order disrupts superconductivity, no matter their metallic or insulating character.

Refer to caption
Figure 4: Pair-pair correlations D(r)𝐷𝑟D(r)italic_D ( italic_r ) as a function of r𝑟ritalic_r on a log-log scale, at t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. Upper panel: Data are reported for the optimal stripe states at different hole-dopings δ𝛿\deltaitalic_δ (circles) and for the (nonoptimal) uniform state at doping δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (blue squares). Lower panel: Data are reported for the uniform state at the critical doping δc=0.29subscript𝛿𝑐0.29\delta_{c}=0.29italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.29 (black squares) along with the (nonoptimal) uniform superconducting state at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (blue squares) and the nonsuperconducting striped state with λ=4𝜆4\lambda=4italic_λ = 4 at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (diamonds).

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 δ𝛿\deltaitalic_δ at which the uniform state becomes energetically favorable and compute the pair-pair correlations. For t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25, as discussed in Table 1, the optimal state at δ=1/4𝛿14\delta=1/4italic_δ = 1 / 4 is a stripe of wavelength λ=3𝜆3\lambda=3italic_λ = 3 while at δ=1/3𝛿13\delta=1/3italic_δ = 1 / 3 we have already reached the uniform state. We then study incommensurate values of δ𝛿\deltaitalic_δ in the range [14,13]1413\left[\frac{1}{4},\frac{1}{3}\right][ divide start_ARG 1 end_ARG start_ARG 4 end_ARG , divide start_ARG 1 end_ARG start_ARG 3 end_ARG ]. Since the wavelengths of the stripes decrease at increasing doping, it is sufficient to compare the striped state with λ=3𝜆3\lambda=3italic_λ = 3 and the uniform state in this doping regime. Their variational energies are presented in Table 2.

Table 2: Energy per site (in units of t𝑡titalic_t) for the best striped state Estripesubscript𝐸stripeE_{\rm{stripe}}italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT with λ=3𝜆3\lambda=3italic_λ = 3 and the uniform state Euniformsubscript𝐸uniformE_{\rm{uniform}}italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT, along with their relative difference ΔE=EstripeEuniformΔ𝐸subscript𝐸stripesubscript𝐸uniform\Delta E=E_{\rm{stripe}}-E_{\rm{uniform}}roman_Δ italic_E = italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT - italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT, as a function of incommensurate δ𝛿\deltaitalic_δ for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25. Data are shown for Lx=48subscript𝐿𝑥48L_{x}=48italic_L start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = 48 for all the stripes and the uniform state. The error bar on the energy is always smaller than 104tsuperscript104𝑡10^{-4}t10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT italic_t.
δ𝛿\deltaitalic_δ Estripesubscript𝐸stripeE_{\rm{stripe}}italic_E start_POSTSUBSCRIPT roman_stripe end_POSTSUBSCRIPT Euniformsubscript𝐸uniformE_{\rm{uniform}}italic_E start_POSTSUBSCRIPT roman_uniform end_POSTSUBSCRIPT ΔEsubscriptΔ𝐸\Delta_{E}roman_Δ start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT
0.26 -0.8799 -0.8784 -0.0015
0.27 -0.8891 -0.8885 -0.0006
0.28 -0.8933 -0.8932 -0.0001
0.29 -0.9015 -0.9016 0.0001
0.31 -0.9086 -0.9087 0.0001

We can identify as the critical doping, the value δc=0.29subscript𝛿𝑐0.29\delta_{c}=0.29italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.29, where the optimal parameters of the stripe state ΔcsubscriptΔ𝑐\Delta_{c}roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and ΔssubscriptΔ𝑠\Delta_{s}roman_Δ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT 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𝛿16\delta=1/6italic_δ = 1 / 6 (blue squares) and the nonsuperconducting striped state at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 with λ=4𝜆4\lambda=4italic_λ = 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, at difference with the striped state.

3.2 t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4

The second set of simulations involved the same search for the optimal state but at a larger value of |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t |, namely t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4. The main effect of a larger value of |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t | is to suppress the stripe pattern and makes the optimal state converge to the uniform one faster. Indeed, while at δ=1/5𝛿15\delta=1/5italic_δ = 1 / 5 the optimal state is a stripe with λ=3𝜆3\lambda=3italic_λ = 3, already at δ=1/4𝛿14\delta=1/4italic_δ = 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|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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 δcsubscript𝛿𝑐\delta_{c}italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT. Following the same reasoning as before, the doping at which the uniform state prevails is for δ𝛿\deltaitalic_δ in the range [15,14]1514\left[\frac{1}{5},\frac{1}{4}\right][ divide start_ARG 1 end_ARG start_ARG 5 end_ARG , divide start_ARG 1 end_ARG start_ARG 4 end_ARG ]. In particular, stripes are suppressed already at δc=0.21subscript𝛿𝑐0.21\delta_{c}=0.21italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.21. In Fig. 5, we show the pair-pair correlations for the case t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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𝛿16\delta=1/6italic_δ = 1 / 6 (empty squares), but still stronger than in the cases where stripe order is established.

Refer to caption
Figure 5: Pair-pair correlations D(r)𝐷𝑟D(r)italic_D ( italic_r ) as a function of r𝑟ritalic_r on a log-log scale. Data are reported at t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4 for the (nonoptimal) uniform state at doping δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (empty squares) and for the optimal striped states at dopings δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (diamonds), 1/5 (circles), as well as for the uniform state at the critical doping δc=0.21subscript𝛿𝑐0.21\delta_{c}=0.21italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.21 (full squares).

To conclude the discussion, we compare the magnitude of the superconducting correlations for the two different values of t/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t. The curves are plotted in Fig. 6. We observe how, when |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t | is larger, superconductivity for the (nonoptimal) uniform state at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 is slightly suppressed. This is equally true for the uniform states at the critical dopings δcsubscript𝛿𝑐\delta_{c}italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT where the uniform state is restored: The effect of a larger value of |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t | in suppressing superconductivity is clearly visible, since superconducting correlations at the critical dopings are smaller at t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4 than at t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25, even if δcsubscript𝛿𝑐\delta_{c}italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is definitely smaller in the first case.

Refer to caption
Figure 6: Comparison of the pair-pair correlations D(r)𝐷𝑟D(r)italic_D ( italic_r ) for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 and t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4, as a function of r𝑟ritalic_r on a log-log scale. Data are reported for the uniform but not optimal superconducting states at δ=1/6𝛿16\delta=1/6italic_δ = 1 / 6 (empty red squares at t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 and empty blue squares at t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4) and for the two uniform states at the critical dopings δc=0.29subscript𝛿𝑐0.29\delta_{c}=0.29italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.29 for t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25 (full red squares) and δc=0.21subscript𝛿𝑐0.21\delta_{c}=0.21italic_δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT = 0.21 for t/t=0.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4 (full blue squares).

4 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/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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)𝐷𝑟D(r)italic_D ( italic_r ) over the last 10 distances is above a threshold value of 3×1043superscript1043\times 10^{-4}3 × 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT  222This threshold value is chosen arbitrarily, but the goal is to show evidences for some residual superconductivity in between the striped states and the uniform nonsuperconducting metal at large doping..

Refer to caption
Figure 7: Phase diagram collecting the optimal states as a function of doping δ𝛿\deltaitalic_δ, for two values of t/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t. The symbol SC denotes superconductivity, while λ𝜆\lambdaitalic_λ indicates the wavelength of the optimal stripe.

By looking for the optimal state for different values of the hole-doping δ𝛿\deltaitalic_δ, 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 δ𝛿\deltaitalic_δ, the wavelength of the stripes shrinks until eventually the uniform state is restored. A larger |t/t|superscript𝑡𝑡|t^{\prime}/t|| italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t | is associated to a shrinking of the wavelength λ𝜆\lambdaitalic_λ and leads to a faster dissolution of the stripes, with the uniform state being the optimal one at a smaller value of δ𝛿\deltaitalic_δ.

The coexistence of superconductivity and stripe order is addressed by looking at the pair-pair superconducting correlations D(r)𝐷𝑟D(r)italic_D ( italic_r ). For both values of t/tsuperscript𝑡𝑡t^{\prime}/titalic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_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 δ𝛿\deltaitalic_δ 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.4superscript𝑡𝑡0.4t^{\prime}/t=-0.4italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.4 all superconducting correlations are weaker than at t/t=0.25superscript𝑡𝑡0.25t^{\prime}/t=-0.25italic_t start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT / italic_t = - 0.25, even if stripes melt at a smaller value of δ𝛿\deltaitalic_δ.

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.

Computational resources were provided by HPC@POLITO (http://www. hpc.polito.it).

References

References