A simple case of d(x-y) pairing: Hubbard ladder

We study the strength and the temperature scale of the d(x-y) pairing correlations in the Hubbard model on a ladder lattice using Quantum Monte Carlo (QMC) simulations. In particular, we present QMC results on the particle-particle interaction and the solution of the Bethe-Salpeter equation for the d(x-y)-wave BCS channel. These data show that there are strong d(x-y) pairing correlations in the Hubbard ladder for certain values of the model parameters.


Introduction
The pairing correlations in the Hubbard ladder have been studied in detail using various techniques of many-body physics [1][2][3][4][5][6]. In particular, the DMRG calculations found that the d(x 2 -y 2 ) pairing correlations are strongest near half-filling for Coulomb repulsion in the intermediate coupling regime and for hopping anisotropy t ┴ /t≈1. 5. In fact, the two-leg Hubbard ladder is probably the only model where it is known from exact calculations that the pairing correlations get enhanced by turning on an onsite Coulomb repulsion in the ground state. Hence, it is important to know the characteristic temperature scale and the strength of the pairing correlations in this model. Here, we present numerical results with this purpose.
In the following, we will study the momentum dependence of the reducible particle-particle interaction Γ in the BCS channel, and show QMC results on the evolution of Γ with temperature T, doping and the Coulomb repulsion U. In order to investigate the strength of the pairing corelations in this system, we will present the solution of the Bethe-Salpeter equation in the BCS channel. These QMC results show that the two-leg Hubbard ladder exhibits strong d(x 2 -y 2 ) type pairing correlations for certain values of the model parameters.
The DMRG calculations have found that the pairing correlations are most enhanced in the intermediate coupling regime and near half-filling for values of t ┴ /t for which the top of the bonding band and the bottom of the antibonding band are near the Fermi level [4]. In this case, the QMC data show that, as the temperature decreases, a large repulsive peak develops in Γ for scattering processes involving q=(π,π) momentum transfer. In addition, Γ becomes strongly attractive for scatterings with q≈0 momentum transfer when the incoming and the outgoing particles are on the Fermi surface. It will be seen that this is due to resonant scattering in the d(x 2 -y 2 ) pairing channel. These QMC data also show that Γ depends not only on the momentum transfer but also on the momentum of both the incoming and the outgoing fermions. For instance, the q=0 resonant scattering is much weaker at the saddle points (0, ±π) and (±π,0). In addition, we will see that at half-filling the q=0 resonant scattering is significantly reduced with respect to the q=(π,π) scattering. The QMC calculations find that, at low temperatures, the leading eigenvalue of the Bethe-Salpeter equation has an eigenfunction which has d(x 2 -y 2 ) type symmetry in the sense that it changes sign as one goes from (π,0) to (0,π) in the Brillouin zone. In addition, the temperature evolution of the d(x 2 -y 2 ) eigenvalue follows the development of the AF correlations. Here, we will see that the d(x 2 -y 2 ) eigenvalue of the Bethe-Salpeter equation becomes as large as 0.75 at T=0.1t for U=4t and 6% doping, while Γ attains values which are an order of magnitude larger than the bare Coulomb repulsion. The results presented here also show that the d(x 2 -y 2 ) eigenvalue and Γ get enhanced when U varies from 4t to 8t.

