Pairing Correlations in the two-layer attractive Hubbard Model

Studies of systems with two fermionic bands with repulsive interaction strength U have a long history, with the Periodic Anderson Model (PAM) being one of the most frequently considered Hamiltonians. In this paper, we use Quantum Monte Carlo to study analogous issues for attractive interactions. As in the Periodic Anderson Model, we focus on a case where one band is uncorrelated (U=0), and focus on the effect of hybridization V between the bands on the pairing correlations. A key difference with the PAM is that there is no sign problem, so that we are able to explore the physics of doped multi-band attractive systems at low temperatures whereas ground state properties of repulsive models can be determined only at half-filling. For small V, pairing in the U<0 layer induces pairing in the U=0 layer. At larger V the ground state of the coupled system loses its superconducting character. The Quantum Monte Carlo data are complemented by results obtained with the Bogoliubov-de Gennes approximation.

that the physics of fermions in multiple layers can equivalently be interpreted in terms of multiple orbitals.
The purpose of this paper is to study the possibility of analogous phenomena in the case when there is an attractive interaction. More specifically, we examine the nature of pairing correlations for coupled planes (orbitals) in situations when an attraction is present only in one plane (orbital). We use the U < 0 Hubbard Hamiltonian as an appropriate simple model, and employ a combination of QMC and Bogoliubov-de Gennes (BdG) mean field theory. We will focus on a case analogous to that of the PAM, namely when U = 0 in one of the layers. We are interested in both how the interlayer hopping affects the s-wave pairing correlations in a layer with U < 0, and also whether pairing can be induced by a 'proximity effect' in the U = 0 layer. Such behavior would be directly analogous to the 'Ruderman-Kittel-Kasuya-Yoshida (RKKY)' spin polarization cloud which is induced in the conduction band by the presence of the local moments in the PAM. A particularly interesting issue is whether the QPT which occurs at commensurate filling is present also in the doped case.
Bilayer Hubbard Hamiltonians with attractive interactions which differ in the two layers have also been studied by Berg et al [8] and Wachtel et al [9]. The motivation is to explore the possibility of enhanced superconductivity, with transition temperatures approaching the mean field pairing scale 0 , by separating a spatial region where there is strong pairing but low phase stiffness, from a metallic region. The emphasis of this earlier work was on the case when the bandwidth W is much larger than 0 , whereas the focus here will be on an 'intermediate coupling' regime where W ∼ 0 . Our numerical work complements these previous methods, since we cannot easily access that part of parameter space which is in the Bardeen-Cooper-Schrieffer (BCS) limit and hence has exponentially small transition temperatures, regions which are however accessible to analytic studies. Experimental studies of superconductivity in systems in which one layer (BaBiO 3 ) provides the electron pairing interaction and the other layer (BaPbO 3 ) is conducting, which test some of the ideas of [8,9] have also been reported. Indeed, the role of inhomogeneity, arising either spontaneously or introduced explicitly into the Hamiltonian, has been studied in a variety of contexts within the repulsive Hubbard model as well [10,11].
The remainder of this paper is organized as follows: in section 2 we write down the precise model and give an overview of the QMC and BdG computational methodologies. Sections 3 and 4 present the BdG and QMC results, respectively. Section 5 recaps our conclusions.

