Metallic and insulating stripes and their relation with superconductivity in the doped Hubbard model

The dualism between superconductivity and charge/spin modulations (the so-called stripes) dominates the phase diagram of many strongly-correlated systems. A prominent example is given by the Hubbard model, where these phases compete and possibly coexist in a wide regime of electron dopings for both weak and strong couplings. Here, we investigate this antagonism within a variational approach that is based upon Jastrow-Slater wave functions, including backflow correlations, which can be treated within a quantum Monte Carlo procedure. We focus on clusters having a ladder geometry with $M$ legs (with $M$ ranging from $2$ to $10$) and a relatively large number of rungs, thus allowing us a detailed analysis in terms of the stripe length. We find that stripe order with periodicity $\lambda=8$ in the charge and $2\lambda=16$ in the spin can be stabilized at doping $\delta=1/8$. Here, there are no sizable superconducting correlations and the ground state has an insulating character. A similar situation, with $\lambda=6$ appears at $\delta=1/6$. Instead, for smaller values of dopings, stripes can be still stabilized, but they are weakly metallic at $\delta=1/12$ and metallic with strong superconducting correlations at $\delta=1/10$, as well as for intermediate (incommensurate) dopings. Remarkably, we observe that spin modulation plays a major role in stripe formation, since it is crucial to obtain a stable striped state upon optimization. The relevance of our calculations for previous density-matrix renormalization group results and for the two-dimensional case are also discussed.


