Superconductivity in the Hubbard model: a hidden-order diagnostics from the Luther-Emery phase on ladders

Short-range antiferromagnetic correlations are known to open a spin gap in the repulsive Hubbard model on ladders with $M$ legs, when $M$ is even. We show that the spin gap originates from the formation of correlated pairs of electrons with opposite spin, captured by the hidden ordering of a spin-parity operator. Since both spin gap and parity vanish in the two-dimensional limit, we introduce the fractional generalization of spin parity and prove that it remains finite in the thermodynamic limit. Our results are based upon variational wave functions and Monte Carlo calculations: performing a finite size-scaling analysis with growing $M$, we show that the doping region where the parity is finite coincides with the range in which superconductivity is observed in two spatial dimensions. Our observations support the idea that superconductivity emerges out of spin gapped phases on ladders, driven by a spin-pairing mechanism, in which the ordering is conveniently captured by the finiteness of the fractional spin-parity operator.


Introduction
Since the discovery of high-temperature cuprate superconductors, it became clear that antiferromagnetic correlations are crucial to the onset of superconducting behavior in low dimensions [1]. The phenomenon was observed upon doping an antiferromagnetic insulator, and led to the proposal of a resonating valence bond state of singlet pairs formulated by Anderson [2]. Already in one-dimensional (1D) strongly-correlated materials, a similar superconducting state can be found. Here, the peculiar spin-charge separation allows the opening of the spin gap in the presence of gapless charge excitations, leading to the so-called Luther-Emery (LE) phase. For an appropriate behavior of the charge channel, singlet superconducting correlations may decay to zero slower than any other correlation; in this case the LE phase is often dubbed "superconducting", meaning that superconducting correlations are dominant. Such phase turns out to be characterized by the hidden ordering of a non-local operator, the spin parity operator [3][4][5], which describes the presence of pairs of electrons with opposite spin, that have a finite correlation length. Non-local operators can be measured via quantum gas microscopy in fermionic cold atom systems [6,7].
Because of the characteristic antiferromagnetic behavior of undoped cuprates, the Hubbard model has immediately become the reference model for their investigation. In the following, we will consider the Hubbard model on ladders with M legs and L x rungs, which is defined by: 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 intrachain (t ) and interchain (t ⊥ ) hopping amplitudes are taken to be equal, t = t ⊥ = t; U is the on-site Coulomb interaction. The electronic density is fixed and given by n = N e /L, where N e is the number of electrons (with vanishing total magnetization) and L = M × L x is the total number of sites.
In 1D, the Hubbard model is known to exhibit a LE phase only for attractive interaction, i.e., U < 0 [8]; by contrast, in two spatial dimensions (2D) the Hubbard model is believed to support d-wave superconductivity for repulsive interaction in a large range of doping values [9][10][11][12][13][14][15][16]. A possible precursor of such phase in quasi-one-dimensional lattices could be the spin gapped phase, which is known to characterize the repulsive case on ladders with an even number of legs. In this respect, recent studies [17] confirm that the insulating state at half filling, which is the counterpart of the LE phase in the charge sector, has the same microscopic structure passing from the one to the two dimensional case by increasing the number of legs M of a ladder. A similar behavior has been also observed within the bosonic Hubbard model [18]. The appearance of the insulating phase is signaled by the hidden ordering of a charge-parity operator, which has to be properly normalized to M in order to remain finite up to the 2D limit [18].
For generic values of t and t ⊥ , ladder systems have been an object of intense research due to the fact that they can be assessed both by analytical tools at small values of the Coulomb repulsion, i.e., bosonization and renormalization group, and by numerical approaches, like density-matrix renormalization group (DMRG). In particular, the two-leg ladder has been thoroughly studied by DMRG [19][20][21]. The system is a spin gapped insulator at half filling, with the spin gap persisting also at finite doping for t ⊥ /t < 2. In the case of interest for this paper, i.e., t ⊥ /t = 1, the spin gap closes exactly at quarter filling, in agreement with the U = 0 one-band to two-band transition. In the spin-gapped region, superconducting correlations are present with a power-law decay. On general grounds, it is known that they are dominant when decaying with a power law with exponent smaller than 1. Previous DMRG calculations [19][20][21] suggested that this was not the case, and density-density correlations are dominant; instead, more recently, it was proved by accurate DMRG investigations on large clusters that superconducting correlations dominate in the slightly doped region [22]. As for the presence of the spin gap, similar results have been obtained by the weak-coupling approach of Ref. [23], where the doped two-leg ladder falls into the LE universality class in a large region of the phase diagram. Within a weak-coupling approach, the existence of a spin gap had been already postulated in Refs. [24][25][26]. Finally, the presence of dominating superconducting correlations in the two-leg ladder has been also suggested by bosonization [27]. The three-leg and four-leg case have been also addressed by weak-coupling approaches [28]. The phase diagram is very rich and with a strong dependence on boundary conditions; however some common features may be observed, like the d-wave nature of the electron pairing and the odd-even effect in the spin gap. The latter one means that a spin gap is present at half-filling (and also at finite doping) only when the number of legs is even, as originally observed for spin models [29,30] and later derived with bosonization techniques for the Hubbard model at weak coupling [31]. Finally, the six-leg case has been investigated by DMRG and by other numerical approaches in order to assess the stability of striped phases with charge and spin modulations in the strong coupling limit [32][33][34].
In this paper, by extrapolating to the 2D limit the variational Monte Carlo results for the fractional spin parity on two-, four-, and six-leg ladders, we provide evidence that the superconducting phase of the 2D repulsive Hubbard model derives from the corresponding LE phase with dominant superconducting correlations on ladders with M even legs. In particular, as resumed in the main panel of Fig. 1, we show that the density range where both spin gap and spin parity are finite on ladders extrapolates to the range in which superconductivity is observed in the two-dimensional model.