Model and computational methods
Our starting point is the two layer (orbital) attractive Hubbard Hamiltonian Here an intralayer kinetic energy term c † ilσ c jlσ describes the destruction of a fermion with spin σ on site j of layer l and its creation on site i of the same layer. A second interlayer kinetic energy term c † ilσ c il σ hybridizes the two layers l, l . Fermions of different spin feel an attractive interaction −|U l | on layer l. The chemical potentials µ l control the filling. Our lattice geometry consists of two coupled square lattice of linear size L. We choose U 1 < 0 in the 'SC' layer and U 2 = 0 in the metallic one. At half-filling, a particle-hole transformation formalizes the similarity of the repulsive and attractive models. The vanishing of AF order with increasing V for U > 0 bilayers immediately implies that ground state pairing (and charge density wave order, with which it is degenerate) must be destroyed as V grows for U < 0. However, whereas AF order occurs only at half-filling (and only at T = 0) in the 2D repulsive Hubbard Hamiltonian, superconductivity appears away from half-filling, and at a finite (Kosterlitz-Thouless) transition temperature [30,31], in the attractive Hubbard Hamiltonian. This opens up fundamentally new issues in the attractive case. The manner is which fermions on sites in adjacent layers lock into local objects and the interplay with LRO when the filling is incommensurate is an interesting question. The absence of a sign problem in the case of the attractive Hubbard model enables the study of low temperatures even when the filling is incommensurate, allowing access to ground state properties.
We perform 'determinant' QMC [14] simulations of the Hamiltonian, equation (1), by writing a path integral for the partition function Z . The exponential of the interaction term U l is decoupled with a Hubbard-Stratonovich field which allows the trace of the remaining exponentials of quadratic forms of the fermion operators to be performed analytically. The result is a sum over configurations of the Hubbard-Stratonovich field with a weight which takes the form of the product of the determinants of two matrices of dimension the spatial lattice size, one for each spin species. In the case of attractive interaction the two matrices, and hence their determinants, are identical, so that the weight is a perfect square and there is no sign problem. In our simulations we have used a discretization β = L τ of the inverse temperature β with τ = 0.125, and an exact exponentiation of the kinetic energy operator. These choices ensure that the systematic 'Trotter' errors are smaller than the statistical error bars for the pair structure factor. Most of our runs represent the averaged results of 40 independent simulations of 1000 measurement sweeps, specifically ten runs for each type of boundary condition (periodic/antiperiodic) in the x and y directions. This boundary condition averaging can be thought of as sampling additional discrete momenta values intermediate to those with just periodic boundary conditions, and is known to reduce finite size effects somewhat.
We present results here for the real space pair correlation functions in the two orbitals and also their associated structure factors as functions of interband hybridization V , on-site attraction U l , density n and temperature T . The long distance behavior of P l ( j) yields the square of the SC order parameter, as does the extrapolation of the structure factor to the thermodynamic limit.
Other quantities of interest are the average near-neighbor intra-and inter-layer hopping and the near-neighbor intra-and inter-layer pairing correlations whereδ ≡ (1, 0) or (0, 1). In order to study the possibility of superfluidity mediated by interlayer pairs, we also define the pairing correlations is the inter-layer pair creation operator [15]. Equation (1) can also be studied in mean field theory via the solution of the BdG equations. In this approach, the four fermion (interaction) term is decoupled in the Hamiltonian itself, Here the gap il = c il↑ c il↓ and density n ilσ = c † ilσ c ilσ are determined self-consistently by diagonalizing the quadratic BdG Hamiltonian and putting together the eigenvalues and eigenvectors appropriately to compute refined values which are inserted back in the Hamiltonian in an iterative process. A related BdG treatment of superconductivity in the presence of randomness (spatially varying µ i ) is given in [16,17].μ il = µ + |U l | n il /2 includes a sitedependent Hartree shift with n il = σ n ilσ .
In the general case, when inhomogeneous terms are present in the Hamiltonian, or are expected to develop spontaneously, (like the charge and spin stripes of the doped, repulsive Hubbard model), the order parameter and densities are allowed to depend on site index and the diagonalization must be done numerically. If translation invariance is present, equation (8) can be diagonalized analytically by going to momentum space. The resulting momentum sums are typically still done numerically, but larger lattices can be studied. We have implemented both approaches in the work presented here.
Although inhomogeneous Hartree-Fock theory allows for spatial variation of the densities and SC order parameters, the expectation values in equation (8) are independent of imaginary time, unlike the fluctuating Hubbard-Stratonovich field in the QMC approach. This leads to a less accurate treatment of interparticle correlations, but a numerically much more simple problem. In particular, larger spatial lattices can be studied, and the BdG solution can exhibit broken symmetry so that l itself can be non-zero, rather than having to be extracted from the asymptotes of the correlation function equation (2). Hence the combination of QMC and BdG serves to provide complementary information.
We present only a brief set of results of BdG calculations, focusing on the values of the order parameter l as a function of the same parameters as varied in the QMC. As mentioned above, the asymptotic value of P l s ( j) is a measure of 2 l . The BdG approach, since it neglects fluctuations, yields larger values of l than those of the QMC, which is an exact method. The main conclusion, is that, despite this difference in size of the order parameter, the qualitative physics-a dome-like induced pairing in the non-interacting layer-is similar in BdG and QMC.