Model
The two-leg Hubbard ladder is defined by where t (t ┴ ) is the hopping parameter parallel perpendicular) to the chains. The operator c + iλσ (c iλσ ) creates (annihilates) an electron of spin σ at site i of chain λ, the electron occupation number is n iλσ =c + iλσ c iλσ , and μ is the chemical potential. In addition, periodic boundary conditions are used along the chains. In the ground state and for U=4t, the d(x 2 -y 2 ) pairing correlations are most enhanced near half-filling for t ┴ ≈1.6t [4]. When U=8t, this occurs for t ┴ ≈1.4t. Here, the QMC data will be presented using these parameter sets for a 2×16 lattice.
In obtaining the data presented here, the determinantal QMC technique described in Ref. [7] was used. The calculation of the reducible interaction Γ follows the procedure given in Refs. [8,9] for the twodimensional case. In this procedure, Γ is obtained from the QMC measurement of the two-particle Green's function Λ defined by where T τ is the τ-ordering operator, x i =(x i ,τ i ) denotes both space and Matsubara imaginary-time variables, and c σ (x i ) (c + σ (x i )) annihilates (creates) an electron of spin σ at x i . By Fourier transforming with respect to x i and τ i , one obtains where p=(p,iω n ), G σ (p) is the single-particle Green's function which is obtained separately by QMC, and Γ(p',k'|p,k) is the reducible particle-particle vertex. The single-particle Green's function is defined by and its Fourier transform by G(p,iω n )= ∫ 0 β dτ exp(iω n τ) G(p,τ) where ω n =(2n+1)πT is the fermion Matsubara frequency. The BCS component of the particle-particle vertex is Γ(p',-p'|p,-p), which will be denoted by Γ(p'|p) in the following. Figure 2 shows an illustration of Γ(p'|p), where the incoming fermions at (p,iω n ) with up spin and (-p,-iω n ) with down spin scatter to states (p',iω n' ) with up spin and (-p',-iω n' ) with down spin by exchanging momentum q=p'-p and Matsubara frequency ω m =ω n' -ω n . The reducible interaction Γ is related to the irreducible interaction Γ I by the t-matrix equation, Feynman diagram for the reducible particle-particle interaction Γ(p'|p) in the BCS channel.
In BCS, Γ I is the effective pairing interaction.
We perform a quantitative study of the pairing instability by solving the Bethe-Salpeter equation in the BCS channel which is given by with p=(p,iω n ) and ω n =(2n+1)πT. Here, λ α is the irreducible eigenvalue in channel α, and φ α (p) is the corresponding eigenvalue. For a three-dimensional infinite system, when the maximum λ α reaches 1, this signals a BCS instability to a state where the pair-wave function at T c has the form of the corresponding eigenfunction φ α (p,iω n ). In this paper, the discussion will concentrate on the momentum structure of the reducible interaction in the singlet channel, and on the temperature evolution of the leading eigenvalue of the Bethe-Salpeter equation.

