Superconducting pairing from repulsive interactions of fermions in a flat-band system

Fermion systems with flat bands can boost superconductivity by enhancing the density of states at the Fermi level. We use quasiexact numerical methods to show that repulsive interactions between spinless fermions in a one-dimensional (1D) flat-band system, the Creutz ladder, give a finite pairing energy that increases with repulsion, though charge quasi-order (QO) remains dominant. Adding an attractive component shifts the balance in favor of superconductivity and the interplay of two flat bands further yields a remarkable enhancement of superconductivity, well outside of known paradigms for 1D fermions.

Systems with flat-band dispersions have recently attracted much theoretical and experimental attention [1][2][3][4][5][6][7][8][9][10][11][12][13], because when the kinetic energy of quantum particles is near zero at some or all momenta, interactions become completely decisive for the many-body ground state. In particular, a major line of work has been devoted to the enhancement of superconductivity in electronic systems with flat bands near the Fermi level [6][7][8][9][10][11][12][13][14]. While significant boost to superconductivity is generally found in flat-band systems due to the diverging density of states, previous works have predominantly inserted the electron pairing, which is at the root of superconductivity, just by hand, i.e., only using an assumed effective attractive interaction between electrons. However, microscopic interactions are primarily repulsive between electrons. In non-flat-band systems, electron pairing, and subsequently superconductivity, might emerge from repulsive interactions, but charge or magnetic orders are also closely competing, making reliable results on unconventional and high-T c superconductivity extremely challenging to obtain [15].
How the competition between superconducting and charge or magnetic orders plays out in flat-band systems with repulsive electron interactions is an almost unexplored area, which the present work addresses. Meanfield results have shown that superconductivity is more robust than charge orders when doping away from the flat band [11], further motivating us to study the competition between different orders using the most reliable approaches possible. For this purpose, 1D systems are ideal targets, as density-matrix renormalization group (DMRG) techniques [16][17][18][19] yield quasiexact results. It also allows us to study the effect of flat bands in 1D systems specifically, a burgeoning field in itself [4,14,[20][21][22]. In particular, flat-band systems might fall outside the widely used low-energy Luttinger liquid (LL) description of 1D fermions, as the necessary linearization of the dispersion breaks down [23].
In this work, we focus on a concrete 1D flat-band fermion system, the Creutz ladder [24] (depicted in Fig. 1(a)), as it has been realized experimentally using for V1/∆ 1/2, with phase (charge) QO for lower (higher) densities marked by three points. Around V1/∆ 1/2 an IP with highly enhanced superconducting correlations exists for n ≤ 0.73. Larger interactions give phase separation. a parametric cavity for bosons [2] and an optical lattice for fermions [25], and for which the physics of repulsively interacting bosons has recently also been investigated [4,[20][21][22]. The Creutz ladder has two flat bands separated by a finite energy gap, and by employing DMRG numerics, we obtain the many-body ground state in the presence of repulsive interactions between spinless fermions. We find that an itinerant-paired (IP) phase of fermions exists for a wide range of densities, with a pairing energy that even increases monotonically with repulsion. At the same time, we find that pure repulsion gives a dominant charge QO. To achieve dominant superconducting phase QO, we tune the interactions by adding some attraction. Here, in the weak in-teraction regime with the physics described by a single flat band, we realize a previously proposed state of a hole-based superconductor [26], forming a pair Luttinger liquid (PLL) [23,26]. Most remarkably, with increasing interaction strength, we find that the two flat bands in conjunction lead to a massive boost in superconducting correlations, far beyond known 1D paradigms.
Setup.-We study interacting spinless fermions on the Creutz ladder, comprised of sites a j and b j in the j-th unit cell, with annihilation operatorsf a,j andf b,j . The Hamiltonian we focus on iŝ +V 1 n b,jnb,j+1 +n a,jna,j+1 wheren a/b,j =f † a/b,jf a/b,j , V 1 denotes the nearest-neighbor (NN) interaction along the legs, and V 2 next-nearest-neighbor (NNN) interaction along the diagonals. We set the overall energy scale by fixing t = 1.
The flat-band dispersion becomes explicit by rewriting Eq. (1) in the basis diagonalizing the kinetic energy, using the maximally localized Wannier basis [4,20,21] In this basis, the kinetic energy terms give two flat bands at energies ± = ±2, with a gap ∆ = 4. Focusing on the lower band (upper band behaves similarly) and including interactions, we arrive at the single-band projected Hamiltonian where , and n −,j =ĉ † −,jĉ−,j . This representation makes it explicit that in the single-band limit the physics is that of densityassisted tunneling, equivalent to tunneling of NN pairs of fermions, with a purely interaction-set amplitude t p , plus remnant short-range interactions U 1 and U 2 .
In terms of interactions, we first study a purely repulsive model (PRM), where we set V 1 > 0, V 2 = 0, generating U 1 , U 2 > 0 in the single-band limit. Later, we expand our focus by adding an attractive interaction V 2 < 0, resulting in U 1 = U 2 = 0, which we name the mixed interaction model (MIM). We solve the full Hamiltonian, Eq. (1), numerically using two DMRG codes [27,28] to quasiexactly obtain the many-body ground states. Here, we exploit the conserved particle number, such that our results are presented as a function of density n = N/L x , where N is the total number of fermions and L x is the ladder's linear size, restricted to 0 < n ≤ 1 due to particle-hole symmetry of the Hamiltonian. Observables, apart from gap extrapolations, are computed at L x = 100 with the default bond dimension χ = 512, but we ascertained convergence using up to χ = 1024. Studying the single-band projected Hamiltonian, Eq. (3), alongside the quasiexact DMRG results, provides us with further insights.
To determine the characteristics of the many-body ground states, we first show that single-particle excitations are gapped and that low-lying excitations consist of fermion pairs by computing the one-and two-particle where E 0 (N ) denotes the ground-state energy of N fermions. These gaps are computed at different L x and then extrapolated to 1/L x → 0. We also obtain important information by focusing on the lower band. Here, we study both the phase (ph) and charge-density (ch) correlation functions particularly extracting their power-law decay exponents α −,ph/ch to determine if phase or charge QO is slowest decaying and thus dominant. We also compute the singleparticle density matrix of the lower band to verify that its single-particle excitations are gapped. As a further marker of fermion pairing [29][30][31][32], we calculate the Fourier transform of the fluctuations of the local density, FT [δ n −,j ]. In the case of strong interactions, we additionally find it useful to analyze the densitydensity correlation function across the diagonal links of the ladder in the original two-band basis Finally, for complementary information, we obtain the bipartite entanglement entropy S(l) = −Tr [ρ(l) ln ρ(l)] for 1 ≤ l < 2L x , where ρ(l) is the reduced density matrix of the first l sites [33]. PRM results.-The complete many-body phase diagram for the purely repulsive model (PRM) with V 1 > 0 and V 2 = 0 is summarized in Fig. 1(b), with a large part of the phase diagram covered by an IP phase. This is consistent with the single-band Hamiltonian, Eq. (3), displaying explicit pairing physics through the non-zero t p term.
(a) Single-particle gap ∆1 as a function of interaction V1 for various densities n in the IP phase. Inset shows same data with the role of n and V1 interchanged, with V1 increasing in direction of the arrow. (b) Power-law decay exponents α −,ph(ch) for phase (charge) correlations as a function of density n in the IP phase at V1 = ∆ /2.
Specifically, we find that by obtaining the ground state of the full Hamiltonian Eq. (1), densities n ≥ 2/5 give ∆ 1 > 0 and G − (r) decaying exponentially with distance r. Thus, single-particle excitations play no role, in strong contrast to repulsively interacting spinless fermions on ladders with non-flat bands [23]. Even more remarkably, the pairing energy ∆ 1 increases approximately linearly with the repulsion V 1 across the entire range ∆ /16 ≤ V 1 ≤ ∆ /2, as summarized in Fig. 2(a). This is also in contrast to the behavior of other 1D spinless fermion systems, such as the two-leg Hubbard ladder [23]. We further find ∆ 2 = 0 in the thermodynamic limit and that the peak of FT [δ n −,j ] occurs at momentum k = 2π(n/2), indicating that the low-energy degrees of freedom are effectively halved and thus those of fermion pairs. These results clearly indicate the existence of an IP phase with fermion pairs, and not any other higher-order bound states, at densities n ≥ 2/5.
Within the IP phase, we find that the strong repulsion results in a dominant charge QO in the lower band for all n ≥ 2/5 . This is quantified in Fig. 2(b), where we plot the power-law decay exponents α −,ph(ch) for phase (charge) correlations and find the charge-density correlations decaying slowest. These results are obtained at V 1 = ∆ /2, but are virtually identical at lower V 1 values. As a consequence, for interactions up to V 1 = ∆ /2, the only energy scale in the system is V 1 , as also evident in the projected Hamiltonian Eq. (3) and in the linear scaling of ∆ 1 in Fig. 2(a). Spot checks at V 1 = ∆ show strong deviations from these lower V 1 results (not included in Fig. 1(b)).
To check whether the IP phase is a LL described by a conformal field theory (CFT) with central charge c = 1 [23,34], we study the entanglement entropy. We find significant fluctuations, which cannot be captured by the Calabrese-Cardy formula [35][36][37]. Still, performing a naive fit yields a central charge either with a large uncertainty or clearly different from c = 1. Thus, we cannot identify the IP phase as a PLL (see the Supplemental Material (SM) [38] for details).
Complementing the IP phase, we find for lower densities a fully or partially localized flat-band behavior that is not sensitive to the repulsive interactions. At n ≤ 1/3, we find that the many-body ground state is completely localized. We ascertain this by running ground-state searches with multiple different initial conditions, falling into three categories: (1) initially localized particles, (2) random initial matrix-product states, and (3) adiabatic transition of Hamiltonian parameters from a closed ring to the Creutz ladder. All three approaches result in different, but completely localized states with the same Wannier orbital energy E 0 = nL x − for n ≤ 1/3, while all three yield the same state for n > 2/5 [39]. This numerical finding is in agreement with the analytical observation of the projected Hamiltonian Eq. (3) that a pair of adjacent fermions with center-of-mass momentum k has an effective dispersion U 1 − 2t p cos k. Thus, with U 1 = 2t p for the PRM, a pair always adds a positive amount of energy to the system. For densities n < 1/3, it is possible to keep this kinetic energy cost to zero by not creating any fermion pairs, while simultaneously making the remnant interaction U 2 zero by putting at most one fermion in every three unit cells. We therefore designate this phase as flat band in Fig. 1(b).
Finally, at n = 1/3, our numerics show a first-order transition into a coexistence regime between the flat band and the IP phases for 1/3 < n < 2/5. At these densities, placing two fermions adjacently cannot be fully avoided, and, as a consequence, it is favorable to create some itinerant pairs. At the same time, low-density patches can still exist intermittently, where both repulsive interactions and pair formation are avoided, resulting in fully localized fermions, overall yielding a mixture of IP and FB phases.
MIM results.-Since we find an IP phase with dominant charge QO within the PRM, we next seek to enhance the superconducting correlations by reducing the remnant single-band interactions U 1,2 . We do this by additionally switching on an attractive interaction V 2 across the diagonals of the Creutz ladder. In particular we investigate the case where U 1,2 are tuned to zero by setting V 2 = −V 1 , which we designate as the mixed interaction model (MIM). The complete many-body phase diagram is summarized in Fig. 1(c), capturing the behavior for both V 1 below as well as when comparable with the gap ∆ . Starting with small interactions, we note that for V 2 = −V 1 , only the pair-hopping survives in the singleband projected Hamiltonian Eq. (3). We additionally note that by performing a particle-hole transformation,ĉ −,j →ĉ † −,j , Eq. (3) reduces to a model recently solved exactly by mapping to a system of hard-core bosons representing two-fermion composites [26]. For V 1 = ∆ /8, ∆ /4, we find that the ground-state energies for the full system (Eq. (1)) numerically agree very well with the energies of this analytical solution at all densities 0 < n ≤ 1 (see SM [38]). Furthermore, based on the exact solution [26] and its LL description [23], G −,ph (r) is expected to scale as ∼ r −1/(2K) . Here K is the Luttinger parameter, equal to 1 in the limit of n → 0 and approaching 1/4 as n → 1, with K ≥ 1/2 giving superconducting phase QO [23]. Our numerics matches this behavior very closely, finding dominant superconducting QO for all densities where K ≥ 1/2, as marked in Fig. 1(c). In addition, we find a finite value for ∆ 1 , an exponential decay of G − (r), and a peak of FT [δ n −,j ] at momentum k = 2π(n/2), all further confirming that the MIM realizes a PLL in the small interaction, single-band limit.
The behavior changes drastically in a regime around V 1 ∆ /2 where we start seeing the upper flat band affecting the physics. Our DMRG numerics find the occupation of the upper band, previously reaching at most 1.3% and typically much less, now boosted to 3.5-8.5%. Most interestingly, we find for n ≤ 0.73 both a strong enhancement of the superconducting correlations and the effective physics being very different from the PLL found at lower V 1 . For example, while still decaying largely with a power law, G −,ph (r) now shows an extremely slow decay, especially at smaller densities as illustrated in Fig. 3(a) (see SM for V 1 dependence [38]). This decay is much slower than typically found for interacting fermions in 1D [23]. We also find the entanglement entropy to be very different compared to smaller V 1 , as shown in Fig. 3(b), saturating to a constant value inside the bulk. Moreover, for n ≤ 0.47 we find no peak at finite momentum in FT [δ n −,j ] due to complete vanishing of real-space density fluctuations. A peak at k = 2π(n/2), indicating pairing, only appears suddenly for n > 0.47. At the same time, we find ∆ 1 > 0 as expected, and ∆ 2 = 0, showing that the attractive component of the interaction does not give rise to multi-particle bound states beyond pairs.
To better understand the ground state in the regime around V 1 ∆ /2, we analyze the correlation function G diag (r), as illustrated in Fig. 3(c). While we find the expected short-range correlation hole for composite hard-core bosons for V 1 = ∆ /8, ∆ /4, the behavior at V 1 = ∆ /2 is very different, with a much deeper and wider correlation hole, explaining why we see no higherorder bound states beyond pairs. There is now also a positive correlation for finding pairs at distance r = 2 which can explain the major boost in superconducting correlations. We interpret these results as the upper flat band providing a strong additional short-range attraction for pairs in the lower flat band through virtual processes, rendering a substantial boost to superconductivity, well beyond the known PLL regime for pairing in 1D fermion systems.
While superconducting correlations are boosted for V 1 ∆ /2 in a broad density range, we eventually find charge QO dominating for n > 0.73. In this highdensity regime the entanglement entropy approaches the Calabrese-Cardy pattern, although never quite obtaining it, pointing to significant remnants of non-LL physics. Finally, for strong interactions V 1 > ∆ /2, we find phase separation where also higher-order bound states are formed.
Discussion.-Using quasiexact many-body methods, we show how flat bands give rise to fermion pairing from purely repulsive interactions, with pairing strength even increasing with repulsion. Although we find chargedensity QO ultimately dominating for pure repulsion, it would be pre-mature to conclude that 2D or 3D ar-rays formed from Creutz ladders will necessarily have a charged-ordered ground state, as weak inter-ladder coupling could help stabilize the superconducting phase. Our work also highlights the great potential of multiple flat bands to strongly stabilize superconductivity, establishing a regime beyond known paradigms for pairing in 1D systems. As Creutz ladders have been realized with neutral ultracold fermionic atoms via optical potentials [25], our results might be tested experimentally, e.g., using dipolar ultracold gases in the near future.
to extract the central charge c from our data. Here, l is the length of the subsystem, L = 2L x is the size of the system, and S 0 is a model-dependent constant. For the purpose of fitting, we also sometimes use an offset l 0 , such that l ∈ [1 + l 0 , Lx]. In Fig. S1, we plot the entanglement entropy as a function of subsystem size l for the PRM in black circles. The entanglement entropy shows strong oscillations with several different wavelengths, which do not fit known Calabrese-Cardy formulas that incorporate oscillations, for example due to open boundary conditions [37]. A naive fit using Eq. (S1) is shown with solid blue line with the extracted central charge given in each subfigure. In Fig. S2, we show the extracted central charge c as a function of density n using several different offsets l 0 . We find that the value of c is often both offset-dependent and suffers from large errors due to the strong oscillations in the entanglement entropy. Usually, whenever the Calabrese-Cardy formula works well, the effect of an offset l 0 should be negligible. In aggregate, while over significant segments of densities c = 1 is not implausible, it is seemingly impossible to conclude so with any degree of certainty. There exist even some densities where c = 1 seems highly unlikely, such as around n = 2/3, and for n > 9/10. For these reasons, we cannot identify the IP phase of the PRM as a pair Luttinger liquid (PLL).

MIM: Accuracy of analytical solution for small interaction strengths
As we state in the main text, the single-band projected Hamiltonian is exactly solvable for the mixed interaction model (MIM) [26]. In this section we compare this exact solution with the DMRG results of the full Hamiltonian. The projected Hamiltonian of the MIM readŝ The solution in Ref. [26] provides an exact analytical expression for the ground-state energy per unit cell We compare this analytical expression with our DMRG results and plot in Fig. S3 the relative difference, δe 0 = | e DMRG 0 −e theory 0 e DMRG 0 |, as a function of density n for several different interaction strengths V 1 . For V 1 = ∆ /8, the largest δe 0 is less than 1% and for V 1 = ∆ /4 the maximum of δe 0 only rises to 3%. These small deviations justify our expectation that the single-band projected Hamiltonian, Eq. (3) in the main text and Eq. (S2) above, captures the physics very well. On the other hand, for larger interaction V 1 = ∆ /2, the deviation from the analytical result increases monotonically with density, to as high as 15%. This is not surprising as we also see a different behavior in this regime for other quantities, such as correlation functions and entanglement entropy, in the main text discussed as being due to the interplay between the two bands present in the full Creutz ladder.

MIM: Phase correlation function for different interaction strengths
In Fig. 3(a) in the main text, we present the exponent of the algebraic decay of the pair correlation function, G −,ph (r), for the MIM at V 1 = ∆ /2. To show the drastically different behavior of G −,ph (r) for this V 1 in comparison to smaller interaction strengths, we plot G −,ph (r) for several different interaction strengths at fixed density n = 0.3 in Fig. S4 . It is evident that G −,ph (r) for V 1 = ∆ /2 decays much slower than for lower values of V 1 .