Bogoliubov-de Gennes results
The results of BdG calculations are summarized in figure 1. We consider two layers with U 1 /t = −8 and U 2 = 0, at a density of n = 0.8 electrons per site and study the evolution of the pairing correlation, l , in the two layers as a function interlayer coupling, V . Pairing in the correlated layer is suppressed by coupling to the non-correlated layer and eventually destroyed completely at V ∼ 3.6 t for β t = 8. Whereas for U 1 > 0 there is a genuine quantum critical point so that the AF is destroyed above V c even at T = 0, the destruction of superconductivity occurs here in the attractive case only because T is finite so that, as V increases, eventually the SC transition temperature T c falls below T = t/8. Finite size effects are small at this temperature, so that a 12 × 12 lattice already captures the thermodynamic limit. For lower T , especially in regimes where the SC transition temperature is expected to be exponentially small, it is necessary to study larger lattices where discrete level spacings do not affect the results. A further interesting feature of figure 1 is that the l = 2 layer, with U 2 = 0, develops induced pairing with the onset of interlayer coupling. The pairing amplitude, P 2 , increases with increasing V up to a maximum before decreasing and eventually vanishing. Crucially, the pairing order parameter in both layers remains non-zero over a finite range of interlayer coupling-the induced pairing reported previously at half-filling extends to finite dopings. In the next section, we present extensive QMC results to confirm and complement these basic BdG results.
The dependence of the correlation functions, P 1,2 ( j), on separation j in a bilayer system comprised of one correlated layer coupled to an uncorrelated layer for two different lattice sizes L. The system is at half-filling, n = 1, with interlayer hybridization V = 0.8t, and on-site interaction strength U 1 = −6 t for the correlated layer. The inverse temperature is β t = 12 and is sufficiently low that the ground state has been reached. The correlation functions converge to non-zero value at large separations, providing clear evidence for LRO. Despite the absence of interactions, there is proximity-effect induced LRO in the uncorrelated layer, l = 2, with an order parameter 2 approximately a factor of three lower than for the correlated layer. As in figure 1, finite size effects are modest.

Half filling
We begin our discussion of QMC results at half-filling, n 1 = n 2 = 1. As discussed, this corresponds to the particle-hole symmetric point µ 1 = µ 2 = 0 where there exists an exact mapping to the repulsive Hubbard model and the results can be benchmarked against previous studies for U > 0 [7,12]. Lattices of the form 2 × L × L, with 4 L 14 were simulated with periodic boundary conditions over a wide range of interlayer hybridization, and on-site interaction for the interacting layer. An inverse temperature βt = 12 was used. In many of the cases we studied this was found to be sufficient for the observables to have converged to their ground state values. However, in some parameter regimes, the transition temperatures assume a BCS form and vanish exponentially, so that the pairing order is not easily accessible. The spatial dependence of the pair correlations P l ( j) at half-filling (n = 1), interlayer hybridization V = 0.75t, and interactions U 1 = −6 t, U 2 = 0 t is shown in figure 2. The separation j follows a trajectory along the x-axis to maximal x separation ( L 2 , 0) on a lattice with periodic boundary conditions, and then to ( L 2 , L 2 ) before returning to separation (0, 0). Results for lattices with L = 8, 12 are shown in the figure. The interacting layer l = 1 exhibits clear LRO with P 1 ( j) nearly independent of j beyond a separation of approximately L 4 . The non-interacting layer l = 2 also exhibits proximity-effect induced LRO, although the correlation function is roughly three times smaller. This corresponds to an order parameter 2 ∼ 1 /3.
where the linear in 1/L term is the expected spin-wave correction to the thermodynamic limit value P l s (∞) in the ordered phase [32], and the quadratic in 1/L 2 term describes the 1/N fall-off of the structure factor which occurs when the real-space correlations are short-ranged (see equation (3)). We found the structure factor data to fluctuate more strongly at small V , precluding the acquisition of data for linear dimension L > √ 72. For larger V , lattices of up to L = 12 were studied.
The extrapolated results are then given in figure 4. Increasing the interlayer hybridization reveals that the pair structure factor for the interacting layer, P 1 s , decreases monotonically and eventually the intra-layer pairing is destroyed at some critical V c . On the other hand, the pairing structure factor for the non-interacting layer P 2 s varies non-monotonically with V . It vanishes (or is exponentially small at low but finite T ) as V → 0. As V increases, long-range pairing correlations become more robust, reach a maximum at an intermediate V (which, as we shall later show, depends on the strength of on-site interaction in the interacting layer) and then decrease continuously to zero at V c . For V > V c , the ground state is dominated by inter-layer singlets. Obtaining a precise value of V c is complicated by the error bars in the raw data for each lattice size and then additional uncertainties in the extrapolation of figure 3. Our data suggest that P 2 s and P 1 s are zero to within error bars for V 2.0 and have certainly vanished for V > 2.5.

