Many-variable variational Monte-Carlo study of superconductivity in two-band Hubbard models with an incipient band

We study superconductivity in two-band models where one of the bands does or does not intersect the Fermi level depending on the parameter values. Applying a many-variable variational Monte-Carlo method for a Hubbard model on two-leg ladder and bilayer square lattices, we show that superconductivity can be enhanced in a parameter regime where the edge of one of the bands is near the Fermi energy, that is, when the band is incipient. The resemblence of the present results to those obtained by a weak coupling method in a recent study suggests that, even in the large $U$ regime, the suppression of the near-zero-energy spin fluctuations and the development of finite energy spin fluctuations are the key factors for the enhancement of superconductivity by an incipient band.


I. INTRODUCTION
Purely electronic mechanism of superconductivity is expected to exhibit extremely high T c due to the large energy scale of the pairing glue originating from quantum fluctuations, such as spin and charge, orbital fluctuations. Especially, spin-fluctuation-mediated pairing is one of the leading candidate mechanisms at work in unconventional high temperature superconductors, namely, cuprates and iron-based superconductors.
In particular, in the early days of the study of the iron-based superconductors, it was considered that the Fermi surface nesting between electron and hole Fermi surfaces, combined with Hubbard U , induces spin fluctuations, which in turn act as pairing interaction around a certain wave vector Q if the gap sign changes across Q. This kind of superconducting gap is referred to as the s ± pairing [1][2][3][4][5][6]. However, the spin-fluctuation theory has been challenged by the discovery of (heavily) electrondoped iron-based superconductors with a relatively high T c where hole-like bands sink below the Fermi level leaving only electron-like Fermi surfaces [7][8][9][10][11][12][13][14][15][16][17][18]. Naively, removing the hole pocket is expected to destroy the spinfluctuation-mediated pairing interaction and suppress T c rapidly.
After these observations, "incipient bands", which sit close to, but do not intersect the Fermi level, have received much attention. Various authors have suggested that the spin-fluctuation scattering of pairs between an electron Fermi surface and an incipient hole band can induce s ± pairing [1,15,17,[19][20][21][22][23][24]. In this context, the bilayer Hubbard model, which has been extensively studied in the past [25][26][27][28][29][30][31][32][33][34] has recently attracted renewed focus. Having hole and electron Fermi surfaces, it can be regarded as a single-orbital analogue of the iron-based superconductors. In fact, it has been found in previous studies that s ± pairing is favored over d x 2 −y 2 pairing by * kuroki@phys.sci.osaka-u.ac.jp increasing the relative strength of the inter-layer nearest hopping to the intra-layer nearest hopping [35]. Further, as one of the bands becomes shallow or incipient, the spectral weight of spin-fluctuation is transferred to higher energies, which can lead to s ± pairing state in which a gap appears on the hole band with the opposite sign to the gap on the electron Fermi surfaces [36][37][38].
Regarding the incipient band situation, it was proposed in ref. [39] by one of the present authors and his coworkers that strongly enhanced superconductivity can take place in a system with coexisting wide and narrow bands when the narrow band sits in the vicinity of the Fermi level. There, the Hubbard model on a two-leg ladder was studied within the fluctuation exchange (FLEX) approximation. In the two-leg ladder model, which is a two-band model with bonding and antibonding bands, one of the bands becomes wide and the other becomes narrow when diagonal hoppings are introduced. In nowadays' terminology, the narrow band in this case is incipient. Quite recently, partially motivated by studies on various lattice models with coexisting wide and narrow (or flat) bands [40][41][42][43][44][45], one of the present authors and his coworkers studied the bilayer model with diagonal interlayer hoppings [46], where one of the bands becomes wide and the other narrow, as in the two-leg ladder. There again, it has been shown using the FLEX approximation that s ± -wave superconductivity is strongly enhanced when one of the bands is incipient. The role of the finite and low energy spin fluctuations on superconductivity, along with the commonalities and differences with the two-leg ladder, has been discussed [46].
The above mentioned studies on the two-leg ladder and bilayer lattices with diagonal hoppings adopted the FLEX approximation [39,43,46], but because FLEX is basically a weak-coupling method, it is not clear whether the method can be applied to regimes where the electronelectron interaction is large. In the present study, we study Hubbard models on the two-leg ladder and bilayer square lattices, using a many-variable variational Monte-Carlo (mVMC) method [47,48], which can be considered as reliable in the strong coupling regime [49]. By compar- ing the results for the two-leg ladder (one dimensional) and the bilayer lattice (two dimensional), and with and without the diagonal hoppings, we discuss how the density of states (DOS) affects superconductivity and antiferromagnetism when one of the bands is close to being incipient.

II. MODELS, METHOD, AND DEFINITIONS OF PHYSICAL QUANTITIES
We study Hubbard models on the two-leg ladder and bilayer square lattices (Fig. 1). The Hamiltonian for the two-leg Hubbard ladder is Here c † imσ /c imσ creates/annihilates a fermion with spin σ(=↑, ↓) on the i th site on the m th chain (m=0 or 1) and n imσ = c † imσ c imσ . The nearest neighbor hoppings in the leg and rung directions are t l and t r , respectively, and the next nearest neighbor diagonal hopping is t ′ . Since t r connects two chains, we will call this the inter-chain hopping. The band structure for this model is where the case of k y = 0 (π) corresponds to the bonding (anti-bonding) band. For t ′ > 0, the bonding band is wide and the anti-bonding band is narrow.
The Hamiltonian for the square lattice bilayer Hubbard model is Here c † imσ /c imσ creates/annihilates a fermion with spin σ(=↑, ↓) on the i th site on the m th layer (m=0 or 1). The intra-layer hopping is t and the inter-layer hopping is t ⊥ , the next nearest neighbor inter-layer hopping is t ′ ⊥ . The band structure for this model is where the case of k z = 0 (π) corresponds to the bonding (anti-bonding) band. For t ′ ⊥ > 0, the bonding band is wide and the anti-bonding band is narrow.
We take N s = 60 × 2 (12 × 12 × 2) sites for the two-leg ladder (bilayer) Hubbard model with the antiperiodicperiodic boundary condition in x (y) direction. The band filling is defined as n = N e /N s where N e = miσ n imσ . Hereinafter, the site index (i, m) is simply written as i.
To study the ground state of these Hubbard models, we employ a mVMC method [47,48], which can describe the strong correlation and various ordering fluctuations accurately. Our variational wave function is defined as where P G , P J are the Gutzwiller and Jastrow correlation factors, respectively. The Gutzwiller factor punishes the double occupation of electrons defined as The Jastrow factor is defined as where n i = σ n iσ . The long-range part of this factor drives the distinction between the metal and insulator. |φ pair is the one-body part defined as where f ij is assumed to have 2 × 2 (2 × 2 × 2) sublattice structure or equivalently 2 × 2 × N s (2 × 2 × 2 × N s ) independent variational parameters for one-body part in the two-leg ladder (bilayer) systems.
To study a possible superconducting state, we consider the following BCS wave function, with where ξ(k) = ε(k) − µ, µ is the chemical potential and ∆(k) is the superconducting gap. The BCS wave function is rewritten in the real space representation as follows: In this study, we employ the BCS partial d (s ± )-wave superconducting state as the initial states for the ladder (bilayer) system, namely, ∆(k) = ∆ 0 cos k y (∆ 0 cos k z ). The variational parameters are simultaneously optimized to minimize the variational energy by using the stochastic reconfiguration method [50].
To investigate the ground state properties of these Hubbard models, we calculate the momentum distribution function and spin-structure factor, equal-time superconducting correlation. The momentum distribution function is defined as and the spin-structure factor is defined as Further, the equal-time superconducting correlations are defined as Superconducting order parameters ∆ α (r i ) are defined as Here f α (r) is the form factor that describes the symmetry of the superconductivity. For the partial d-wave superconductivity in two-leg ladder systems, we define where δ ij denotes the Kronecker's delta. For the s ± -wave superconductivity in bilayer systems, we define f s ± (r x , r y , r z ) = δ rx,0 δ ry,0 δ rz,1 .
To reduce stochastic errors, we calculate long-range average of the superconducting correlation, which is defined as where r max is 30 (6 √ 2) for the present two-leg ladder (bilayer) models. M is the number of vectors satisfying 2 < r < r max . Here, we eliminate the short range part of the superconducting correlation since it does not reflect the off-diagonal ordering nature of superconductivity in order to reduce the effect of the boundary condition.

A. Two-leg Hubbard ladder
We begin with the two-leg ladder. Figure 2a shows the inter-chain hopping dependence of several physical properties for t ′ /t l = 0 and U/t l = 4; peak value of the spin structure factor S(q max ), which is the square of the anti-ferromagnetic ordered moment, and average value of superconducting correlation P d at long distance with the partial d symmetry, corresponding to the square of the superconducting order parameter. We also plot the momentum distribution function at the anti-bonding band minimum n(0, π), which monitors whether or not the anti-bonding band intersects the Fermi level. For 1.2 ≤ t r /t l ≤ 1.6, n(0, π) decreases rapidly and S(q max ) is strongly suppressed as t r /t l increases. Therefore, the incipient-band regime is estimated to be in a range of 1.2 ≤ t r /t l ≤ 1.6. In the incipient-band regime, P d is maximized. Further, the partial d-wave superconducting phase exhibits the dome structure around t r /t l ∼ 1.5, which is consistent with a previous FLEX study on the two-leg ladder Hubbard model without t ′ [51]. For t ′ /t l = 0.4, where there are wide and narrow bands, the inter-chain hopping dependence of several physical properties are similar to those for t ′ /t l = 0 as shown in Fig.  2b. For a larger interaction value of U/t l = 8, the variation of n(0, π) against t r /t l becomes broad due to correlation effects as shown in Figs. 2c and 2d. On the other hand, we find a clear suppression of S(q max ), which indicates the Lifshitz transition. Thus, superconductivity is optimized when a band becomes incipient also in the strongly correlated regime.
We also study t r /t l dependence of the superconducting correlation P d for various values of U/t l as shown in Fig. 3. Both for t ′ /t l = 0 and t ′ /t l = 0.4, the regime of enhanced superconducting correlation extends to smaller t r /t l as U/t l increases, presumably due to band narrowing caused by U . In general, Hubbard U can narrow bands near the Fermi level and induce the Lifshitz transition [52].  2. (color online). Inter-chain hopping tr/t l dependence of the averaged partial d-wave superconducting correlation P d and the peak value of the spin structure factor S(qmax) (upper panels), the momentum distribution function of the anti-bonding band minimum n(0, π) (lower panels) for the two-leg ladder model with (a) t ′ /t l = 0 and U/t l = 4, (b) t ′ /t l = 0.4 and U/t l = 4, (c) t ′ /t l = 0 and U/t l = 8, (d) t ′ /t l = 0.4 and U/t l = 8. The band filling is n = 0.97. The yellow region denotes the incipient-band regime. In the present plots and the plots in the later figures, the error bars indicate the estimated statistical errors of the Monte Carlo sampling.

B. Bilayer Hubbard model
We next move on to the bilayer model. Figure 4a shows the inter-layer hopping dependence of physical properties for t ′ ⊥ /t = 0 and U/t = 8; the peak value of the spin structure factor S(q max ) and the average value of the superconducting correlation P s ± at long distances with the s ± symmetry. We also plot the momentum distribution function at the anti-bonding band minimum n(0, 0, π). For t ⊥ /t > 1.8, n(0, 0, π) decreases steeply, and S(q max ) is strongly suppressed as t ⊥ /t increases. Thus, the incipient-band regime is estimated to be in Inter-layer hopping t ⊥ /t dependence of the averaged s ± -wave superconducting correlation P s ± and the peak value of the spin structure factor S(qmax) (upper panels), the momentum distribution function of the anti-bonding minimum n(0, 0, π) (lower panels) for the bilayer model with (a) t ′ ⊥ /t = 0 and U/t = 8, (b) t ′ ⊥ /t = 0.6 and U/t = 8. The band filing is n = 0.94. a range of 1.8 ≤ t ⊥ /t ≤ 2.4. Around the incipientband regime, P s ± is enhanced. The s ± -wave superconducting correlation exhibits a dome structure around t ⊥ /t ∼ 2.0 [53], which is consistent with FLEX [30] and fRG [34], DCA [35] studies. For t ′ ⊥ /t = 0.6, t ⊥ /t dependence of the physical properties are basically similar to those for t ′ ⊥ /t = 0 as shown in Fig. 4b.
As in the two-leg ladder, we also study t ⊥ /t dependence of superconducting correlation P s ± for various values of U/t as shown in Fig. 5. Both for t ′ ⊥ /t = 0 and t ′ ⊥ /t = 0.6, the regime where the superconducting correlation develops extends to lower t ⊥ /t as U/t increases in the same way as the two-leg Hubbard ladder.
We note that n(0, 0, π) of the U = 8t bilayer Hubbard model varies steeper than n(0, π) of the U = 8t l two-leg Hubbard ladder, and actually resembles that of U = 4t l ladder. Since the broadness of the momentum distribution variation around the Lifshitz transition presumably originates from the correlation effect, the present result indicates the strength of the electron correlation is roughly determined by U/W , where W is the band width.

IV. DISCUSSION
The present results show that superconductivity is enhanced in the incipient band regime regardless of whether the system is one-or two-dimensional, or whether one of the bare bands is narrow or not. This is in fact consistent with the recent FLEX studies [43,46,54]. In this section, we further discuss the relation between the bilayer and two-leg ladder models, based on observations made in previous studies. Refs. [36,37,46] pointed out the important role of the finite-energy spin-fluctuations played in the enhancement of superconductivity in the bilayer Hubbard model. Especially, quite recently, in ref. [46], the role played by the spin fluctuations in different energy ranges in two-band systems has been discussed as follows. There is a critical frequency ω c in the spin- fluctuation spectrum, which should be smaller than a pairing cutoff energy ε c . In general, the low-energy spinfluctuation with ω < ω c leads to strong renormalization and hence are "pair breaking" whilst the finite-energy spin fluctuations with ω > ω c enhance T c [55]. Thus, when the low-lying spin fluctuations are suppressed while the finite energy spin fluctuations are enhanced, superconductivity can be enhanced. In multi-band systems, as one of the bands moves away from the Fermi level, the spin-fluctuation spectral weight is transferred to higher energies. When the spin-fluctuation spectral weight is away from the critical frequency of spin-fluctuations, but is within the paring cutoff energy (ε c > ε > ω c ), the pairing interaction can be strong without strongly renor-malizing the quasiparticles. On the other hand, when the spin fluctuations are concentrated at very low or too high energies, superconductivity is degraded. From this viewpoint, we further discuss the relation between the superconducting correlation obtained by mVMC and the shape of the DOS of the antibonding band.
One difference between the bilayer and the two-leg ladder observed in the present study is the U dependence of superconductivity. Namely, in the two-leg ladder, the superconducting correlation is enhanced even for U/t l = 2 when the antibonding band is incipient, while such an enhancement is not obtained for the bilayer model for U/t = 4. Note that here we compare two cases where U normalized by the bare band width are the same. Actually, a similar result was obtained in the recent FLEX calculation [46]. There, it has been pointed out that in the bilayer model, the correlation effect reduces the width of the incipient band [52], which makes more spin fluctuation weight lie within the energy regime effective for pairing. In the two-leg ladder, such effect is not necessary when the antibonding band is incipient because in a one-dimensional system, DOS is diverging at the band edge (see Fig. 6).
Another point that we notice, if we look closely, is that in the two-leg ladder, the superconducting correlation, maximized around the incipient-band regime, is reduced as t r /t l decreases and the antibonding band intersects the Fermi level, but does so rather mildly and smoothly especially for U = 8t l , whereas the reduction of the superconducting correlation in the bilayer model upon reducing t ⊥ occurs rapidly after the antibonding band intersects the Fermi level. If we compare in more detail the two cases for the bilayer model, the reduction of the superconducting correlation is more abrupt for the case with t ′ ⊥ = 0.6. These differences can again be understood from the shape of the DOS (see Figs. 6 and 7). Namely, in the two-leg ladder, where the DOS at the band edge is diverging, the DOS at the Fermi level decreases (rapidly, especially for large U because the band width shrinks due to renormalization) as t r is reduced after the antibonding band intersects the Fermi level, so that the pair-breaking lowenergy spin fluctuations become weaker. By contrast, in the bilayer model, where the DOS is diverging around the middle of the antibonding band, the DOS at the Fermi level increases upon reducing t ⊥ after the antibonding Fermi surface is formed, resulting in an increase of the pair-breaking spin fluctuations. The diverging DOS of the antibonding band approaches the Fermi level "faster" when t ′ ⊥ is finite, so that the superconducting correlation is rapidly suppressed for the case of t ′ ⊥ = 0.6 as t ⊥ is reduced. A similar analysis has been performed in ref. [46], not for the t r , t ⊥ variation, but for the t ′ (t ′ ⊥ ) variation of superconductivity.

V. SUMMARY
To summarize, we have studied superconductivity in the Hubbard model on the two-leg ladder and bilayer square lattices. In both systems, superconductivity can be optimized in a region around the Lifshitz point, where one of the bands is (nearly) incipient. The parameter dependence of the superconducting correlation function is reminiscent of the FLEX results obtained in ref. [46], which supports the scenario that superconductivity is enhanced by an incipient band due to the suppression of the near-zero-energy spin fluctuations and enhanced finite energy spin fluctuations working as an effective pairing glue. We stress that the consistency between the two approaches is highly nontrivial because the two approaches are totally different; FLEX is based on a weak coupling perturbational theory, which takes into account the spin fluctuations (in momentum space) explicitly in the effective interaction, whereas the present mVMC method takes into account the electron correlation effect in a realspace-based manner, which is expected to be more appropriate in the strong coupling regime. Since it has been shown that incipient bands enhance superconductivity in other models [40][41][42][43][44][45]52], it is an interesting future problem to study those models using the mVMC method.