Introduction
Despite a long time since the first experimental evidence of superconductivity in copper-oxide materials [1], the phase diagram of high-temperature cuprate superconductors still contains many unsolved questions. One of them is the possible coexistence or competition of superconductivity with inhomogeneous orders in the charge and spin sectors [2][3][4]. Starting from the neutron scattering observation of half-filled stripes in La 1.48 Nd 0.4 Sr 0.12 CuO 4 around doping 1/8 [5], different inhomogeneous orders have been proposed with a variety of experimental probes, ranging from scanning tunneling microscopy to x-ray scattering [6][7][8][9].
From a theoretical point of view, striped states were first considered in the strong-coupling limit, i.e., within the so-called t − J model. Here, contradictory results have been obtained with different numerical methods. In particular, half-filled stripes are found to be stable by the density-matrix renormalization group (DMRG) [10], the infinite projected entangled-pair states (iPEPS) [11], and variational Monte Carlo [12] approaches, while the absence of stripes has been reported by improvements of variational Monte Carlo, including the application of a few Lanczos steps and an imaginary-time projection performed within a fixed-node approximation [13]. Even though the t − J model has been widely used in the past to describe the low-energy physics of cuprate materials, neglecting local charge fluctuations is not completely justified; in addition, antiferromagnetic correlations are usually overestimated, because of the presence of a direct super-exchange term. These two ingredients may have an important impact on the possible stabilization of charge and/or spin modulations at finite dopings. In this respect, a better description can be obtained by considering the Hubbard model: where c † R,σ (c R,σ ) creates (destroys) an electron with spin σ on site R and n R,σ = c † R,σ c R,σ is the electronic density per spin σ on site R. In the following, we indicate the coordinates of the sites with R = (x, y). The nearest-neighbor hopping integral is denoted as t, while U is the on-site Coulomb interaction. The electronic doping is given by δ = 1 − N e /L, where N e is the number of electrons (with vanishing total magnetization) and L is the total number of sites. For ladder systems, which will be considered in the paper, we define M as the number of legs and L x as the number of rungs, the total number of sites being given by L = M × L x .
For positive U and an even number of legs M , the formation of short-ranged antiferromagnetic electron pairs is responsible for the opening of a spin gap, with dominant superconducting correlations for M ≥ 4 [14] in a finite doping range 0 < δ ≤ δ c . This is signaled by the hidden ordering of a spin-parity operator [15,16] and is ultimately responsible for the d-wave superconductivity observed in the uniform state of the two-dimensional Hubbard model [17][18][19][20][21][22][23][24]. Such gapped spin phase is possibly accompanied by the formation of stripes, which may suppress the superconducting order. Indeed, early DMRG studies suggested the formation of stripes on the 6-leg Hubbard model [25,26]; these results have been corroborated by using constrained path auxiliary-field Monte Carlo [27] and density-matrix embedding theory [28]. More recently, a combination of advanced numerical methods [29] focused on doping 1/8, where maximal inhomogeneity is observed in many high-temperature superconductors. This work highlighted the fact that the ground state is a bond-centered linear stripe, with periodicity λ = 8 in the charge (which defines the wavelength of the stripe) and periodicity 2λ = 16 in the spin, due to the so-called π-phase shift; in this case, the stripe is filled and commensurate with the doping, since it accommodates two holes in a one-dimensional cell of length 16. This striped state is reported to have an energy lower of about 0.01t with respect to the uniform state. Other states, such as diagonal stripes or checkerboard order have been tested, but they were not found to be competitive with the linear stripe. Even if different stripe wavelengths are close in energy, the most favorable one is the filled one with wavelength λ = 8, in agreement with older Hartree-Fock calculations [30][31][32][33], but unlike the stripes in real materials, that are half-filled with wavelength λ = 4. The presence of nearest-neighbor d−wave coupling coexisting with the stripe has been also reported in this study [29]. A striped nature of the ground state, also away from the 1/8 doping, has been suggested by both cellular dynamical mean-field theory [34] and variational Monte Carlo calculations [35,36]. Both methods report the presence of stripe in the doping range 0.07 δ 0.17, with the ground state being uniform for higher dopings; the stripe wavelength is shown to decrease when the doping increases, with the optimal stripe wavelength at doping 1/8 being slightly different from what reported in Ref. [29]. Interestingly, a finite superconducting order parameter at large distances has been reported to coexist with stripes in Ref. [36]. Concerning the wavelength of the stripe, it has been suggested by a finite-temperature determinant quantum Monte Carlo study that the presence of a finite next-nearest-neighbor hopping t reduces the stripe wavelength, from λ = 8 at t = 0 to λ = 5 at t /t = −0.25 [37].
In this paper, we employ the variational Monte Carlo method with backflow correlations to study the presence of striped ground states in the doped Hubbard model. First, we consider doping δ = 1/8 on ladders with M = 2, 6, and 10, in order to extrapolate to the twodimensional limit. While no evidence of static stripes is seen on the 2-leg ladder, a filled stripe of wavelength λ = 8 is stabilized on a larger number of legs, as well in the two-dimensional extrapolation. This state is compatible with an insulator, with charge and spin order and no superconductivity. On a 6-leg ladder, we have investigated also the dopings away from δ = 1/8. In particular, at doping δ = 1/6, a striped state of wavelength λ = 6, compatible with an insulating state, occurs, while, at doping δ = 1/12, the best striped state is not filled, having wavelength λ = 8. However, this state is only weakly metallic, since, on a 6-leg ladder, it is commensurate with the doping, having an even number of holes per spin modulation. On the contrary, away from these commensurate dopings, stripes can be metallic and coexisting with superconductivity, even if electron pairing is suppressed with respect to the homogeneous case. In this respect, we show that in the whole range between doping δ = 1/12 and δ = 1/8, the best ground-state approximation is a striped state with wavelength λ = 8, but with a metallic behavior and superconducting correlations with a power-law decay. Furthermore, we report that spin modulation plays a major role in stripe formation, since it is crucial to obtain a stable striped state upon optimization. The role of spin modulation can be also observed in the formation of a peak in spin-spin correlations, that is much stronger than the one in density-density correlations. In summary, our results show that striped states are a common feature of the doped repulsive Hubbard model. They coexist with superconductivity away from selected dopings where they are commensurate, their formation being favored by the presence of magnetic order with π shift.