Spin parity
In 1D, nonlocal order parameters play a fundamental role in identifying gapped phases of interacting fermions, when spontaneous symmetry breaking by a local order operator, like magnetization, is absent. In the presence of spin-charge separation, it has been shown that each spin/charge gapped phase of correlated electrons has a hidden long-range order, i.e., a finite value of a spin/charge parity or string operator [3,4]. In particular, a finite spin gap signals the emergence of hidden order in the spin channel, which in the LE liquid phase is detected by the finiteness of the expectation value of a spin parity operator O 1D This describes a liquid of holons and doublons in which single electrons appear as bound states (with opposite spins).
Since such state has no reason to disappear when moving to ladders, it is desirable to generalize the definition of spin parity to a higher dimensional geometry. The natural extension to a ladder would be a brane of parities, that includes a further summation over the sites of the rung. Similarly to what has been proposed for describing the Mott phase on ladders in the charge channel [17,18], we can define: where ∆S(M ) = M y=1 S z x,y is the spin fluctuation on the brane with fixed x (here, it is necessary to indicate explicitly the x and y components of the site R in the label of the spin operator, i.e., S z R ≡ S z x,y ). As in the 1D case, the asymptotic limit with L x → ∞ must be considered. The quantity defined in Eq. (2) equals 1 in case no single electrons are present and is also unaffected by those electrons which form pairs with opposite spin within the brane; instead, it changes sign whenever one of such pairs crosses the boundary. The longer the length of the boundary (in our case 2M ) the larger is the number of such defects which destroy the presence of order. In particular, in the 2D limit M → ∞ (in addition to L x → ∞), the expectation value of O s (M, L x ) would become vanishing, despite the persisting order. Indeed, within a Gaussian approximation one has: where . . . indicates the expectation value over the ground-state wave function. In this case, we have that lim M →∞ [∆S(M )] 2 = ∞. The problem can be solved upon properly renormalizing the prefactor at the exponent of Eq. (2) with the number of legs, and introducing the fractional spin parity order as: We mention that, with respect to Ref. [18], for simplicity here we only consider the case in which the exponent of M is equal to 1. For convenience, we also denote as the L x → ∞ limit of the fractional spin parity. With the definitions of Eqs. (4) and (5), we obtain that the fractional spin parity remains finite also in its 2D extrapolation C s (∞) within the spin gapped phase. Notice that the vanishing of C s (∞) in the disordered, i.e., spin gapless, phase is not affected by the fractional definition of the spin parity [18,35]. To this end, we emphasize that, when extrapolating within the gapless phase, the limit for L x → ∞ has to be performed before the limit M → ∞.
In summary, the fractional parity remains finite also in the 2D limit whenever the pairs of electrons with opposite spin have finite correlation length. As anticipated in the inset of Fig. 1, we will see that this latter feature is connected to the presence of d-wave superconductivity.

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 [36]. In particular, we consider the so-called Jastrow-Slater wave functions that extend the original formulation proposed by Gutzwiller to include correlations effects on top of uncorrelated states [37,38]. Our variational wave functions are described by: where |Φ 0 is an uncorrelated state that corresponds to the ground state of the following uncorrelated Hamiltonian [39,40]: which includes a free-band dispersion ξ k , defined as in the Hubbard Hamiltonian of Eq. (1), and a BCS coupling ∆ k , with pairings ∆ x and ∆ y along the x and the y direction, respectively. The parameters ∆ x and ∆ y are optimized to minimize the variational energy (while t = 1 sets the energy scale of the uncorrelated Hamiltonian). The presence of the BCS pairings allows us to generate a finite spin gap in the variational state on ladder systems, as shown in the next section, while it leads to superconductivity in the 2D case [13,16,41]. The effects of correlations are introduced by means of the so-called Jastrow factor J [42,43]: 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 |. In particular, the Jastrow factors are crucial to describe the Mott insulator at half filling. While the model is insulating only at half filling and metallic at any finite doping, a spin gap is present in a large region of doping close to half filling. In order to assess the gapped or gapless nature of the spin excitations, we calculate the spin-spin structure factor S(q), defined as: where . . . indicates the expectation value over the variational wave function. In analogy with the relation between the density structure factor and the charge gap [44,45], we have that a spin gap ∆ s is present whenever S(q) ∝ |q| 2 for |q| → 0, while gapless spin excitations are present for S(q) ∝ |q|, since ∆ s ∝ lim q→0 |q| 2 /S(q). In the following, we consider the quantity: to detect the presence/absence of the spin gap. 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). As before, here we explicitly denoted the coordinates of the site R, i.e., c † R,σ ≡ c † x,y,σ . Our simulations are performed with periodic boundary conditions along the x direction, while, in the y direction, they are taken to be open for M = 2, antiperiodic for M = 4, and periodic for M = 6.