Incommensurate filling
We now turn to the case of doped layers. As mentioned in the introduction, this is interesting for two reasons. Firstly, there is a possibility of a finite temperature Kosterlitz-Thouless transition because the charge density wave (CDW)-SC symmetry is broken. Secondly, simulations of coupled, doped repulsive layers at low T are not possible. Hence there is no determinant QMC (DQMC) information available for the nature of magnetism in coupled layers away from halffilling. Our simulations address this issue, albeit for attractive on-site interactions. The results in figures 5 and 6 confirm that induced pairing in the uncorrelated layer extends to finite dopings away from half-filling; indeed the pair correlation functions in the two layers are qualitatively similar to those at half-filling. Pairing in the correlated layer decreases monotonically to zero at a finite V c . In the non-interacting layer, pairing increases rapidly from zero as the interlayer hybridization is turned on, reaches a maximum value at intermediate V and then decreases to zero. We expect however, that unlike the half-filled case where the vanishing of pairing at large V occurs even at T = 0 (reflecting a QPT), in the doped case it is most likely that pair correlations become short ranged because the transition temperature falls below the simulation temperature as V increases. Because the transition temperature is suppressed exponentially in the BCS regime, it is not possible to perform the simulations at the extremely low T needed to settle this issue.
At intermediate to strong values of the on-site interaction (|U | 6 t), the ground state of a system of coupled correlated and uncorrelated layers away from half-filling consists of intra-layer pair formation in both layers over a finite non-zero range of interlayer coupling. Eventually, at sufficiently strong interlayer hybridization, the system is characterized by singlet formation between the layers, and pairing occurs only at a (much) reduced temperature. We find here that induced pairing is found to be absent for |U | 6t at an inverse temperature βt = 12. The simulation results are summarized in figure 6 where the extrapolated pairing structure factors for the two layers are shown as a function of the interlayer hybridization V for some representative values of the on-site interaction in the correlated layer, and at a fixed value of the density, n = 0.8. The maximum of the pairing structure factor in the uncorrelated layer and the critical value of interlayer coupling for suppression of pairing increases with |U |. Our results are consistent with the disappearance of pairing in both layers occurring simultaneously in these finite T simulations. The error associated with determining the extrapolated pairing structure factor in the thermodynamic limit constrains the accuracy with which V c can be determined to about ±0.2t. In a marked departure from the BdG results, the QMC data show that a finite nonzero interlayer coupling is required for the onset of induced pairing in the uncorrelated layer at βt = 12.  Figure 6 shows that P 1 s increases systematically for V > 0.5 as U changes from U = −4 to −10. In the case of the repulsive Hubbard Hamiltonian at half-filling the magnetic response first increases with U but then reaches a maximum at U ∼ 8t before falling. This behavior is understood qualitatively from the fact that the superexchange J = 4 t 2 /U which provides the large U magnetic energy scale declines with U . Thus, for example the Néel temperature of the 3D Hubbard model at half-filling is maximized at U ∼ 8t. The analogue in the attractive Hubbard model is that the pairs become heavy, with an effective hopping t eff ∼ t 2 /U as a consequence of the (temporary) pair-breaking that must occur as a pair hops. This analogy suggests the pairing response might also be maximal at an intermediate U . Because of the sign problem, it is not known from DQMC whether such a non-monotonic behavior exists in the repulsive model away from half-filling, or, if it does, at what value of U the maximum occurs. In any case, the issue of an optimal U is a bit subtle. At a fixed non-zero temperature, the t 2 /U energy scale can fall below T , leading to a decrease in order at strong coupling. However, strictly at T = 0 the AF order parameter in the repulsive Hubbard model, for example, increases monotonically to the large U Heisenberg value [13].
In some cases the data are more well behaved (smaller error bars) at ρ = 0.8 (figures 5 and 6) than at ρ = 1.0 (figures 3 and 4). We believe this is a consequence of the degeneracy between charge density wave and SC order at half-filling. There is a tendency, for ρ = 1 and especially at low T , for the system to get stuck in one of the competing ordered phases, leading to larger fluctuations between runs. This degeneracy is lifted away from half-filling. Throughout these calculations an important question is whether the temperatures are low enough for the ground state, and its associated order, if any, to have been reached. Figure 7 provides one measure of this by showing the pair structure factors versus β for n = 0.8, U = −10 and V = 1.2. From such plots we conclude that β t ≈ 12 is sufficiently cold. Hence this is the value of inverse temperature chosen throughout this paper. It is important to emphasize, however, that the required value of β will depend on filling, lattice size and the energy scales of the model. Small lattices generally reach their ground states at higher temperature, since it is easier for the correlation length to exceed the lattice size. More problematic than large lattices, however, are situations where an exponentially small BCS pairing temperature scale is possible. We expect this to occur for doped lattices at small V . Hence it is possible that some situations where the extrapolated structure factor vanishes arise from the pairing scale falling below the temperature at which the simulation is performed. An exponentially small scale is not possible to pursue numerically.
Having established the existence of induced pair formation at finite doping and the evolution of the associated structure factor with V , we investigate the dependence of the pairing phenomenon on the density of electrons in the layers. Our results (figure 8) show that induced pairing occurs generically at all densities and that the pairing structure factor behaves in a qualitatively similar manner. In the single layer (2D) attractive Hubbard Hamiltonian, T c rises abruptly from zero at half-filling where charge density order competes with pairing, but then shows a gradual decline below an optimal doping [30]. Our data in figure 8 suggest the increase in pairing might be saturating in the range ρ ∼ 0.6-0.8. If this two layer model were to behave like the single layer case, this saturation would be followed by a decline with further doping away from ρ = 1.