Variational Monte Carlo method
Our numerical results are obtained by means of the variational Monte Carlo method, which is based on the definition of suitable wave functions to approximate the ground-state properties beyond perturbative approaches [38]. In particular, we consider the so-called Jastrow-Slater wave functions that include long-range electron-electron correlations, via the Jastrow factor [39,40], on top of uncorrelated states, extending the original formulation proposed by Gutzwiller [41,42]. In addition, the so-called backflow correlations will be applied to the Slater determinant [43,44], in order to sizably improve the quality of the variational state. Our variational wave functions are described by: where J is the Jastrow factor and |Φ 0 is a state that, starting from an uncorrelated wave function obtained from an auxiliary Hamiltonian, redefines the orbitals on the basis of the many-body electronic configuration, incorporating virtual hopping processes, as discussed in Refs. [43,44]. Unless otherwise specified, results are obtained considering backflow corrections. We consider two different kinds of wave functions, the first one that is appropriate to study homogeneous superconducting states (with possible additional Néel antiferromagnetism) and the second one that describes bond-centered striped states with charge and spin modulation. The latter one can be further supplemented with either uniform or modulated pairing.
Let us start with the homogeneous superconducting state that is the ground state of the following auxiliary Hamiltonian: which includes a free-band dispersion k , defined as in the Hubbard Hamiltonian of Eq. (1), a BCS coupling ∆ k , with nearest-neighbor pairings ∆ x and ∆ y along the x and the y direction, respectively, and a chemical potential µ. Néel antiferromagnetism can be further included by All the parameters ∆ x , ∆ y , ∆ AF , and µ are optimized to minimize the variational energy (while t = 1 in the definition of k , sets the energy scale of the uncorrelated Hamiltonian). Bond-centered striped states can be included in the variational wave function by means of charge and spin modulations along the x direction. The auxiliary Hamiltonian reads as Here, the term describes a charge modulation in the x direction with periodicity λ = 2π/Q; the term describes an antiferromagnetic order that is the standard Néel one along the y direction, while it is modulated along the x direction with a periodicity that is doubled with respect to the charge one, i.e., 2λ = 4π/Q; indeed, after one period in the charge modulation, spins are reversed, as illustrated in Fig. 1 for λ = 8. In the presence of stripes, ∆ c and ∆ AF are further variational parameters to be optimized.  The BCS term that is included in Eq. (5) can be either uniform with pairings ∆ x and ∆ y along the x and the y direction, respectively, or modulated with the following periodicity, that has been named "in phase" in Ref. [12]: wherex andŷ are the versors along the x and the y directions, respectively, while ∆ x and ∆ y are variational parameters to be optimized. In Ref. [12], a different modulation in the pairing, named "antiphase" has been introduced, without the absolute values that are present in Eq. (8). However, in our simulation, we always found that the "antiphase" modulation does not lead to any energy gain with respect to the "in phase" case.
The Jastrow factor J is defined as: where n R = σ n R,σ is the electron density on site R and v R,R (that include also the local Gutzwiller term for R = R ) are pseudopotentials that are optimized for every independent distance |R − R |. The presence of charge disproportionations can be analyzed by means of the static structure factor N (q) defined as: where . . . indicates the expectation value over the variational wave function. In particular, the presence of a peak at a given q vector denotes the presence of charge order in the system. Moreover, the static structure factor allows us to assess the metallic or insulating nature of the ground state. 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 [44,45]. Analogously, spin order in the system is associated to a peak in the spin-spin correlations defined as: where S z R is the spin operator along the z direction, i.e., S z R = 1/2(c † R,↑ c R,↑ − c † R,↓ c R,↓ ). Finally, pair-pair correlations can be computed by calculating a correlation function between singlets on rungs at distance x, defined as where is a vertical singlet located on the rung between sites of coordinates (x, 1) and (x, 2). Here, we explicitly denoted the coordinates of the site R, i.e., c † R,σ ≡ c † x,y,σ .