Results
First of all, we report in Fig. 2 the optimal value of the variational BCS parameters ∆ x and ∆ y , as a function of the electron density n. The symmetry of the two parameters resembles the d-wave one, being ∆ y −∆ x , as expected for a ladder system that becomes a square lattice when the number of legs equals the number of rungs.
Then, we investigate the presence of a spin gap, by looking at the behavior at small momenta of the spin-spin correlations defined in Eq. (9). As an example, we plot in Fig. 3 the quantity q 2 x /S(q x ), as a function of q x at q y = 0, for the six-leg case. If this quantity extrapolates to a finite number, the system has a spin gap, otherwise it is gapless. Then, in Fig. 4, we report the extrapolation of q 2 x /S(q x ) to the q x = 0 limit, denoted as E s , for ladders with two, four, and six legs. The lattice sizes are large enough to be at the thermodynamical limit. We remark that the quantity E s is only proportional to the spin gap, thus not providing a quantitative estimate of it. In particular, the proportionality constant depends on the geometry of the lattice, thus not allowing for a direct comparison of the gap size between lattices with a different number of legs. Nonetheless, our results show that the gap is finite in the region where the variational parameters ∆ x and ∆ y are finite, while it vanishes when the pairing terms are not present in the optimal variational state. Indeed, a finite spin gap is associated to a BCS gap with no gapless points in the wave function. Given the optimal dwave nature of the pairing, this condition is satisfied with the boundary conditions described in the previous paragraph, that do not include the points (±π/2, ±π/2) in reciprocal space. Finally, note that while the spin gap decreases monotonically to zero for increasing doping in the four-leg and six-leg systems, it is non monotonic in the two-leg case.
In Fig. 5, we plot the pair-pair correlations D(x) defined in Eq. (11) for two-, four-, and sixleg ladders at n = 0.875. Data are shown on a log-log-scale in order to highlight the power-law decay typical of ladder systems: D(x) ∝ x −γ . The exponent γ changes from approximately 1.45 for M = 2 to approximately 0.34 for M = 6. Even if a precise determination of the exponent γ is hard, as highlighted by DMRG in the two-leg case [21,22], nonetheless our results indicate that pair-pair correlations become stronger and stronger as M increases, in agreement with a finite value of lim x→∞ D(x), that is expected in the 2D limit [16].
Let us now consider the fractional spin parity of Eq. (4). Our results, shown in Fig. 6, indicate that the fractional spin parity C s (M, L x ) increases when approaching half filling, following the behavior of the spin gap, with a value that does not depend on L x . On the contrary, a flat behavior in C s (M, L x ) occurs in the region where no spin gap is present, with this offset getting smaller and smaller when L x → ∞. The gray areas in Fig. 6 indicate the doping region where the behavior of the fractional spin parity changes and are in agreement with the opening of the spin gap shown in Fig. 4. In Fig. 7, we show, for the two-leg case, that the almost constant value of C s (M, L x ) in the spin gapless region indeed scales to zero at increasing linear size L x , with the power-law decay illustrated in Ref. [18]. We would like to point out that, even if there are strong finite size effects in the fractional spin parity, one can still clearly distinguish between two different regimes: A flat and size dependent behavior of C s (M, L x ) in the spin gapless region and a value of the fractional spin parity that is proportional to the spin gap and that does not depend on the lattice size in the spin gapped region. Our results indicate then that the fractional spin parity can be considered a good indicator for the presence of a spin gap in ladder systems. Moreover, we observe that C s (M, L x ) is approximately equal in the spin gapped regions of the four-leg and of the six-leg systems, supporting the fact that this quantity remains finite in the 2D limit.
We would like to mention that our picture is still valid when the variational wave function includes the stripe order of Ref. [34]. Indeed, we have optimized at doping 1/8 (i.e., n = 0.875), a variational state that combines BCS pairing with stripe order, on different lattice sizes [46]. Even if in such a state both the variational BCS parameters and the fractional spin parity are reduced with respect to the uniform case, they are size independent, thus indicating that the fractional spin parity remains finite in the thermodynamical limit.
The previous results are resumed in the phase diagram anticipated in Fig. 1 fractional parity C s (M ) is finite, is reported as a function of the number of legs. Our results show that the density at which the spin gap opens becomes closer to half filling as the number of legs increases, approaching the value where superconductivity develops in the 2D case (n ≈ 0.75) [16]. Moreover, in the inset of Fig. 1, we present a tentative extrapolation of the fractional spin parity C s (M ) of Eq. (5) to the 2D case. In the shaded region with n 0.75, i.e., where superconductivity is present, C s (∞) is extrapolated from the results for the two-, four-, and six-leg cases, that are shown in Fig. 6. In the region with n 0.75, we expect that the fractional spin parity vanishes in the thermodynamical limit, according to the results presented in Figs. 6 and 7.