QMC data
In this section, we will show QMC data for U=4t and t ┴ ≈1.6t, and U=8t and t ┴ ≈1.4t. In the following, the QMC data on Γ s (p',iω n' |p,iω n ) will be shown at the lowest Matsubara frequency ω n =ω n' =πT for these sets of model parameters. These QMC results will show that the structure in Γ s (p',iπT|p,iπT) depends on where p' and p are located in the Brillouin zone. In particular, Γ s depends strongly on where p and p' are located with respect to the Fermi surface or the saddle points. In order to display these features, Γ s (p',iπT|p,iπT) will be plotted as a function of p while p' is kept fixed at various points in the Brillouin zone. On the 2×16 lattice for t ┴ =1.6t and <n>=0.94, the tight-binding band structure ɛ p =-2tcos(p x )t ┴ cos(p y )-μ has Fermi surface crossing points near (±3π/4,0) and (±π/4,π). These momentum are represented by the filled circles in Fig. 3. The single-particle spectral weight A(p,ω) for the two-leg Hubbard ladder, which was calculated by the maximum-entropy analytic continuation of the QMC data find that due to the Coulomb correlations A(p,ω) gets redistributed and renormalized. However, for the 2×16 lattice with the parameters U=4t, t ┴ =1.6t, <n>=0.94 and T=0.1t, the 00006-p.2 Fermi surface crossing points are close to those of the noninteracting case, and they occur near (±3π/4,0) and (±π/4,π). In the Section 3.2, Γ s (p',iπT|p,iπT) versus p will be plotted for p' located at (0,π) and at (π/4,π).

Single-particle spectral weight
We begin presenting the QMC data by first showing the momentum and frequency dependence of the singleparticle spectral weight A(p,ω), because it plays an important role in determining the solution of the Bethe-Salpeter equation. We obtained A(p,ω) by using the maximum-entropy procedure from the QMC data on G(p,τ). Figure 4 shows A(p,ω) versus ω at different p=(p x ,p y ) for parameters T=0.1t, U=4t, t ┴ =1.6t and <n>=0.94. Here, the solid curves denote results for the bonding (p y =0) and the dotted curves are for the antibonding (p y =π) band. When interpreting the maximum-entropy images of A(p,ω), it is necessary to be cautious, since this technique has a finite resolution which worsens away from the Fermi level. In this respect, we note that the QMC data on G(p,τ) which were used in this procedure had an error covariance matrix which exhibited a continuous eigenvalue spectrum. In addition, the classical and the Bryan's maximum-entropy algorithms yielded identical results for A(p,ω). Here, we will use the maximum entropy results mainly for determining the Fermi wavevectors. In Fig. 4, for p y =0, we observe a peak which moves to higher energy as p x goes from 0 to 3π/4. At p=(3π/4,0), most of the weight in A(p,ω) is contained in the peak which is centered ≈0.3t below the Fermi level. As p goes from (3π/4,0) to (π,0), the main peak in A(p,ω) does not cross the Fermi surface but remains pinned right below it, while some of the spectral weight gets transferred to a second peak above the Fermi level. At p=(0,π) the spectral weight again consists of two peaks. The one which is lying lower in energy is located near the quasiparticle position of the noninteracting system. The second peak is about 1.8t higher in energy. As p moves from (0,π) to (π/4,π), the peak lower in energy has some of its spectral weight transferred to the second peak. In addition, at p=(π/4,π) it is seen that the lower peak develops a suppression at ω≈0. At higher temperatures, this suppression goes away and a single peak is recovered at the Fermi level. Since the maximum-entropy technique has finite resolution, this suppression of A(p,ω) at the Fermi level might be an artifact of the analytic continuation. However, it is also possible that A(p,ω) at the Fermi level is getting suppressed because of the strong scattering which develops for these wave vectors at low T, as we will next see in the QMC data on Γ. Finally, as p goes from (π/4,π) to (π,π), the lower peak loses more of its weight to the peak at higher energy, and they both disperse away from the Fermi level.

Pairing instability
In order to show the d(x 2 -y 2 ) symetry of the leading pairing instability, we present QMC data on the Bethe-Salpeter eigenfunction φ d (p,iω n ). In Figure 5, the momentum dependence of φ d (p,iω n ) is plotted for the lowest Matsubara frequency iπT. Here, we note that φ d (p,iπT) changes sign as p goes from ≈(π,0) to ≈(0,π) and, hence, φ d has d(x 2 -y 2 ) type symmetry. We also note that, for the antibonding band, the peak in φ d (p,iπT) occurs at p x =π/4. At higher temperatures, this peak moves to (0,π). On the other hand, for the bonding band we find that φ d (p,iπT) has maximum amplitude at (π,0) at these temperatures.

QMC data on Γ s for U=4t
Next, we discuss the QMC results on the momentum dependence of reducible particle-particle interaction in the singlet channel. Figure 6 shows Γ s (p',iπT|p,iπT) for <n>=0.94 at various temperatures. In Fig. 6(a), Γ s (p',iπT|p,iπT) is plotted as a function of p while p' is kept fixed at (π/4,π). In the left panel of Fig. 6(a), p x is scanned from -π to π while p y is kept fixed at 0, and in the 00006-p.3 right panel p x is scanned while p y =π. Here, we observe that sharp structures develop in Γ s as T decreases from 0.25t to 0.1t. In the left panel of Fig. 6(a), we see that a peak develops in Γ s when p≈ (-3π/4,0). The location of this peak corresponds to a scattering process involving momentum transfer q=(π,π) (backward scattering). In addition, in the right panel of Fig. 6(a), we see that a dip develops in Γ s when p is located near p'=(π/4,π), corresponding to zero momentum transfer (forward scattering). This dip is because of resonant scattering in the d(x 2 -y 2 )-wave BCS channel, which will be further discussed below. In Section 4, we will see that when a d(x 2 -y 2 )-wave superconducting instability is approached in a three-dimensional infinite system, Γ s (p',iπT|p,iπT) on the Fermi surface diverges to +∞ for scatterings involving large momentum transfer, and to -∞ for q=0 scatterings. In Fig. 6(a), we observe that Γ s is developing large repulsive and attractive peaks for scatterings on the Fermi surface for T≤0.125t. The development of the dip in Fig.  6(a) means that resonant scattering in the d-wave BCS channel is taking place already at these temperatures.
In Fig. 6(b), Γ s (p',iπT|p,iπT) versus p is plotted in the same way as in Fig. 6(a) but for p' located at the saddle point (0,π). In this case, Γ s (p',iπT|p,iπT) develops a peak for p≈(±π,0), corresponding to scatterings with q=(π,π) momentum transfer, and the magnitude of this peak is comparable to that seen in the left panel of Fig.  6(a). However, the behavior for q=0 momentum transfer is different. In the right panel of Fig. 6(b), we observe that Γ s (p',iπT|p,iπT) for q=0 scattering remains pinned at zero for T down to 0.125t, and becomes attractive only at T=0.1t. Hence, in this case, the resonant scattering in Γ s for q=0 momentum transfer is significantly weaker. Figure 7 shows the temperature dependence of the backward and the forward scattering components of Γ s . Here, Γ s is plotted as a function of the inverse temperature β for momentum transfers q=(π,π) (filled circles) and (π/8,0) (open circles), while p' is kept fixed at (π/4,π). In this figure, we observe that the backward scattering increases almost linearly with β in this temperature range. We also note that Γ s for q=(π/8,0) becomes negative for β>4/t and exhibits rapid growth as Fig. 6. Reducible particle-particle interaction in the singlet channel Γ s (p',iπT| p,iπT) versus p x for p y =0 (left panel) and p y =π (right panel). In (a), p' is kept fixed at (π/4,π) and in (b), p' is fixed at (0,π). These results are for U=4t, t ┴ =1.6t and <n>=0.94. Fig. 7. Reducible particle-particle interaction in the singlet channel Γ s (p',iπT| p,iπT) versus the inverse temperature β for backward (q=(π,π)) and forward (q=(π/8,0)) scatterings. Here, p=p'+q and p' is kept fixed at (π/4,π). These results are for U=4t, t ┴ =1.6t and <n>=0.94. β increases further. These results on Γ s show that d(x 2 -y 2 ) pairing correlations develop rapidly at low temperatures in this system.

QMC data on Γ s at half-filling
In Figure 8, we present results on the momentum dependence of Γ s (p',iπT|p,iπT) at half-filling. Here, we observe that these data are qualitatively different than those for <n>=0.94 seen in Fig. 6. For example, as T decreases, we see that the q=(π,π) scatterings grow faster. On the other hand, the q≈0 resonant scatterings in the dwave channel are significantly weaker in comparison to the q=(π,π) scatterings. Hence resonant scattering in d(x 2y 2 )-wave channel does not take place at half-filling. The growth of Γ s for q=(π,π) scatterings follows the development of the AF correlations, which grow rapidly in this T range. However, the AF correlations saturate at lower T, since the ground state of the half-filled Hubbard ladder is an insulator with short-range AF correlations.

QMC data on Γ s for U=8t
Next, we investigate how Γ s and λ d changes as U increases from 4t to 8t. We will see that the overall features of the momentum dependence of Γ s (p',iπT|p,iπT) are similar for U=4t and 8t, however, the magnitudes of Γ s and λ d are larger for U=8t. Figure 9 shows the temperature evolution of Γ s (p',iπT|p,iπT) for U=8t and <n>=0.94. In this figure, we see that, at T=0.25t and for q≈(π,π), Γ s becomes as large as 100t, which is an order of Fig. 9. Reducible particle-particle interaction in the singlet channel Γ s (p',iπT| p,iπT) versus p x for p y =0 (left panel) and p y =π (right panel). Here, Γ s is shown for p'=(π/4,π), U=8t, t ┴ =1.4t and <n>=0.94. magnitude larger than the bare Coulomb repulsion. Here, we also observe that the q≈0 scatterings become attractive with decreasing T, when p' is located at (π/4,π). However, for p'=(0,π), the q≈0 scatterings are weaker (not shown here). Hence, the momentum structure in Γ s for U=8t is similar to that for U=4t, except for the magnitude.

QMC data on the Bethe-Salpeter eigenvalue λ d
The T dependence of λ d is shown in Figure 10(a) for U=4t and t ┴ =1.6t. Here, it is seen that, as T decreases, λ d grows monotonically eaching 0.75 at T=0.1t for <n>=0.94. Upon doping to <n>=0.875, λ d decreases. Also shown in Fig. 10(b) is λ d for U=8t and t ┴ =1.4t. At the temperatures where these calculations were performed, we find that λ d is larger for U=8t.

Fermion sign problem
The lowest temperatures at which these QMC calculations can be performed are restricted by the fermion sign problem [10]. In particular, the statistical error in the QMC data grows exponentially as the average sign of the fermion determinants <sign> decreases from unity. In the determinantal QMC method, the operator expectation values are evaluated at finite temperatures using the grand canonical ensemble. For the Hubbard model within this approach, first the Trotter breakup is applied.
Afterwards the Hubbard-Stratonovich transformation is performed introducing an Ising spin variable S iℓ at each lattice site i and Matsubara time slice τ ℓ . In the following step, the fermion degrees of freedom are integrated out. The resulting expectation values are evaluated by averaging over the configurations of the spin variables {S iℓ } with the Monte Carlo technique. For instance, the partition function Z=Tr e -βH is obtained, within this approach, from (11) ).
Hence, when Γ(k|p) is sharply peaked for k≈p+(π,π), Γ(p|p) becomes attractive. This is the origin of the strong attractive scattering in Γ s (p',iπT|p,iπT) found for q≈0 momentum transfers at the Fermi surface. In order to demonstrate how the repeated scatterings in the BCS channel affect the particle-particle interaction, it is also useful to consider the case of an irreducible interaction which is independent of frequency and separable in momentum. In this simple case, the irreducible interaction has the form where α denotes the various pairing channels, V α is the component of Γ I in channel α and g α (p) is the corresponding form factor. Through repeated scatterings in the BCS channel, the reducible becomes where the bare pair-field susceptibility  P in channel α is defined by In general, the irreducible interaction has both attractive and repulsive components. Here, it is clear that for repulsive components the momentum structure in Γ I gets suppressed by particle-particle scatterings in the t-matrix. However, if V α is attractive, then this component gets enhanced. If the d(x 2 -y 2 )-wave component is attractive and gets sufficiently enhanced, then the reducible Γ(p'|p) in Eq. (13) becomes attractive for q=p'-p≈0. The resonant scattering for q≈0 momentum transfer on the Fermi surface, which was observed in the QMC data on Γ s , is due to this process. From Eq. (13), it is also seen that, if a d(x 2 -y 2 )-wave BCS instability is approached, Γ(p'|p) diverges to +∞ for backward scattering and to -∞ for forward scattering. It is useful to note that, in these QMC data, strong resonant scattering is observed on the Fermi surface at T=0.1t for U=4t and 6% doping, whereas for U=8t the resonant scattering is taking place already at T=0.25t.

Summary and conclusions
In this paper, the magnitude and the characteristic temperature scale of the d(x 2 -y 2 ) pairing correlations have been investigated for the two-leg Hubbard ladder. For this purpose, QMC calculations have been carried out for the reducible particle-particle interaction in the singlet BCS channel, Γ s, and the corresponding Bethe-Salpeter equation has been solved. In particular, the momentum dependence of Γ s has been studied as a function of the model parameters. An important result of these QMC calculations is the observation of resonant scattering in the d(x 2 -y 2 ) BCS channel. At the Fermi surface, the particle-particle scatterings with q≈(π,π) momentum transfer become strongly repulsive, while scatterings with q≈0 become strongly attractive as the temperature is lowered. It was seen that the backward and the forward scattering amplitudes can attain values which are large compared to the bare bandwidth or the bare Coulomb repulsion depending on the model parameters. Another finding of these calculations is that the forward resonant scattering is anomalously weak at the saddle points. The momentum structure of Γ s is also reflected in the momentum dependence of the d(x 2 -y 2 )-wave eigenfunction of the Bethe-Salpeter equation. These features are a signature of the underlying pairing interaction in the two-leg Hubbard ladder.
In order to discuss quantitatively the strength of the d(x 2 -y 2 ) pairing correlations, results were shown for the d(x 2 -y 2 ) eigenvalue of the Bethe-Salpeter equation in the BCS channel. For a three-dimensional system, the BCS instability occurs when the leading Bethe-Salpeter eigenvalue reaches one. Here, it was found that, for U=4t, t ┴ =1.6t and <n>=0.94, λ d increases monotonically as T is lowered and becomes as big as 0.75 at T=0.1t. In these 00006-p.7 calculations it was also seen that λ d increases as U varies from 4t to 8t. The QMC data presented in this paper provide useful information about the magnitude and the temperature scale of the pairing correlations in the twoleg Hubbard ladder.
In summary, we have investigated the d(x 2 -y 2 ) pairing correlations in the Hubbard ladder by presenting QMC data on the T, U and the doping dependence of the d(x 2 -y 2 )-wave eigenvalue of the Bethe-Salpeter equation. These data show that there are strong d(x 2 -y 2 ) pairing correlations in the Hubbard ladder for certain values of the model parameters.