Results
Let us focus on the possibility that a bond-centered stripe may be stabilized at dopings δ = 1/p, with p even. These cases can, in principle, accommodate a filled stripe of wavelength λ = p (the full periodicity of the stripe is 2p due to spin), as recently suggested in Ref. [29], where the case with δ = 1/8 has been considered in detail. In order to fit the stripe with our cluster with L = M × L x sites, we take L x to be multiple of 2p. First of all, in Fig. 2, we plot the energy gain ∆E = E stripe − E homogen between a homogeneous superconducting state and a striped one for the 2-, 6-, and 10-leg ladders. In order to properly scale to the two-dimensional case, we preserve the aspect ratio of the ladder considering L = M × (8 × M ). We observe that the energy gain of the striped state with respect to the uniform one increases linearly with the system size. We mention that the less accurate variational wave function without backflow correlations has a larger energy gain with respect to the case in which these correlations are included. Still, even in the most accurate calculations, there is a finite (and relatively large) energy gain when charge and spin modulations are considered. Even though a reliable extrapolation to the thermodynamic limit is not easy to perform, we verified that a striped state can be stabilized also in a square-lattice geometry, as we checked for a L = 16 × 16 lattice, observing an energy gain of about 5 × 10 −3 . In Fig. 3, we report the density-density correlations N (q) and the spin-spin correlations S(q) at doping δ = 1/8, for both the 2-leg and the 6-leg case, comparing the results for a homogeneous wave function and for a striped state. For the 2-leg case, the correlation functions are the same for both wave functions, in agreement with the energy results of Fig. 2, where it is shown that almost no energy gain comes from the explicit inclusion of stripes in the variational state. On the contrary, for the 6-leg case, the striped wave function is characterized by a peak in N (q) at Q = (π/4, 0), in agreement with a periodicity λ = 8 in the x direction and no charge modulation along the y direction. A peak appears also in the spin-spin correlations S(q) at a vector Q = (7π/8, π). Indeed, the spin modulation described in Eq. (7) is antiferromagnetic along the y direction, while it is antiferromagnetic, with an additional modulation of periodicity 2λ = 16, along the x direction. We remark that the peak in S(q) is much stronger than the one in N (q), suggesting that spin modulation plays a major role with respect to charge modulation. We would like to emphasize that, by setting ∆ c = 0 in Eq. (6) and optimizing only ∆ AF , we recover the striped state, with peaks in both S(q) and N (q), while, setting ∆ AF = 0 in Eq. (7) and optimizing only ∆ c , the stripe is not stable and the optimized state converges to the homogeneous one. This fact implies that spin modulation plays a major role in stripe formation. As suggested in Ref. [29], our results show that the filled striped state with periodicity λ = 8 in the charge sector is the appropriate one for the Hubbard model at doping δ = 1/8. The half-filled striped state with periodicity λ = 4, which has been previously suggested for the t − J model [10,12], is not stable upon optimization in the Hubbard model. This observation remains valid also for a larger value of the Coulomb repulsion, i.e., U/t = 16. In addition to charge and spin modulation, our variational state can also include pairing, either homogeneous or with the modulation of Eq. (8). Both states can be stabilized in the variational state even if neither of them leads to an energy gain with respect to the striped state with no pairing that is shown in Fig. 2. We have then compared homogeneous and striped wave functions at dopings δ = 1/p, with p = 12, 10, and 6, that are also relevant for the phase diagram of cuprate superconductors. We observe that at dopings δ = 1/12 and δ = 1/10, the best striped state has wavelength λ = 8, like in the δ = 1/8 case. Striped states of wavelength λ = 12 and λ = 10 can be stabilized as local minima, with an energy difference with respect to the best stripe of the order of 10 −3 . On the contrary, the striped state of wavelength λ = 8 cannot be stabilized at doping δ = 1/6, where, instead, a stripe of wavelength λ = 6 characterizes the ground state. In Table 1, we present a summary of the energies of homogeneous and striped states. Our results show that the energy gain due to striped phases with respect to homogeneous phases is Table 1: Energies for the homogeneous state, the best striped state, as well the energy gain ∆E = E stripe − E homogen between the best striped state and the homogeneous one, at dopings 1/p with p = 6, 8, 10, and 12. Within the error bar, results are independent from the lattice size (L = 6 × L x with L x = 48, 72, and 96 at p = 6; L = 6 × L x with L x = 48, 64, and 96 at p = 8, 12; L = 6 × 80 at p = 10).  (1) becoming larger when approaching half filling. Furthermore, we observe that at all densities we can stabilize a finite BCS pairing, both uniform and modulated, within the striped state, even if this pairing is reduced with respect to the uniform case. The insulating or metallic nature of the ground state can be seen in the small-q behavior of the static structure factor N (q). In Fig. 4, we show N (q)/q x for homogeneous states (empty symbols) and for the striped ground states (full symbols) at doping 1/p, with p = 12, 10, 8, and 6. We observe that in the homogeneous cases, as well as for the optimal striped state at p = 10 (with λ = 8), N (q) q x at small q x , clearly indicating that the state is metallic. The latter result can be explained by the fact that the optimal stripe is not commensurate with the doping, since an even number of holes cannot be accommodated in each of the L x /16 periods in the spin modulation. The results for the striped ground states at p = 8 and 6 are instead compatible with an insulating ground state, with the behavior of N (q)/q x extrapolating to zero when the lattice size increases. These striped states accommodate 2M holes in each of the L x /2p periods in the spin modulation. At p = 12, the situation is less clear, but more compatible with a metallic behavior for the optimal striped state with λ = 8 and with an insulating behavior for the filled stripe with λ = 12. In this case, even if the optimal striped state is not filled, it is still commensurate with the doping since, on a 6-leg ladder, an even number of holes can be accommodated in each of the L x /16 periods in the spin modulation (with an unequal distribution of holes in different legs). The ordered nature of these striped states can be clearly seen in the emergence of peaks in the charge and in the spin structure factors. Given the wavelength λ of the optimal stripe, N (q) has a peak at Q = (2π/λ, 0) and S(q) has a peak at Q = (π(1 − 1/λ), π), as reported in Fig. 5 for p = 12, 8, and 6. The divergence of these peaks with the system size is also reported in Fig. 6, with both lim L→∞ N (Q)/L and lim L→∞ S(Q)/L decreasing upon increasing doping.
We mention that the stripe wavelength increases when approaching half filling, as we have verified at doping δ = 1/16, where the optimal striped state is metallic with λ = 12. Here, different stripe wavelengths with a long periodicity are close in energy, as expected when the system is prone to phase separation. Indeed, at the accuracy level of variational Monte Carlo, a region of phase separation is predicted close to half filling [36,46,47].
Generic dopings are more difficult to treat, because they cannot be exactly reproduced on different lattice sizes; however, some predictions can be formulated also in this case. In particular, for the density range between dopings δ = 1/12 and δ = 1/8, where we can assume that a possible striped state has to have a wavelength λ = 8, since this is the optimal one for the two extremal cases. Given this assumption, we can consider a generic incommensurate doping: for example, doping δ = 0.104 that can be approximately obtained on both L = 6×48 and L × 64 lattice sizes. Here, the ground state is metallic in analogy with the p = 10 case, with an energy gain of the order of 6 × 10 −3 when stripes are included. A similar analysis can be performed also for a generic doping, in the range between δ = 1/8 and δ = 1/6, leading to analogous conclusions on the metallic nature of the ground state. In this case, we cannot exclude that, for a sufficiently large lattice size, a suitable striped state leading to an insulating ground state may be found, since the optimal stripe wavelength shifts from λ = 8 to λ = 6. Finally, we observe that for dopings δ 0.20, no stripe order can be stabilized in the wave function, with the ground state being homogeneous up to δ c 0.27, as in Ref. [14].
Finally, we show in Fig. 7 the superconducting correlations D(x). From the analysis of the results, we deduce that superconductivity, in the presence of stripes, disappears for all the three dopings δ = 1/p with p = 12, 8, and 6, with D(x) 0 for large enough x. This result is in agreement with the insulating (or weakly metallic) nature of the striped states. Once we move away from these dopings, the system is metallic and superconducting correlations are present, even if suppressed with respect to the homogeneous case. In particular, in Fig. 7, we present results for p = 10 and for two generic dopings, δ = 0.104 and δ = 0.146. Our results in the presence of stripes are shown with uniform superconductivity, but they would be very similar also for the modulated one of Eq. (8).