Conclusions
By means of a variational quantum Monte Carlo analysis of the repulsive Hubbard model, we have investigated the possibility that at the origin of both the LE spin gapped phase with dominant superconducting correlations, observed on ladders with even number M of legs, and the 2D d-wave superconducting region, there is the same mechanism of formation of pairs of electrons with opposite spin, that have a finite correlation length. This aspect is captured by the hidden ordering of a spin parity operator, which remains finite also in the 2D limit upon an appropriate normalization to M .
Our results show a good agreement between the presence of a finite spin gap and a finite value of the fractional spin parity. In particular, the value of the fractional spin parity is size independent and proportional to the size of the gap in the spin gapped phase, while it becomes a size dependent constant in the gapless phase, extrapolating to zero as a power law in the linear dimension of the system. The spin gapped region is characterized by a finite value of the BCS pairing terms in the variational state, with d-wave symmetry, and by superconducting pair-pair correlations with a power-law decay. The region where the spin gap (and the fractional spin parity) is finite shrinks with the number of legs, with an extrapolation to the 2D case that matches with the region where superconductivity is detected in the 2D Hubbard model, i.e., where pair-pair correlations remain finite at large distances.
It must be stressed that the nature of the LE superconducting phase changes when passing from the attractive to the repulsive case discussed here. In fact, in the attractive case, the pairs of correlated single electrons represent minority fluctuations with respect to the majority of uncorrelated holons and doublons. Whereas in the repulsive case the electron pairs with up-down spin are themselves the majority of electrons. The difference is seen for instance in the dependence on the number of legs of the average squared spin fluctuations on the brane, which according to Eqs. (3) and (4) appears to be proportional to M 2 , at variance with the M dependence typical of the U < 0 case. In addition, while in the U < 0 case the superconducting pairs are on site, for U > 0 they are on neighboring sites. At half-filling, they coexist with the minority of localized holon-doublon pairs in the insulating phase. Upon doping, the holon-doublon pairs break and the electron pairs with up-down spin delocalize. In the 2D limit, this originates a collective behavior characterized by spontaneous symmetry breaking and local order, as the finiteness of pairing correlation proves. Thus, we provide an example of how hidden orders may evolve into standard local orders in higher dimension.