Nature of the phase at V > V c
In section 4 we saw that when the interlayer hybridization, V , is larger than a threshold value V c , the intra-and inter-layer pair structure factors are strongly suppressed at fixed, non-zero T . While the value of V c depends on the filling and the interaction strength in the correlated layer, we found it to be typically V c ∼ 2t. The question then arises as to the nature of the correlations for V > V c .
To address this question, we show in figure 9 the average intra-and inter-layer hopping, equation (4), and the short range (near neighbor) pair correlations, equation (5) for the half filled system. We see that as V increases, the interlayer near neighbor hopping, T 12 , and pairing correlations, P 12 , increase monotonically and saturate when V > V s with V s ≈ 4 t. At the same time, the corresponding intra-layer quantities decrease and become negligible V s ≈ 4 t. This value of V s is quite insensitive to the value of U 1 .
This behavior can be understood by recalling that in the limit of U 1 = 0, there is a sharp phase transition from a metal to a Mott insulator when a finite gap opens in the dispersion relation, E ± (k x , k y ) = ±V − 2t ( cos k x + cos k y ) for V = V s = 4 t. One way of interpreting the results of figure 9 is that, as |U 1 | is increased from zero, a crossover line continues above this non-interacting critical point at U 1 = 0. Our results indicate that the trajectory of this crossover line is fairly vertical in the V -U 1 layer, i.e. V s is only weakly U 1 dependent.
Our results indicate that, for half-filling, in the region above, but close to, V c ≈ 2t, long range pairing order vanishes, but short range intralayer correlations are still present. These decrease as V grows. When V > V s ≈ 4 t, a crossover occurs in which short range intralayer correlations are completely suppressed-the particles in interlayer nearest neighbor sites form independent pairs. Inset: on-site and near neighbor interlayer density-density correlations. The onset of interlayer pairing corresponds to the vanishing of interlayer nearest neighbor up-up density-density correlation at V ≈ 3t. Error bars for these local quantities are of the order of, or smaller than, the symbol size.
When the system is doped, it exhibits some similarities to the half-filled case. The results for T ll and P ll are shown in figure 10 for n = 0.8. Densities n = 0.6, 0.4, 0.25 and 0.2 (not shown), all display the same qualitative behavior. The main difference between half-filled and doped lattices is that when T 12 reached its saturation value for V > V s , T 11 and T 22 do not vanish. The reason is that the interlayer pairs which form due to strong interlayer hybridization are still mobile, the band is not full due to the doping. To study the possibility that these mobile interlayer pairs might condense to a superfluid, we measured the interlayer pair correlations equation (6) for several fillings, n, and values of V . We find no evidence for such superfluid behavior. This is likely to be caused by the finite T at which the QMC is performed. If the pairs are in the BCS regime, the temperature at which SF will develop will be exponentially low.