Conclusions
In this work, we studied the possibility to stabilize charge and spin modulations (stripes) in the Hubbard model in a non vanishing doping range, and their coexistence with superconductivity. To this end, we employed variational wave functions that contain both Jastrow and backflow correlations. Even though the calculations have been performed on ladder geometries with L = M × L x sites (with M = 2, 6, and 10, and L x M ), the results are expected to be relevant also for the thermodynamic limit. Indeed, as discussed in Ref. [29], the DMRG calculations performed on ladders with M = 4 and 6 already contain the qualitatively correct features of the two-dimensional limit. In addition, our results on the 16 × 16 cluster (not shown) are compatibles with the ones obtained on ladders.
The main results are summarized in Fig. 8. At doping δ = 1/8 and 1/6, we obtain insulating stripes of wavelength λ = 1/δ (the full periodicity is doubled due to spin). The stripe with λ = 8 appears to be particularly stable, since it corresponds to the lowest-energy state also for δ = 1/10 and δ = 1/12. In the former case, the ground state shows a metallic behavior, which is compatible with the mismatch between the doping level and the periodicity of the stripe; furthermore, we have verified that also for the 10-leg case the best striped state at δ = 1/10 is metallic with a wavelength of λ = 8. In the latter case, instead, it is not so clear whether the ground state is metallic or insulating, since N (q)/q x may extrapolate either to zero (insulating) or to a small finite value (metallic); moreover, superconducting correlations decay rapidly to zero with the distance x. A weakly metallic behavior with no superconducting correlations may be related to the proximity to the phase separation region that, within the variational Monte Carlo accuracy, is present close to half filling. Intermediate (incommensurate) dopings are much more difficult to assess; however, we obtain that a metallic stripe with λ = 8, coexisting with superconductivity, can be stabilized in the whole doping range between δ = 1/12 and δ = 1/8. Also for δ > 1/8 we have indications of superconducting stripes, even if we cannot exclude that, for a sufficiently large lattice size, a suitable striped state leading to an insulating ground state may be found. Finally, for 0.20 δ δ c = 0.27, uniform superconductivity is present. Our results indicate that the best place to observe a suppression of superconductivity due to stripes is indeed 1/8 doping. In addition, we foresee that sizable effects should be also visible at δ = 1/6.
One important aspect of our calculations is that stripe order is driven by spin and not by charge. Indeed, as discussed in the text, if we only include modulations in charge, but not in spin, no stripe order can be stabilized in the variational optimization; instead, a striped ground state is reached when only spin modulations are included in the variational state.