Conclusions
In this paper we have used QMC to elucidate the properties of a multi-orbital attractive Hubbard Hamiltonian. Our key conclusions are the quantification of the induced pairing in a noninteracting orbital by superconductivity in a correlated orbital, and, conversely, the suppression of pairing in an orbital with an attractive interaction through its coupling to a non-interacting orbital, in the intermediate coupling regime U 1 ∼ W = 8t. The general trends of the evolution of the order parameters are similar to those found in BeG mean field theory.
As we have emphasized in the introduction, analogous issues for repulsive models have a long history. In fact, over the last 4-5 years, the question of 'orbitally selective Mott transitions', has been extensively explored. The fundamental objective has been an examination of the density of states to determine whether or not one subsystem can be metallic (zero energy gap at E F ) and the other insulating (non-zero gap), or whether the coupling forces the transition to occur simultaneously in all orbitals (layers) [18][19][20][21][22][23][24][25][26][27][28][29]. Our work suggests that for attractive interactions the suppression of superconductivity at large inter-orbital hybridization V seems typically to occur simultaneously. It is to be noted, however, that it appears that at weak V pairing can exist in the correlated orbital before it is induced in the non-interacting one. That is, we find that for the attractive case at small interlayer coupling, SC exists only in the interacting layer despite the fact that at larger V the vanishing of SC appears at a common V c . This paper has focused on the evolution of the pair structure factor with the hybridization V between the attractive and non-interacting layers at low temperatures. At weak V there is no LRO in the U = 0 layer, however there is an onset of induced pairing at intermediate V ∼ 1.5t if U is sufficiently large. Further increase in V returns the layer to a phase in which long range pairing is suppressed. This non-monotonic behavior is qualitatively similar to the evolution of T c reported in [8,9], despite the rather different parameter regimes being studied there (|U | ∼ t and no hopping, or only one dimensional hopping, between the attractive U sites). This induced SC dome found in QMC is similar to that seen in the BdG treatment. The doping dependence of the critical temperature in a system with localized attractive orbitals coupled to a conduction band has also been explored [34].
Our work is related to several other recent determinant QMC and Lanczos studies of the attractive Hubbard model. Assman et al [33] found that in an attractive Hubbard model with a smoothly varying chemical potential (such as would be present in a confined ultracold atomic cloud) pairing is significantly increased in the half-filled portion of the lattice through its contact with doped regions. That is, the suppression of superconductivity caused by the appearance of a degenerate charge density wave phase at n = 1 is eliminated. Paiva et al [35] explored a one-dimensional model of borocarbides in which U < 0 and U = 0 sites alternate. Superconductivity is found to be possible only above a critical density.
We note that it was recently shown in repulsive models with more than two layers that there can be further non-trivial evolution of magnetic properties with the hopping across an interface between a Mott insulator and a metal [36]. It would be interesting to pursue related questions in the attractive model and, in particular whether the formation of a layer of local pairs at the boundary can lead to a shielding of penetration effects.
Finally, we investigated the nature of a state emerging at high interlayer hybridization V and find that it creates interlayer pairs, which at half-filling are stationary in the sense of a vanishing (or very small) intralayer kinetic energy T ll . At incommensurate filling these pairs remain mobile at large V , with finite T ll . We do not find evidence for the presence of any associated long-range pairing order, but that is most likely due to the finite T at which the QMC simulations are performed. One expects that for doped bands in two dimensions, pairing will cause superfluidity at T = 0.