Superconductivity in a topological lattice model with strong repulsion

I. INTRODUCTION

Superconductors are remarkable phases of matter wherein charge-e fermions-which microscopically repel one another-instead bind into charge-2e Cooper pairs that then condense. Traditionally, these phases are best understood in the context of weak interactions [1][2][3]. For example, in the paradigm put forth by Bardeen, Cooper, and Schrieffer (BCS), phonons mediate pairing between electrons near the Fermi surface of a weakly correlated Fermi liquid [2]. Despite its widespread applicability, the limitations of the phonon-based theory are suggested by the experimental discovery of "unconventional" superconductors, which display a plethora of exotic phenomena ranging from high transition temperatures, to nonstandard pairing symmetries, and parent states with anomalously low carrier densities [1,4]. * These authors contributed equally to this work. Such superconductors arise in surprisingly distinct settings, e.g., the cuprate compounds [5][6][7][8][9], iron pnictides [10][11][12], and possibly even in twisted bilayer graphene (TBG) [13][14][15][16][17][18] and other graphene multilayers [19][20][21][22][23][24][25]. This diversity motivates a fundamental question that is especially relevant in the present era of highly tunable solid-state and cold atom synthetic quantum materials. Namely, what microscopic ingredients aid the realization of unconventional superconducting ground states?

For the past several decades, investigations into unconventional superconductivity have brought to the forefront the ingredient of strong interactions. This is encapsulated in the celebrated Fermi-Hubbard model [26][27][28][29][30], which consists of spinful fermions occupying tightly localized orbitals coupled by strong repulsion. These minimal ingredients underlie several proposed mechanisms for unconventional superconductivity [6][7][8][9] and give rise to numerous other exotic phenomena (see, e.g., Refs. [31][32][33][34][35][36]). As a consequence, the Fermi-Hubbard paradigm has been the subject of decades of sustained interest in condensed matter physics, pushing the forefront of theory and motivating new experimental paradigms such as cold atom quantum simulation [37][38][39][40][41][42][43][44].

Despite their success, the simple ingredients in the Fermi-Hubbard model are insufficient to describe strongly interacting phenomena reliant on band topology, such as the celebrated quantum hall effects [45,46]. This setting of topological bands induced by strong magnetic fields, with accompanying time-reversal breaking, is naively unfavorable for superconductivity in the presence of repulsive interactions. However, a number of works have suggested regimes where superconductivity may emerge in topological Hofstadter-type models, such as [47][48][49][50][51]. Recently, experiments on TBG and other graphene moiré multilayers-in which topology is widely believed to play a crucial role-present a related time-reversal-invariant context [52,53] in which superconductivity can and does appear [13]. As such, various studies have sought to demonstrate the presence of an all-electronic route to superconductivity in models of TBG and, more broadly, in electronic Hamiltonians with the ingredients of time-reversal symmetry, strong repulsive interactions, and topological bands. However, owing to the complexity of the many-body problem, theoretical evidence for such a route has so far only derived from approximate or perturbative methods, e.g., mean-field theory [54,55], large-N expansions [56], and effective continuum descriptions [56][57][58]. Elucidating the existence and character of superconductivity in a fully microscopic model with these ingredients, therefore, remains an important and outstanding challenge.

In this work, we introduce a minimal model-the "double Hofstadter" model-that incorporates precisely these ingredients and employ unbiased matrix product state techniques to study its phase diagram on finite-circumference cylinders. Our study reveals that their interplay leads to a plethora of exotic phenomena and culminates in the numerical demonstration of an unconventional superconducting ground state at a finite circumference, with further evidence that this behavior may persist on larger cylinders.

II. SUMMARY OF KEY RESULTS

This work studies the double Hofstadter model, a simple lattice model with the following properties:

1. short-range repulsive interactions (akin to the Hubbard model), 2. narrow or "flat" topological bands (akin to quantum Hall and some moiré systems), 3. time-reversal symmetry (facilitating superconducting phase coherence).

To bake all of these ingredients into a single Hamiltonian, we start with the familiar Harper-Hofstadter model of spinful fermions on a square lattice with a magnetic flux = 0 /q per plaquette [59,60]. This is a lattice model whose lowest band carries a Chern number and limits to the lowest Landau level as q → ∞, but is not time-reversal symmetric. We therefore consider a bilayer Hofstadter model (Fig. 1) with opposite magnetic fields in opposite layers, whereupon the model becomes time-reversal symmetric. Finally, we allow tunneling between the layers and add strong local repulsive interactions, as in the Hubbard model. We analyze the FIG. 1. The Double Hofstadter Model. This work considers a simple lattice model with three essential ingredients: (1) narrow topological bands with (2) strong repulsive interactions and (3) timereversal symmetry. These properties are realized in a bilayer lattice model of spinful fermions made by joining together two copies of the Harper-Hofstadter model with opposite magnetic fields in opposite layers. Bonds represent electron hoppings with in-layer magnitude t and interlayer magnitude g. Peierls phases of e inπ/2 , graphically depicted via n arrows along the direction of hopping, yield a magnetic flux of ± = ±π/2 per plaquette in the top and bottom layers, respectively. Finally, there are repulsive density-density interactions (wavy lines) with magnitude U 1 onsite and U 2 between the layers. resulting double Hofstadter model and argue that it contains an unconventional p-wave superconductor of purely repulsive origin. After defining the model precisely in Sec. III A, our key results are fivefold and culminate in connections to contemporary experiments.

Flat bands and fragile topology. First, in Sec. III, we present the bandstructure of our model (see Fig. 2). Crucially, the lowest four bands are narrow and possess "fragile topology" [61][62][63][64][65]. Even though their total Chern number vanishes, fragile topology implies that no symmetric, exponentially localized Wannier basis for these bands exists (without the addition of trivial filled bands to the model). This precludes any tight-binding description for the flat band subspace consisting of symmetric and exponentially localized orbitals.

Correlated insulator at half-filling. Next, in Sec. IV, we investigate the ground state phase diagram of the double Hofstadter model at half-filling of the flat band subspace. We first demonstrate analytically that the ground state at strong in-layer coupling is a generalized quantum Hall ferromagnet, which is a spontaneously generated quantum spin Hall insulator. We refer to this state as a layer antiferromagnet ("LAF") as it is characterized by ferromagnetism within each layer but antiferromagnetism between layers [see Fig. 3(a)]. We confirm this analytical understanding numerically using infinite cylinder density matrix renormalization group (cylinder iDMRG) methods. Our simulation of the model is enabled by recent advances in the compression of matrix product operators [66].

Unconventional superconductivity. Third, in Sec. V, we provide numerical evidence that hole doping the LAF produces an unconventional superconducting ground state. In this purely repulsive setting, our evidence is based on stateof-the-art parallelized iDMRG numerics and is strongest at circumference L y = 5. There, we conclusively demonstrate both (algebraically) long-ranged superconducting correlations (a) depicts the single-particle energy bands of the model at a flux fraction 1/q = 1/4 and interlayer tunneling amplitude g = 0.1. Including spin, there are 16 total bands. In this work, we focus on the lowest four (highlighted in blue), which have narrow bandwidth and are well-isolated from the remaining bands. (b) These narrow bands, shown here for g = 0.3 in a smaller energy window, have w 2 = +1 fragile topology and possess two Dirac cones that wind in the same direction. In (c), we show a schematic of the layer-polarized basis states, which span the flat band subspace and carry opposite Chern number C = ±1 in opposite layers. When interactions are large, the energy difference between flat bands is irrelevant; in Sec. III D, we show that a maximally layer polarized basis is more natural for studying ground state physics in the low-energy subspace. and a finite charge gap. Despite the exponential scaling of simulation complexity in L y , we attempt a similar analysis on larger cylinders where we find preliminary-but ultimately inconclusive due to bond dimension limitations-evidence for a superconducting ground state. We conclude by verifying that the superconducting phase at L y = 5 is robust, with stability to several perturbations including additional repulsive interactions, and always coexists with LAF order.

In Sec. VI, we devise a novel and widely applicable method to probe the superconducting gap function in chargeconserving cylinder iDMRG, which we use to systematically characterize the putative superconductor. We find that it is well-described by the weak pairing of quasi-holes above the strongly correlated LAF insulator (i.e., it is "BCS" like) with a p-wave pairing symmetry. We emphasize that this weakpairing phenomenology emerges in a setting where electronic repulsion dominates the kinetic energy and is, moreover, unearthed using nonperturbative numerics that do not presuppose a mean-field description. This observation then enables us to estimate the superconducting transition temperature.

Optical lattice realization. We conclude by discussing connections between the double Hofstadter model and contemporary experiments in both cold atoms and the solid state. First, drawing on recent developments in the quantum control of alkaline-earth atoms in optical lattices [67][68][69][70][71][72][73][74][75][76][77], we develop a concrete experimental blueprint to realize the double Hofstadter model by encoding its degrees of freedom in the clock and hyperfine states of such atoms.

Connection to solid-state experiments. Next, we present a link between the double Hofstadter model and models of moiré materials. Explicitly, we define a mapping between the symmetries, Hamiltonian, and even some ground states of our model and those of chiral twisted bilayer graphene [78][79][80][81][82].

A. Model

The double Hofstadter model consists of spinful fermions in a bilayer square lattice (see Fig. 1). At each lattice site r = (x, y) ∈ Z 2 , we define a row vector of many-body creation operators,

We use notation ℓ ∈ {T, B} for layer and s ∈ {↑, ↓} for spin, with associated Pauli matrices ℓ 0,x,y,z and s 0,x,y,z . The Hamiltonian has a few simple components. First, it contains nearest-neighbor hopping within each layer. Second, fermions in the top (bottom) layer experience an upwards (downwards)-pointing magnetic field that induces a 2π/q flux per plaquette, where q is a fixed positive integer. This is implemented in a tight-binding model via the Peierls substitution with a choice of electromagnetic gauge field whose value in each layer is A T/B = (0, ±(2π/q)x, 0). Concretely, the in-layer hopping Hamiltonian consists of two time-reversalrelated copies of the Harper-Hofstadter model with a 1/q flux fraction:

where we use the convention that the top and bottom layers correspond to ℓ z = +1 and -1, respectively. Next, electrons may tunnel vertically between the two layers:

Finally, the Hamiltonian contains both an in-layer and interlayer Hubbard interaction, given by

where nrℓ = s nrℓs and where the colons indicate normalordering of the fermion operators. With these three ingredients, the full Hamiltonian is (see Fig. 1):

where ĥ is the noninteracting portion of the model and V consists of interaction terms. In what follows, we set the magnitude of the in-layer hopping strength to t = 1 and scale the other couplings to be in these units. Furthermore, throughout the rest of this work, we will fix the flux fraction to 1/q = 1/4 (the case depicted in Fig. 1). We remark that spinless variants of this model have been explored in previous works for both bosons and fermions [50,51,[83][84][85][86]. However, we emphasize that most of the physics revealed in this work crucially depends on the presence of both the layer and spin degrees of freedom.

B. Symmetries

We now discuss the symmetries of the model beyond charge U(1) Q and spin SU(2) S , which are manifest. We will find that the coupling of two opposite Hofstadter layers makes both the magnetic point group and the translation symmetries of the model slightly nontrivial.

We first note that the model, while not directly invariant under time reversal (which would flip the magnetic field penetrating either layer), is invariant under a modified symmetry M z T . This operator consists of the usual time-reversal operator T , that flips time and electronic spin, followed by an in-plane mirror flip M z that flips the layer. Furthermore, the point group of Ĥ can be determined by considering the model's embedding into 3D space (see Fig. 1). In particular, the model is invariant under the dihedral group of twofold rotations around all three cartesian axes (though not fourfold, see below). This group is generated by π rotations around the ẑ and x axes, which we denote by C 2z and C 2x , respectively. The latter naturally exchanges the two layers and is chosen to also flip spin. The combination of the modified time-reversal symmetry with the point group defines the "magnetic point group." As its generators, we take

where we have de ned a spacetime inversion symmetry Î = Ĉ2z • Mz T , which has the convenient property that it preserves the electronic crystal momentum, and where C 2x r = (x, -y).

We remark that the spacetime inversion symmetry is antilinear ( Îi Î-1 = -i) while the ordinary point group symmetries are linear.

Next, the translation symmetries of Ĥ are also governed by the magnetic eld. De ne the operators

These are not symmetries (unless g = 0), but will be used to de ne the magnetic translation group. Their nontrivial commutator t-1 y t-1

x ty tx = e 2πi q ℓ z encodes the layer-dependent AharonovBohm phase accrued by an electron hopping counter-clockwise around a plaquette.

The translation symmetry of the Hamiltonian is substantially reduced by the interlayer tunneling, Eq. (3). While the magnetic ux through each layer remains spatially uniform, the tunneling induces a ux pattern in the region between the layers that is nonuniform;foot_0 see the shaded vertical plaquettes shown in Fig. 1. For even q, this pattern repeats every q/2 sites,foot_1 and thus the magnetic translation symmetries of Ĥ are generated by mx = ( tx ) q/2 and my = ty , which anticommute:

To de ne the Bravais lattice and Bloch states, we consider the largest abelian subgroup of translations. Its generators

specify a Bravais lattice L with primitive lattice vectors a x = qx and a y = ŷ. The corresponding reciprocal lattice is then spanned by G x = 2π x/q and G y = 2π ŷ. Each unit cell contains q evenly spaced sites. We describe positions in the underlying square grid by specifying a unit cell and sublattice:

and relabel the fermion operators as ψ † σ ℓs (u).

Layer-polarized basis for flat bands

When the interaction strength U greatly exceeds the bandwidth W of the low-lying bands, the single-particle energy eigenbasis is no longer a physically useful way to organize the narrow-band subspace. Drawing inspiration from the g = 0 limit in which layer is a good quantum number, we instead work in a basis of states that are maximally polarized into either the top or bottom layer. Crucially, these states will carry a layer-dependent Chern number, a fact that will underpin our analytical and numerical understanding of the many-body phenomena. Concretely, we de ne a "layer-polarized" basis for the at-band subspace as an eigenbasis of the projected layer operator:

and demand that the (real) eigenvalues carry a consistent sign across the Brillouin zone. Here, lz (k) = ψ † (k)ℓ z ψ (k) is dened in terms of the Fourier transformed microscopic fermion operators,

In this basis, the creation operators for the layer-polarized orbitals take the form of Bloch-like states:

For numerical convenience, we x the residual gauge freedom by requiring that these layer-polarized orbitals be maximally localized in x when Fourier transformed, resulting in hybrid Wannier states [92][93][94][95].

We further remark that the two eigenvalues in Eq. ( 11) are opposite in sign and equal in magnitude. For all parameters of interest, their magnitude is bounded below as |λ T/B (k)| > 0.95, very near the maximum possible value of unity. Like the Harper-Hofstadter bands of the layer-decoupled Hamiltonian ĥt , the layer-polarized basis states carry a Chern number C = ±1 for l = T/B respectively. We emphasize that this holds in spite of the fact that layer is not a good quantum number in the presence of interlayer tunneling. We therefore refer to the label l as the "dressed layer" index and remark that it is indistinguishable from the true layer index ℓ, within the lowenergy subspace, in the layer-decoupled limit.

Flat band projected Hamiltonian

We now write the projected Hamiltonian in terms of the layer-polarized basis for the at-band Hilbert space:

This decomposition of the Hamiltonian is chosen so that the interaction part is positive semi-de nite:

and δh(k) consists of the remaining single-particle terms. The latter can be further decomposed as δh 0 (k) + δh g (k), where δh g (k) ∼ O(g) is purely off-diagonal in dressed layer and δh 0 (k) ∝ l 0 is controlled by the bandwidth of the bare Hofstadter bands, up to a constant offset and O(g 2 ) corrections. We refer to V as the "strong-coupling" form of the interaction. Explicitly, δ ρqℓ = ρqℓρqℓ Î denotes the layer-resolved charge density relative to a translationally invariant reference charge background that evenly lls the two layers, where (cf. Refs. [79,80])

and ρqℓ = 1

For purely on-site interactions, Eq. ( 4), the interaction kernel V (q) is a matrix in the layer indices that can be expressed as

which is crucially independent of q. Moreover, we note that the momentum q lives in an extended Brillouin zone, or "EBZ," which is de ned relative to the microscopic square lattice as opposed to the Bravais lattice (see Ref. [90] for more detail). The symmetries of the double Hofstadter model strongly constrain its form factors. Namely, spacetime inversion symmetry gives q (k) = ℓ qℓ (k) the following functional form:foot_2

with real-valued functions F S/A and S/A . In the "decoupled" limit g → 0, the U(2) symmetry of the double Hofstadter model is enlarged to U(2) × U(2), corresponding to independent spin rotations and charge conservation in each layer. This implies F A vanishes when g = 0, and in fact F A = O(g) is always perturbatively small. The decoupled limit, therefore, serves as the starting point for the strong coupling theory of the double Hofstadter model that we develop in the next section.

In-layer ferromagnetism

We start by working in the limit with vanishing interlayer tunneling and interlayer interactions, g = U 2 = 0, in which the two layers are completely decoupled. We further take the limit where the remaining in-layer dispersion δh 0 (k) [Eq. (14)] vanishes (up to an inconsequential momentum independent shift). At the lling fraction ν = 2, the in-layer interaction Hamiltonian favors states that occupy the two layers equally. All that remains, therefore, is to determine the behavior of spin in each layer. For generic narrow, half-lled Chern bands, strong repulsive interactions will generally favor ferromagnetic spin alignment in a phenomenon known as quantum Hall ferromagnetism [45,97]. Intuitively, this behavior arises from the interplay between strong interactions and the spread of Wannier functions. In a topological band, Wannier functions cannot all be atomically localized. It is therefore advantageous for the electrons to spin polarize, which forces the many-body wave function to vanish when two electrons are brought together and results in a favorable direct exchange energy. We remark that this is different from the paradigm of the Fermi-Hubbard model. There, electrons can be perfectly localized to atomic lattice sites, resulting in a vanishing direct exchange energy and a dominant antiferromagnetic "superexchange" interaction.

Consequently, the layer-decoupled double Hofstadter model should have bilayer quantum Hall ferromagnets as favorable ground states, namely,

Here, each layer is associated with an order parameter characterizing its spontaneous magnetization:

and electronic states whose spins are aligned to this axis:

where φ ℓ and θ ℓ are the azimuthal and polar angles of n ℓ , respectively. We can place this intuition on rigorous analytical footing by recognizing that the term δ ρqℓ δ ρ-qℓ appearing in the strong-coupling interaction V [Eq. ( 15)] is positive semide nite for each q. This implies that any state annihilated by every δ ρqℓ is an exact ground state of this interaction.

In particular, this is true of the ferromagnetic trial states in Eq. ( 20):

Analytically, this follows from the Pauli exclusion principle and the form factors qℓ (k) being layer-diagonal in the g = 0 limit, Eq. ( 19), that is itself a consequence of layer U(1) symmetry. This shows that the trial states (20) are exact ground states of the strong-coupling Hamiltonian (15) and are therefore approximate ground states of the full Hamiltonian.

Cylinder iDMRG numerical method

Let us now describe the cylinder iDMRG numerical method that we will use to obtain the phase diagram at ν = 2 and, upon doping, superconductivity. Since we are interested in the physics at partial lling of the isolated at bands, we restrict to matrix product state (MPS) wave functions that live entirely within these bands, thereby reducing the number of degrees of freedom from all 4q bands to just 4, for layer and spin. Moreover, since the fragile topology of the at band subspace forbids an exponentially localized basis that respects all the symmetries [61][62][63][64][65], we choose a computational basis of mixed position-momentum space states. Such states, commonly referred to as "hybrid" Wannier states [92][93][94][95]98,99], are exponentially localized along the cylinder and are delocalized around the cylinder. A similar topological restriction arises in TBG [100,101], where a line of work by some of us has established the use of such states in cylinder iDMRG [102][103][104].

This basis, which we encode into the MPS tensors, takes the following form:

and is depicted graphically in Fig. 3(b). Here, j labels "rings" along the cylinder length, k y = (2πn + y )/L y labels discrete momentum cuts through the BZ that can be shifted by tuning the ux through the cylinder, y ∈ [0, 2π ). The basis is chosen such that its Fourier transform φls (k) is a layer-polarized basis (as described in Sec. III D), permitting φls ( j, k y ) to be labeled by dressed-layer and spin. Moreover, the phase of each φls (k) is chosen such that φls ( j, k y ) are maximally localized in the x direction [90,98,105], giving hybrid Wannier states. Each ring of this computational cylinder, therefore, has a total of 2 × 2 × L y tensors whose physical legs encode the occupations of fermions at a given dressed-layer, spin, and momentum. Our numerics explicitly conserve electric charge U(1) Q , spin U(1) S ⊂ SU(2) S , and Z/L y Z momentum around the cylinder.

Carrying out iDMRG on this model requires "MPO compression" [66]. When the Hamiltonian is resolved in the basis (27), both the dispersion and interaction become long-ranged. A naive matrix product operator (MPO) representation of this Hamiltonian has bond dimension D ≈ 10foot_5 , far beyond what is computationally feasible. We therefore apply an MPO compression technique [66,102] to create a faithful MPO representation of our Hamiltonian with bond dimension D ≈ 300 and operator norm errors of δ < 10 -3 t or better. We consider cylinder circumferences L y = 5, 6, 7. Since the MPS bond dimension required to accurately describe ground states increases exponentially with L y , we employ parallelized iDMRG simulations to reach MPS bond dimensions as large as χ = 12288 using the open source tenpy library [106].

The LAF at strong in-layer interactions

We now numerically verify the presence of LAF order at ν = 2 in the regime of large in-layer repulsive interactions U 1 ≫ t. First, however, we must address a subtlety associated with the cylinder geometry used for iDMRG. Since the underlying in nite cylinder is quasi-1D and the candidate LAF state spontaneously breaks continuous spin rotational symmetry, Hohenberg-Mermin-Wagner (HMW) considerations imply that it must disorder due to quantum uctuations [107][108][109]. The resulting quasi-1D state is symmetric with a spin gap that decays exponentially in the cylinder circumference [110]. Moreover, its LAF polarization, though nite at nite bond dimension, decays to zero as χ → ∞. To estimate the extent of the LAF phase in the 2D limit, we take an approach analogous to a susceptibility experiment. In particular, we add a weak Zeeman eld to the Hamiltonian:

and take B z = 10 -2 , which is much smaller than all other microscopic scales in the problem. Such a Zeeman eld explicitly breaks SU(2) S spin symmetry down to U(1) S but leaves all other symmetries unbroken. 5 We then approximately identify the phase extent of the LAF in the 2D limit by evaluating

in iDMRG, which effectively measures the susceptibility to LAF ordering. Figure 3(c) shows n z LAF as a function of g and U 1 at xed U 2 = 0, with a heuristic phase boundary identied where ∂n z LAF /∂U 1 has the largest slope. In Ref. [90], we numerically verify that the LAF region of the phase diagram is insulating by computing the charge-1e correlation length, ξ 1e , and demonstrating that it converges to a nite value as χ → ∞. In summary, we have numerically identi ed the presence of a LAF insulating state at ν = 2 that inhabits a broad portion of the strong-coupling regime, con rming our analytical expectations.

Competing phases at weak interactions

Next, we explore the weakly interacting regime. We rst remark that the noninteracting band structure at ν = 2 exhibits two distinct phases as a function of g. The rst occurs at g ≪ 1, where the nearly degenerate, weakly hybridized Hofstadter bands give rise to a metal with a one-dimensional Fermi surface. In contrast, suf ciently large g gaps out these bands except at the two Dirac cones protected by I and SU(2) S (see Sec. III C), yielding a Dirac semimetal (DSM).

Using iDMRG, we con rm that these phases persist in the presence of weak interactions. We diagnose the metallic phase by plotting the fermionic occupations

Since we work at nite cylinder circumference, we access only a nite number of cuts through the BZ. Nonetheless, we can gain continuous momentum resolution by threading ux through our cylinder [111], as depicted in Fig. 3(b). The results of this are shown in Fig. 3(d) for g = 0.15 and U 1 = 0.05. We nd that the occupancies obtained using iDMRG quantitatively match the band structure occupancies of the noninteracting limit.

Next, we con rm the DSM by two independent means. 6 A de ning characteristic of this state is that it possesses two Dirac cones that wind in the same direction within each spin sector. This can be detected via the winding of arg P TB (k), where P is the correlation matrix [102]:

Figure 3(c) plots arg P TB (k) over the ux-threaded BZ at g = 0.05 with nite interactions U 1 = 0.05. We nd that it indeed has the same winding around the points k ± = (π/4, ±π/2). That these points are Dirac cones-even at nite interaction strength-can be con rmed by computing the maximum charge-1e correlation length ξ 1e (k y ) as a function of the momentum k y around the cylinder, made continuous by threading ux. In a weakly interacting system, this quantity parameterizes the decay of correlations of the form ⟨ φ †

x,k y φ0,k y ⟩ ∼ e -x/ξ 1e (k y ) . As such, when charge-1e excitations are gapless, ξ 1e should diverge. Indeed, we nd that it strongly peaks and grows rapidly with bond dimension precisely when k y passes through the isolated values ±π/2, indicative of a semimetallic gap closure [see Fig. 3(e)].

Superconductivity in cylinder iDMRG

Let us discuss how superconductivity manifests in in general in cylinder iDMRG (see also [57,112]). First, since iDMRG uses a quasi-1D cylinder geometry with nite circumference L y , we must understand how to diagnose the quasi-1D analog of the superconductor-the Luther-Emery liquid-and characterize its approach to a superconductor as L y → ∞. The Luther-Emery liquid is characterized by the presence of fractionalized excitations-namely gapless "chargeons" and gapped "spinons"-that independently carry the electron's charge and spin (see, e.g., Ref. [113]). As a consequence, charge-2e excitations (with vanishing total Ŝz ) have correlations that exhibit algebraic long-range order:

As the cylinder circumference is increased, the realization of true long-range superconducting order in 2D is tied to the behavior η(L y ) → 0 as L y → ∞. In the Luther-Emery liquid, furthermore, correlation functions of electrons (gapped states comprised of both the chargeon and spinon) exhibit exponential decay:

At nite circumference, this behavior crucially distinguishes the Luther-Emery liquid from the Luttinger liquid, where such correlations are algebraic.

A second important subtlety concerns the nite bond dimension of the MPS, which only permits the expression of exponentially decaying (connected) correlations, and introduces a length scale ξ χ associated to the largest nontrivial eigenvalue of the MPS transfer matrix. As a consequence, power-law correlations of any charge-Q operator are cut off at a distance ξ Q,χ ξ χ that scales with bond dimension according to Refs. [114][115][116]:

However, in practice, this scaling behavior can also be found in gapped systems when the true correlation length ξ , which dictates the exponential decay and is inversely related to the energy gap, greatly exceeds ξ χ . This makes it dif cult to distinguish the Luttinger and Luther-Emery liquids discussed above, especially in scenarios where the latter has a very small spinon gap, as both may appear to show algebraic scaling of ξ 1e,χ . Nevertheless, even in such cases, the two could be distinguished by making use of standard arguments in the nite entanglement scaling theory of matrix product states [114][115][116]. In particular, when the MPS is attempting to capture a fully gapless state like the Luttinger liquid, ξ χ is the only length scale present in the correlations of the state. As such, we expect correlation functions of charge Q operators will decay with a correlation length ξ Q,χ = a Q ξ χ , where a Q (at asymptotically large χ) is a constant of proportionality that does not depend on bond dimension. For the Luttinger liquid, this leads to the sharp prediction that ξ 2e /ξ 1e approaches a constant as χ → ∞. This contrasts with the Luther-Emery liquid, in which ξ 2e /ξ 1e is expected to diverge.

1e: ?

FIG. 4. Evidence for superconductivity. Working at fixed parameters (g, U 1 , U 2 ) = (0.15, 1.5, 0) and doping ν = 2 -1/L y , we provide our numerical evidence for the existence of a superconducting state at L y = 5 and algebraically decaying pair correlations at all three circumferences. In (a), we examine the superconducting correlation function [Eq. (38)] as a function of bond dimension χ for each circumference, finding that they each approach a power law as χ → ∞. The power law obtained via a scaling collapse [see panel (b)] is shown with an orange dashed line. For ease of viewing, the correlations for L y = 6, 7 are offset vertically by an arbitrary amount. (b) Power law behavior is quantitatively demonstrated by performing a scaling collapse [Eq. ( 39)], shown here for L y = 5. To differentiate this power-law behavior from that of a gapless metallic state, in the inset we show that the ratio ξ 2e /ξ 1e increases as a function of bond dimension, a hallmark of a superconducting state [See Sec. V A 1]. (c) Power-law exponents of the superconducting correlations for all three channels as a function of cylinder circumference L y . The powers η(L y ) decrease as L y increases, consistent with the survival of superconductivity in the 2D limit. At the bottom of the panel, we summarize our confidence in the requisite gaplessness of 2e excitations and gap to 1e excitations at each circumference (see Sec. V A 3 for more details).

In summary, to demonstrate the existence of superconductivity in cylinder iDMRG, one must in principle perform a double scaling analysis. At each fixed circumference, the state must approach a Luther-Emery liquid, characterized by algebraic correlations of charge-2e operators and the divergence of ξ 2e /ξ 1e with bond dimension. In addition, in principle, one must show that the exponent of algebraic decay for charge-2e operators goes as η(L y ) → 0 as L y → ∞. However, since increasing the cylinder circumference results in exponential scaling of the simulation complexity, cylinder iDMRG studies can typically only demonstrate convergence to a Luther-Emery liquid for a few, relatively small circumferences [117][118][119][120][121].

Bond dimension scaling

We begin by providing strong numerical evidence for Luther-Emery liquid physics at L y = 5 by performing a careful scaling in the bond dimension. We work at fixed hole doping ν = 2 -1/5 and at parameters (g, U 1 ) = (0.15, 1.5), choosing for the moment to set U 2 = 0, though we emphasize that our results are robust to the inclusion of further repulsive interactions (as we illustrate in Sec. V B). As discussed above, we must demonstrate that the ground state exhibits algebraically decaying correlation functions of charge-2e operators. To this end, we work in the computational basis of fermions defined in Eq. ( 27) and introduce the following set of zero y-momentum charge-2e operators:

where x is the center of mass position of the holes, δx is their relative position, 7 and α, β ∈ 0, x, y, z. We then consider correlation functions of such operators:

which we choose to normalize by the average number of available hole pairs N 2h = |2 -ν|/2. This correlation function is shown in Fig. 4(a) for k L y =5 y = 4π/5 and αβ = xx. As the bond dimension increases, the correlation length increases and C xx approaches a power law (the dashed straight line on the log-log scale), signaling the onset of algebraic order. To gain further quantitative confirmation, we perform a scaling collapse of the form

where C (χ ) xx is the correlation function Eq. ( 38) at bond dimension χ and ξ 2e,χ is its correlation length. This scaling collapse, shown in Fig. 4(b), is achieved with a power-law exponent η(L y = 5) = 1.3 ± 0.2. This data confirms the existence of algebraic long-range order in 2e pair correlations at fixed L y .

In addition, we demonstrate that the behavior of the 1e correlations is consistent with a Luther-Emery liquid and not a Luttinger liquid (see Sec. V A 1 for the distinction). To do so, we investigate the ratio ξ 2e,χ /ξ 1e,χ as a function of bond dimension, which we plot in the inset of Fig. 4(b). We find that the ratio increases rapidly with bond dimension which, per the discussion of Sec. V A 1, is inconsistent with the behavior of a Luttinger liquid and suggests a gap in the 1e sector. In fact, employing the recently developed methods of Refs. [116,122], we obtain an extrapolated, and clearly nite, value ξ -1 1e,∞ = 0.11 in units of inverse 4a, the length of the unit cell. In contrast, |ξ -1 2e,∞ | < 0.003 is nearly zero, consistent with our above evidence for algebraic order in this charge sector. This demonstrates that the ground state at L y = 5 is a Luther-Emery liquid, the quasi-1D analog of the superconducting phase.

We conclude with some important phenomenological observations about the superconducting correlations. In particular, algebraic long-range order is most dominantly seen in the following pairing channels:

appearing with nearly equal strength, as well as the s x l 0 channel, which is notably of much weaker strength. As such, the pairing is between opposite spins and predominantly between opposite layers. The approximately equal strength of the s x l x and s y l y correlations are a consequence of an approximate U(1) symmetry generated by s z l z . This symmetry also explains the dominance of interlayer pairing over s x l 0 , which we henceforth ignore. In summary, the dominant condensing pairing channels imply that the superconducting order is due to the condensation of interspin, interlayer hole pairs. We return to this point in Sec. VI C in our analysis of the superconducting pairing symmetry.

Circumference scaling

Having established the presence of a Luther-Emery liquid state at nite circumference L y = 5, we now investigate whether such a state tends to a bona de superconductor in the 2D limit by investigating it at larger circumferences. This investigation is challenging for three distinct reasons. First, as with all DMRG studies, the computational resources required increase exponentially with the 2 × 2 × L y tensors that wrap around the cylinder. To partially compensate for this scaling, we use parallelized iDMRG for circumferences L y = 6, 7 with bond dimensions up to χ = 12 288. As a remark, compared to the Fermi-Hubbard model at the same L y , our projected four-avor model has twice as many degrees of freedom, making it vastly more expensive to simulate. Second, since the electronic lling must be commensurate with the periodicity of the MPS, we are limited to the circumferencedependent densities ν = 2 -1/L y . 8 This makes it dif cult to directly compare the physics at each circumference. Finally, we will later show in Sec. VI that the superconductor pairs holes at an emergent Fermi surface (or more accurately Fermi points at nite circumference). Crucially, the number and location of Fermi points vary signi cantly for small and accessible L y and, unlike in the Hubbard model, are nonmonotonic in cylinder circumference (see Ref. [90]). This makes studying the approach to the 2D limit particularly challenging in our model. Nevertheless, in what follows we analyze the superconducting correlations for each accessible circumference, ensuring that the MPS truncation error is 3 × 10 -5 or better at each bond dimension. Ultimately, we nd signatures of large pair correlations at L y = 6, 7 but are unable to rule out the possibility of a Luttinger liquid at these circumferences.

We rst determine whether the states we nd are gapless in the charge-2e sector by examining correlation functions of the form Eq. (38). These are plotted in Fig. 4(a) at all three circumferences in the αβ = xx channel, with the larger two L y vertically offset for display, and taken at k L y =6 y = -2π/3 and k L y =7 y = -2π/7. Evidently, at each circumference, these correlations become algebraically long-ranged with increasing bond dimension. The exponents of algebraic decay are then extracted by collapsing under a scaling hypothesis analogous to Eq. ( 39) and are shown in Fig. 4(c) for each circumference and for each of the dominant superconducting channels. We nd a quantitative decrease in η(L y ) as L y increases for the xx and yy channels, consistent with the superconducting expectation that η(L y ) → 0 as L y → ∞. For the x0 channels, the superconducting correlations at L y = 6 are found to have a quantitatively different exponent and behavior as compare with the other; we comment in Ref. [90] that this coincides with it having a different pairing symmetry in the x0 channel at L y = 6.

Next, we investigate the behavior of the 1e-correlations, which we nd to be more subtle and ultimately inconclusive. We provide the core results here and give further details in Ref. [90]. Recall that our previous analysis of L y = 5 found that these charge-1e correlations remained exponentially decaying even as χ → ∞. This was evidenced by the growth of ξ 2e /ξ 1e and by performing a systematic extrapolation [116,122] to in nite bond dimension that revealed a nonzero inverse correlation length ξ -1 1e,∞ (L y = 5) ≈ 0.11. The charge-2e and 1e correlations therefore both matched the expectations for a Luther-Emery liquid, Eqs. ( 34) and (35), at L y = 5. This conclusively identi ed the L y = 5 ground state as a quasi-1D superconductor.

At larger circumferences, however, the situation is less clear. At L y = 6, we nd that the behavior of the 1ecorrelations does not conclusively point to either a Luttinger liquid or a Luther-Emery liquid. On the one hand, we nd that the magnitude of ξ 2e is nearly 2.5 times larger than ξ 1e at χ = 12 288, indicating large superconducting pair correlations in the state. This is suggestive of a Luther-Emery liquid. On the other hand, we nd that ξ 2e /ξ 1e is roughly constant with bond dimension, which fails to exclude a Luttinger liquid. Generally, at this circumference, we conclude that the nite-bond dimension state is too far from a 'scaling limit' to properly differentiate the two possibilities. This is evidenced by extrapolating ξ -1 1e/2e,∞ (L y = 6) using the method detailed in Ref. [122] and nding that we are too far from convergence to obtain sensible (e.g., non-negative) values of these quantities.

We conclude by discussing the case of L y = 7. Here, unlike the previous two cases, we nd that ξ 2e is only slightly larger than ξ 1e with a ratio ξ 2e /ξ 1e that remains relatively constant across the largest bond dimensions we can access. As such, at this circumference, both correlation lengths appear to diverge as a function of bond dimension, though much larger bond dimension would be required for a more de nitive determination. Hence, there is little evidence in direct support of a Luther-Emery liquid over a Luttinger liquid. Indeed, the two may exhibit similar behavior when the nite-χ We investigate the extent of the L y = 5 superconducting phase by examining the ratio ξ 2e /ξ 1e (see Sec. V B for more details). In (a), we depict the phase diagram at ν = 1.8 in the absence of interlayer repulsion. A robust LAFpolarized superconducting region appears where the LAF order was found at integer filling, suggesting the LAF may be the parent state for the superconductor. We mark, with a star, the point at which the scaling analysis of Sec. V A was performed and remark that the phase competes with a band metal and LAF-polarized metal, the latter of which retreats with increasing bond dimension (see Sec. V B 1). In (b), we fix g = 0.1 and U 1 = 1.8 and depict our heuristic diagnostic, the growth of the ratio ξ 2e /ξ 1e vs χ , as a function of the electronic density ν ∈ [1.4, 2]. Superconductivity is found across the span 1.7 ν < 2, with metallic LAF behavior beginning at ν = 1.6. In Ref. [90], we show that the system remains LAF-polarized throughout this range, further cementing the LAF as the superconducting parent state. In (c), we examine the robustness of our superconductor to additional repulsive interactions. In particular, fixing (g, U 1 ) = (0.1, 1.4), we chart the phase diagram as a function of interlayer repulsion U 2 and interaction range d [see Eq. ( 41)]. We find the superconductor is remarkably robust, giving way initially to a LAF metallic phase when U 2 is sufficiently large, and then to a layer-polarized metallic phase when U 2 exceeds U 1 . In Ref. [90], we performed a scaling analysis similar to that detailed in Sec V A for the point marked with a red star. correlation length is too small to resolve the charge gap (see Sec. V A 1). Altogether, this highly circumference-dependent behavior prevents us from drawing firm conclusions about the 2D phase.

Superconducting phase: LAF coexistence

We now chart the portions of parameter space that host superconductivity at L y = 5. We begin by studying the (g, U 1 ) plane at a fixed electronic density ν = 1.8 and in the absence of interlayer interactions. Using the diagnostic described above, superconductivity is found to occupy a broad portion of the strong-coupling regime that hosts the LAF insulator at integer filling [see Fig. 5(a)]. Intriguingly, the superconductor itself is always found to have sizable coexistent LAF polarization. Moreover, its correlations are most long-ranged in the channels Eq. ( 40), suggesting a consistent pairing symmetry throughout the superconducting phase. On the other hand, in the weak-coupling regime with small U 1 , we identify a metallic state with nearly vanishing LAF polarization, i.e., a band metal. In between these two phases, we observe a metallic state with significant LAF polarization.

Similar to our identification of the LAF phase boundary at ν = 2, we demarcate the transition between the weakcoupling and LAF metallic phases by a dotted boundary corresponding to the largest slope ∂n z LAF /∂U 1 . In contrast, the boundary between the LAF metallic and superconducting phases is more sensitive to bond dimension, even at the sizable value χ = 5793. This exemplifies a pattern that is pervasive in our numerical data: LAF metallic states commonly undergo a finite-χ transition into the LAF superconducting phase at a bond dimension that depends sensitively on parameters, whereas the reverse is exceedingly rare. As a result, we expect that the set of points that meet the criteria for superconductivity will continue to grow with bond dimension, and that the LAF metallic phase may purely be an artifact of the finite accuracy of our simulations.

The intrinsic coexistence between superconductivity and LAF polarization is further elucidated by tuning the electronic density in the range ν ∈ [1. 4,2]. To access density increments of 1/10 of a band, we quadruple the MPS unit cell. In Fig. 5(b), we plot our diagnostic for the superconducting order-the ratio ξ 2e /ξ 1e as a function of bond dimension-and nd it that indicates the persistence of superconductivity in the range 1.7 ν < 2. Moreover, in Ref. [90], we demonstrate that the system remains substantially LAF-polarized across this parameter regime, with n z LAF decreasing only linearly in the hole doping. At the remaining parameter points 1.4 ν 1.6, we nd evidence of metallic correlations and possible translation symmetry-breaking along the x direction, though we caution that a careful circumference scaling analysis would be required to identify the precise nature of this nearby phase. Across the entire range of llings, and in particular across the superconducting region, we nd a persistent LAF polarization that nearly saturates the maximal value ν, i.e., the total density of electrons. This ubiquitous coexistence with LAF order leads us to hypothesize that the LAF insulator is the "parent state" of the superconducting state. We further corroborate this claim in Sec. VI where we identify the superconductor as a weak-pairing instability of holes at the LAF's mean-eld Fermi surface.

Stability to further repulsive perturbations

We now show that the superconductor is robust to both large interlayer repulsion and a nite interaction range. To this end, we replace the microscopic on-site interactions with a longer-range Gaussian, parameterized as

where d speci es the interaction range. In Fig. 5(c), we show the extent of the superconducting phase in the (d, U 2 ) plane at the xed parameters (g, U 1 ) = (0.1, 1.4). Remarkably, in spite of the pairing being predominantly interlayer (see Sec. V A), we nd the superconductor survives in the presence of substantial interlayer repulsion U 2 ∼ O(U 1 ) and is further resilient to simultaneously increasing the range of the interactions. While the superconducting phase is globally identi ed by using the heuristic diagnostic introduced at the beginning of Sec. V B, we cement its existence at L y = 5 in the presence of additional interlayer repulsion by performing a detailed scaling analysis of the superconducting correlations, analogous to the one performed in Sec. V A. We also see indications, via a slow but monotic increase in the ratio ξ 2e /ξ 1e , that this behavior may persist at L y = 6. These analyses are detailed fully in Ref. [90] and are performed at the parameter point (d, U 2 ) = (0.358, 0.2) [marked with a star on Fig. 5(c)]. Despite the resilience of the superconducting phase, we nd that when d is increased suf ciently, the superconductor gives way to a translation-invariant, LAF-polarized metallic state. Moreover, we remark that increasing U 2 in excess of U 1 induces a transition out of the LAF phase into a layerpolarized, translation-invariant metal that originates from the ν = 2 QAH insulator.

A. Numerically probing the superconducting pair wavefunction

The de nition of the pair wave function (42) implicitly assumes that the ground state wave function is a superconducting condensate in which the number of Cooper pairs in the ground state wave function is uncertain. This presents an obstacle to extracting this quantity in our iDMRG calculations, as U(1) particle number conservation causes (k) to vanish uniformly. Here, we provide a brief account of how to circumvent this dif culty, leaving a more complete and precise explanation to Ref. [90]. Broadly, we do so by introducing a proxy pair wave function αβ (k) that can be extracted within iDMRG and is proportional to αβ (k). Readers more interested in our physical results may skip to the next subsection where this quantity is used to probe the structure and pairing symmetry of the superconductor.

In the L y → ∞ limit, the de nite U(1) Q charge of the iDMRG ground state implies that it will fail to satisfy cluster decomposition, de ned as [123] lim

As a result, long-range superconducting order [de ned in Eq. ( 33)], does not imply a nite expectation value of the Cooper pair operator. Nevertheless, let us assume that there exist cluster decomposition-satisfying ground states |ϕ⟩ of the 2D superconducting manifold for which

The pair wave function can then be extracted numerically as follows. First, we de ne the proxy pair wave function by the following correlation function:

Here ˆ αβ (x; 0, p y ) is as de ned in Eq. ( 37), while the ˆ † term is de ned by a Fourier transform in the relative position coordinate:

which, when averaged over the center of mass x, is equal to the Cooper pair operator de ned in Eq. ( 42).

The proxy pair wave function-which need not vanish in U(1)-conserving iDMRG since it involves the expectation value of a charge-neutral operator-should be thought of as a function of k (with xed parameters x and p y ). Provided that we consider translationally invariant states and take x to be larger than both the characteristic size of the Cooper pair (the "coherence length") ξ SC and the connected correlation length ξ , then for suf ciently large L y we expect that:

where the prefactor is a k-independent constant and |ϕ⟩ is a cluster decomposition-satisfying superconducting ground state of the 2D system. As such, by examining the proxy pair wave function (k) as a function of momentum, we gain direct insight into the structure of the true pair wave function (k).

B. Emergent weak pairing at the LAF Fermi surface

In this section, we demonstrate that our superconductor arises from the momentum space "BCS" pairing of holes above the LAF insulator. To be explicit, we will provide evidence that this emergent weak-pairing picture occurs as a two-step process. First, strong interactions select LAFpolarized ground states, even in the presence of hole doping. Such "normal" states are well-described by a LAF bandstructure that strongly renormalizes the four noninteracting bands into two valence bands and two conduction bands. Second, holes-initially occupying a normal LAF metallic state-pair at the renormalized Fermi surface, triggering a superconducting instability. Before continuing, we emphasize once again that this weak-pairing picture emerges in a setting where repulsive interactions are the dominant energy scale.

As preliminary evidence for this picture, we plot the density of electrons in the superconducting state as a function of momentum [see Fig. 6(a)], nding that the density is depleted in localized "pockets" of momentum space but is otherwise uniform. To understand the origin of these disconnected Fermi surfaces, we examine the mean-eld band structure of the LAF insulator computed using standard (see e.g. [80]) selfconsistent Hartree-Fock (SCHF). We show these bands in Fig. 6(b), with the lower two degenerate bands corresponding to the occupied electronic states of the LAF insulator and the upper two degenerate bands corresponding to the spectrum of gapped single-electron excitations above the LAF. At lowest order in g, these two sets of bands are merely the bare Hofstadter bands, but massively split by the interaction scale. In Figs. 6(a) and 6(b), we mark the Hartree-Fock Fermi energy at lling ν = 1.8 with orange lines, nding that their locations coincide closely with drops in the iDMRG electronic density.

To further corroborate the BCS picture described above, we examine the structure of the pair wave function, Eq. ( 42). We rst lay out our expectations. If BCS theory was applicable, we would expect the pair wave function to be related to the Bogoliubov gap δ αβ (k) by

where E HF (k) is the dispersion of holes above the normal LAF state-corresponding to the Koopman spectrum of hole excitations in self-consistent Hartree-Fock-and E F is the Fermi level. Generally, we expect the superconducting gap δ αβ to be smaller than the characteristic scale of the dispersion, namely the bandwidth. Within the validity of BCS mean eld theory, we would therefore expect Eq. ( 48) to vanish suf ciently far from the Fermi surface. On the other hand, the pair wave function would peak at the Fermi surface, at which it would be proportional to the phase of the gap function:

These expectations are borne out in our iDMRG numerics. In particular, in Fig. 6(c,d) we plot the proxy pair wave function (45) and nd that it peaks precisely at the Fermi surface and decays quickly away from it. In Fig. 6(c), we zoom into one Fermi pocket and observe excellent agreement between the Fermi surface predicted from Hartree-Fock, the location at which the charge density drops off in iDMRG, and the peak in the proxy pair wave function. In Fig. 6(d), we plot the latter over the Brillouin zone and nd that its magnitude is maximal at each of the four disconnected Hartree-Fock Fermi surfaces, providing further evidence that the superconductor arises as an instability of the LAF Fermi surface. We leave the analysis The electronic density of the putative superconducting iDMRG groundstate for the L y = 6 cylinder, shown here at χ = 12 288 and (g, U 1 , ν) = (0.15, 1.5, 2 -1/6), is uniform in the majority of the Brillouin zone but sharply depleted in isolated regions. The locations of these drops in the density match the mean-field Fermi surface of noninteracting holes above the LAF insulating state, marked in orange (in all subplots) and obtained using self-consistent Hartree-Fock (SCHF). In (b), we plot the SCHF band structure computed at ν = 2 on 48 × 48 unit cells. The energy E HF of the twofold degenerate valence band is shaded in blue and defined relative to the valence band edge. For small hole-doping away from ν = 2, the Fermi surface (shown in orange) consists of four disconnected pockets. As SCHF enforces electronic U(1) (disallowing pairing), this strongly interacting effective bandstructure is a metallic "parent state" for superconductivity. In (c), we zoom into one particular pocket-specified by blue and magenta boxes in (a) and (d), respectively -and plot both the electronic density (left axis) and the proxy pair wave function Eq. ( 45) (right axis) obtained from iDMRG. The former drops off rapidly, while the latter sharply peaks, precisely at the location of the SCHF Fermi surface, strongly suggesting that the superconductor arises from a weak-pairing instability upon hole-doping the ν = 2 insulating LAF state. (d) gives a global view of the proxy pair wave function, whose real part is orders of magnitude larger than its imaginary part. From its sign structure, we deduce the unconventional pairing symmetry of the superconductor, namely that it is p-wave and is phenomenologically consistent with a Bogoliubov gap function δ(k) ∼ sin(k y )(s + l + + s -l -). and interpretation of the relative complex phase of αβ (k) to the next subsection.

C. Pairing symmetry

In what follows, we characterize the pairing symmetry of the superconducting state by using the proxy pair wave function (k). To be precise, the pairing symmetry corresponds to the irreducible corepresentation-the generalization of a representation to incorporate anti-unitary symmetries-of the magnetic space group under which the pair wave function transforms [1,4,90,124]. In our case, since the Cooper pairs that condense to form the superconductor carry zero total momentum, we need only specify the transformation properties under the discrete rotations C 2z and C 2x , the nontrivial magnetic translation m x , and the anti-unitary spacetime inversion symmetry I (see Ref. [90]). Since these operators each square to the identity and mutually commute in the zero-momentum sector, these generators act on Cooper pairs as the group (Z 2 ) 3 × Z K 2 with anti-unitary Z K 2 . Given that the co-representations of Z K 2 are trivial [90], specifying the irreducible representation amounts to identifying three Z 2 characters for each condensed superconducting channel s α l β .

Let ĝ be the unitary operator corresponding to a symmetry g ∈ {C 2x , C 2z , m x }. We define the corresponding characters ν αβ g by [1,4,90]

where ν αβ g = ±1 and αβ (no implicit summation) correspond specifically to the dominant condensed channels specified in Eq. (40). Intuitively, the characters capture how the pair wave function transforms under the application of symmetries on the ground state. More generally, the role of these characters is played by matrices ν g that form a co-representation of the magnetic space group (see Ref. [90] for a more formal and general treatment [90]).

As an illustrative example, we compute the index for g = C 2z . Its action on our Cooper pairs is given by

In the case αβ = xx, whose proxy pair wave function we have plotted in Fig. 6(d), we readily observe that xx (-k) = xx (k). From this and Eq. ( 49), we immediately find that:

We also find that that ν yy g = ν x0 g = -1, indicating that all condensed channels are p-wave in nature.

A similar analysis can be performed for the other symmetries. In particular, inspection of the proxy pair wave function likewise leads to the conclusion

for both dominant pairing channels αβ = xx, yy. Finally, for the anti-unitary symmetry I, our choice of layer-polarized basis guarantees that Î ˆ † αβ (k) Î-1 = ˆ † αβ (k). This implies that the pair wave function has a constant phase, assuming that I remains unbroken in the ground state. Indeed, our numerically obtained proxy pair wave function is found to be purely real.

We conclude by remarking that our numerical observations are consistent with a superconductor whose Bogoliubov gap function δ αβ (k) has the same symmetry properties as sin(k y ) for αβ = xx, yy. Moreover, in Ref. [90], we provide numerical evidence that the pair wave function in the αβ = xx channel appears with a relative minus sign compared to the αβ = yy channel. As such, the symmetry properties of a putative gap function can be summarized as

Consistent with a parent LAF state, the gap function corresponds to nonunitary pairing with δ † (k)δ(k) ∝ (1 + s z l z )/2 = P LAF , the LAF projector. Furthermore, while this form of the gap function typically implies the presence of symmetry-enforced nodes along the lines k y = 0, ±π , our superconductor instead remains gapped. This is because the Fermi surface at which pairing occurs is disconnected and does not intersect either of these lines in momentum space (see Fig. 6).

D. Estimate of transition temperature

The weak-pairing picture for the superconducting state, advocated for above, provides a heuristic framework for estimating the superconducting transition temperature. To this end, we recall that the nite temperature transition out of a (2 + 1)d superconductor occurs either due to (i) the loss of phase coherence and algebraically long-ranged phase correlations or (ii) the unbinding of Cooper pairs. The former is described by a Berezinskii-Kosterlitz-Thouless (BKT) transition corresponding to the proliferation of vortices in the superconductor and has a transition temperature that is typically upper bounded by the ground state super uid stiffness ρ s as [1,4]

In contrast, the Cooper pair unbinding transition occurs at a temperature scale set by the maximum Bogoliubov gap δ SC at the Fermi surface [1,4,90]:

With input from self-consistent Hartree-Fock, we can estimate both scales. In a BCS-like superconductor, the super uid stiffness at T = 0 is set by the density of paired dopants divided by their effective mass, which in two dimensions scales as the Fermi energy [125]. As such,

On the other hand, the approximate magnitude of the Bogoliubov gap δ SC can be extracted from the pair wave function αβ (k). In particular, for momenta around the Fermi surface k ≈ k F , Eq. ( 48) implies that

Here, we made the linear approximation

, with v F being the Fermi velocity, and took δ αβ (k) ≈ δ αβ (k F ). As such, around the Fermi surface, 4 2 αβ (k) is predicted to have a Lorentzian pro le whose half-width-half-maximum is equal to δ αβ (k F )/v F . Maximizing over the dominant channels αβ = xx, yy and taking v F = 0.17 obtained from self-consistent Hartree-Fock at L y = 5, we obtain the following estimate for the Bogoliubov gap:

In particular, we nd that the pair-breaking temperature T p c is comparable to or less than the BKT temperature. Though we expect the superconducting transition temperature to therefore be set by the energy scale δ SC , we carefully note the possibility of a BKT-type thermal transition out of the superconducting phase. This is because of a reduction in the condensate fraction (equivalently, the Cooper pair density) at nite temperature due to thermal uctuations, which can signi cantly decrease the super uid stiffness ρ s as we approach the pairbreaking scale. Consequently, as T → T p c , the small value of ρ s can decrease the cost of vortices and cause them to proliferate, triggering a BKT transition that preempts T p c .

Optical engineering of double Hofstadter hopping

To engineer the in-layer double Hofstadter hopping, consider placing atoms in a two-dimensional "magic wavelength" FIG. 7. Optical Lattice Implementation. The double Hofstadter model is engineered by placing fermionic alkaline-earth atoms in a two-dimensional optical lattice. Layer and spin are encoded in the clock and nuclear spin states of the atom, respectively (see boxed levels). Furthermore, complex hoppings are generated by first creating a period-two staggered on-site potential (red and green sites; potential sketched on the right) to suppress bare hopping in the y direction, labeled by n in Eqs. ( 60)- (62). This period-two potential crucially has opposite signs for the opposite layers. Second, resonant hopping is restored separately on the even and odd bonds using two pairs of lasers (blue and orange arrows) that form a running wave in the x direction and are retroreflected (shown in lighter colors) to form a standing wave in the y direction. These restored hoppings are then imbued with the spatially dependent phases of the laser drives. In this scheme, atoms see an artificial magnetic flux ℓ = ℓπ /2 that is opposite in opposite layers, as desired. optical lattice with lattice spacing a [130,131], created by superimposing two orthogonal standing waves. To realize the requisite complex hopping, our strategy is to energetically suppress bare nearest-neighbor tunneling in the y direction, and then restore it using Floquet driving [132,133], thereby imprinting the spatially dependent phase of the drive onto the restored hoppings. This approach combines two experimental protocols implemented by one of us [84,134], where the key ingredients are (1) a homogeneous artificial magnetic field and (2) a layer-dependent direction of the field, which is realized using a state-dependent potential at the "anti-magic wavelength" of the two clock states.

To suppress hopping, we superimpose on the optical lattice an additional standing wave in the y direction with lattice spacing 2a at the antimagic wavelength, inducing opposite on-site potentials for the two clock states (see Fig. 7). Neighboring lattice sites will therefore have a site-alternating energy mismatch of ± , which is always opposite between atoms in the |e⟩ and |g⟩ states (i.e., e = ± =g ). This energetically suppresses nearest-neighbor hopping in the y direction. To restore it, we use the two-laser optical driving scheme introduced in Ref. [134], which resonantly modulates the on-site potential on even and odd bonds in the y direction independently. Intuitively, since the different clock states see opposite staggered potentials, one will absorb energy from the drive while the other will emit energy into it. As a result, their hoppings acquire conjugate complex phases. This intuition is borne out at zeroth-order in a Floquet-Magnus expansion [90], wherein this driving protocol gives an effective average Here ψ † m,n is the fermionic creation operator, which is a row vector in the spin and layer [e.g., see Eq. ( 1)], at lattice site m, n ∈ Z 2 . The layer-dependent phases are

where ℓ = +1 (-1) for the top (bottom) layer. 9 These phases produce the same flux pattern as in Eq. ( 2), namely, a uniform π/2 (-π/2) flux through each plaquette in the top (bottom) layer (see Fig. 1).

Hubbard interactions

While the single-particle terms did not depend on the particular alkaline-earth atom used, the specific values of the Hubbard-U couplings depend sensitively on this choice. The coupling between two alkaline-earth atoms at the same lattice site is proportional to one of four independent scattering lengths a gg , a ee , a + eg , and a - eg , depending on the internal states of the atoms [135]. The first two of these specify the interaction when the two atoms are both in the |g⟩ and |e⟩ clock states, respectively, while the latter two correspond to symmetric or anti-symmetric occupation of the |g⟩ and |e⟩ states. Crucially, these scattering lengths do not depend on the nuclear spin due to an approximate SU(2I + 1) symmetry [67,135,136], making them independent of the nuclear spin states chosen for the encoding in Eq. (59). Translating the atomic labels to the language of the double Hofstadter model, the realized interacting Hamiltonian takes the form

where U T/B 1 ∝ a gg/ee controls the strength of two independent in-layer interaction strengths, U 2 ∝ (a + eg + a - eg ) controls the interlayer interaction strength, and J ex ∝ (a + eg -a - eg ) is an additional SU(2)-symmetric "Heisenberg" spin-exchange interaction between the layers which conventionally is ferromagnetic (antiferromagnetic) for J ex > 0 (J ex < 0).

To obtain the LAF phase at integer lling and the superconductor upon doping, the numerical results of Secs. IV and V show that it is favorable to be in an interaction regime where in-layer interactions U T/B 1 are large while interlayer interactions U 2 are small. Furthermore, if the additional J ex term is ferromagnetic, it will disfavor LAF states in favor of spinpolarized states and hence will need to be smaller than the perturbatively generated "superexchange" term that selects out the LAF [cf. Eq. ( 25)]. These requirements on the interactions place constraints on the possible atomic species that can be used to realize our proposal, as their scattering lengths must realize a favorable ratio between U T/B 1 , U 2 , and J ex . However, we now show that simply using an additional clock-state dependent optical potential can systematically suppress J ex , U 2 (see Ref. [90] for more details), thereby ameliorating these constraints.

To see this, recall that the lattice interaction parameters in Eq. ( 63) arise from the overlap between the Wannier functions of the optical lattice potential. In particular, J ex , U 2 ∝ r∈R 3 |w e (r)| 2 |w g (r)| 2 where w e/g (r) are the Wannier functions of the atoms in either clock state. By using a clock-state-dependent potential to displace atoms in different clock states by a distance z δ , such overlaps are suppressed exponentially in z δ , enabling one to effectively eliminate these two couplings. 10 Indeed, under an approximation that the Wannier functions are the wave functions of the lowest harmonic oscillator (which we justify in Ref. [90]), these overlaps fall off as ∼e -z 2 δ /ξ 2 /ξ , where ξ ∼ 1/V 1/4 z is a localization length. Crucially, this leaves the strength of the in-layer couplings unchanged as they are set by on-site Wannier overlaps for each clock state independently.

B. Relationship to moiré materials

Recall that we formulated the double Hofstadter modeland then developed its direct experimental implementation with cold atoms-in an effort to explore, in a minimal setting, the interplay between a set of ingredients commonly found in moiré materials (see Sec. I). Nevertheless, the presence of both correlated insulating states and unconventional superconductivity in this model, both appearing independently in experimental and theoretical studies of moiré graphene platforms, invites us to make precise the connection to such systems. In this section, we reveal a close structural correspondence with the chiral model of magic-angle twisted bilayer graphene 11 [78,80,82,140].

At the level of their band structures, the similarities between the two models are already striking. In chiral TBG, the low-energy single-particle subspace also consists of energetically narrow bands with w 2 = 1 fragile topology in each avor sector (i.e., valley and spin) and, consequently, hosts two Dirac cones of the same chirality per avor [52,53,62,101,[141][142][143]. Furthermore, just as the low-energy bands of the double Hofstadter model are naturally expressed in a Chern basis (11) of layer-polarized orbitals, the at bands of chiral TBG are conveniently described by a basis of graphene sublattice-polarized Chern bands [80].

Going beyond single-particle physics, we will construct an explicit dictionary between the interacting models. To do so, let us rst recall some details of the strong-coupling theory of TBG. The narrow bands of TBG are modeled by electrons carrying electronic spin s and valley pseudospin τ , with each electron occupying one of two narrow sublatticepolarized bands labeled by σ , for a total of eight avors. At the magic angle, these sublattice-polarized bands are exactly at and the chiral model has a large U(4) × U(4) symmetry generated by independent spin and valley rotations in each Chern sector [80,82]. In Fig. 8 (top), we depict the division of the eight at bands into two U(4)-symmetric subspaces with opposite Chern numbers. Since experiments on TBG suggest some degree of avor polarization in the vicinity of the superconductor (see, e.g., Refs. [144][145][146]), many theoretical works reduce the eight degrees of freedom of the U(4) × U(4) model to four active flavors. The result is a model with a U(2) × U(2) continuous symmetry (corresponding to charge conservation and independent rotations of a pseudospin) within each Chern sector γ z = σ z τ z [80,147]. These flavor-polarized models can be mapped explicitly to the double Hofstadter model.

Here, we consider two particular routes to flavor polarization studied in the literature [see Fig. 8 (middle)]. In the first route, the electronic spin is assumed to be fully polarized, and the resulting pseudospin η is related to the valley. Each U(2) symmetry is then generated by the electric charge, along with a pseudospin [139] η = (σ x τ x , σ x τ y , τ z ).

In the second route, recently posited in Ref. [147], electronphonon coupling causes the electrons to polarize into particular intervalley-coherent orbitals [see Fig. 8 (middle, right)]. The resulting pseudospin is then just the electronic spin

whose components, along with electric charge, similarly serve as generators for each U(2) symmetry.

To map these two flavor-polarized models to the double Hofstadter model, we recall from Sec. III D 2 that the latter has an enhanced U(2) × U(2) continuous symmetry in the decoupled (g = 0) limit. We may then identify the Chern index and pseudospin of the chiral models with the layer and electronic spin of the decoupled double Hofstadter model, respectively:

Maps 1 and 2, summarized in Fig. 8 (bottom), are bijections of both the degrees of freedom and continuous symmetries between the double Hofstadter model and flavor-polarized models of chiral TBG. This identification extends to the Hamiltonian level. Since the dispersion vanishes at the magic angle in the chiral limit, the projected Hamiltonian ĤTBG = q V q : ρq ρ-q : with ρq = k ĉ † k TBG q (k)ĉ k+q is entirely specified by the form factors TBG [c.f. Eq. ( 16)]. Symmetries constrain their matrix structure to be TBG q (k) = F q (k)e i q (k)γ z [139]. For the double Hofstader model, the form factors (19) in the decoupled limit are q (k) = F S q (k)e i S q (k)l z , which therefore have exactly the same structure. Moreover, departure from the magic angle in chiral TBG introduces a nonzero dispersion hγ x/y that hybridizes the two Chern sectors (see Fig. 8) [56,80,139]. Such a term maps directly onto the interlayer tunneling gl x/y of the double Hofstadter model. 12 Maps 1 and 2 therefore give an almost term-by-term correspondence in matrix structure between the Hamiltonians.

The similarities between the Hamiltonians suggest that their ground states may also be analogous. To check this, we apply the above mappings to the LAF ground state at ν = 2. Recall that the LAF state occupies the layer-polarized bands in the pattern: (68) Under map 1, the LAF becomes either the Valley Hall (VH) state or the Kramer's Intervalley Coherent State (KIVC) (which are symmetry-equivalent in chiral TBG). Strikingly, these are the ground states of the strong coupling model of chiral TBG at |ν| = 2 [80,82]. Under map 2, the LAF instead corresponds to the "TIVC-QSH" state, which has intervalley coherence and forms a spontaneous quantum spin Hall insulator in the microscopic electronic spin. Studies of TBG where electron-phonon terms favor the relevant flavor polarization of route 2 indeed find a TIVC-QSH ground state [147]. The LAF therefore maps to ground states of flavor polarized models of TBG. This is remarkable. Indeed, this suggests that the double Hofstadter model-which minimally incorporates the ingredients introduced in Sec. II-may yield valuable conceptual insight into the interacting phenomena present in chiral TBG.

We conclude this section by highlighting notable differences between chiral TBG and the double Hofstadter model. Manifestly, the unit cells are quite different: hexagonal with C 6z symmetry for unstrained TBG, versus rectangular in the double Hofstadter model. Chiral TBG also has an exact particle-hole symmetry [139], whereas particle-hole is only approximately satis ed in the lowest-band-projected double Hofstadter model, becoming exact only in the q → ∞ limit. Another obvious difference lies in the interactions, which are local in the double Hofstadter model but longer ranged (screened Coulomb) in models of TBG. While the range of interactions is expected to affect the energetics of excitations above the aforementioned insulating ground state [148], the relevance of the different discrete symmetries is less obvious and requires further investigation.

However, the key difference between the two models lies in the gauge-invariant properties of the Bloch wave functions [149]. These include the Berry curvature, or more generally the "quantum geometry" (reviewed in Refs. [139,150]), and determine numerous physical quantities like the interactiongenerated (Hartree-Fock) dispersion [80]. A key feature of TBG is that all these quantities are sharply concentrated near the point, which is known to strongly affect the ground state phenomenology (see, e.g., Refs. [151][152][153][154]). In the double Hofstader model, by contrast, these quantities are distributed almost uniformly across the Brillouin zone, and the quantum geometry is much closer to an ideal "vortexable" limit [155] (comparable to TBG at a chiral ratio w 0 /w 1 ∼ 0.3; details Ref. [90]). Incorporating such an inhomogeneity into the double Hofstadter model could help gain a "bottom-up" view of the phenomenology of TBG and is an interesting direction for future research. Nevertheless, given that experiments in magic-angle TBG [156] and its multi-layer generalizations [22] point towards nodal pairing [157], the emergence of pwave unconventional superconductivity in the related double Hofstadter model makes it a valuable departure point to explore pairing in more realistic moiré models.

VIII. DISCUSSION

In this paper, we introduced the double Hofstadter model, a simple lattice model featuring strong repulsive interactions, narrow topological bands, and time-reversal symmetry. By using state-of-the-art cylinder iDMRG, we investigated the ground state phase diagram of this model. We found strong numerical evidence for the existence of a robust p-wave superconductor at circumference L y = 5 upon hole-doping a quantum spin Hall (QSH) insulator appearing at half-lling of the at band subspace. This correlated insulator, which we termed the layer antiferromagnet (LAF), has opposite spin polarization in opposite, time-reversal-related Hofstadter layers. Using a novel technique for probing the gap function of the superconducting condensate, we deduced both its pairing symmetry and characterized it as a "BCS"-like superconductor comprised of Cooper pairs near the mean-eld LAF Fermi surface. We also presented an experimental blueprint for implementing the double Hofstadter model in near-term optical lattice setups and elucidated a close correspondence to chiral models of twisted bilayer graphene (TBG).

Before discussing our work in a broader context, we comment on the extent of our evidence for superconductivity and routes toward a more de nitive veri cation. Recall that our results were enabled by exploiting recent developments in the compression of matrix product operators [66], extensive parallelization, and working at the largest accessible bond dimensions. Although this enabled us to nd de nitive superconducting signatures at L y = 5, the story at the larger circumferences, though promising, remains unclear. Indeed, this is a consequence of a fundamental limitation of the cylinder iDMRG method, namely the exponential complexity of scaling the circumference. An important and interesting direction for future work is to more conclusively demonstrate superconductivity in the 2D limit, either by developing new numerical methods to enable larger scale exploration within cylinder iDMRG, or by developing and utilizing intrinsically 2D numerical methods [158].

Caveats aside, one of the most intriguing features of this superconducting state is the persistence of superconductivity down to a very small LAF Zeeman pinning eld, Eq. ( 28). While we introduced this as a numerical conveniencenamely for stabilizing the LAF phase-its numerical value in this study is smaller than all other microscopic energy scales. It is therefore likely that superconductivity is stabilized in the absence of any Zeeman eld. In this limit, the LAF is a spontaneous symmetry-breaking state whose various lowenergy uctuations may play an essential role in the pairing mechanism.

This phenomenological feature of the model suggests intriguing connections to previous literature on superconductivity in doped, spontaneous QSH systems [159,160] and important qualitative difference from prior numerical studies of superconductivity in models where the spin rotational symmetry is broken explicitly by strong spin-orbit coupling, e.g., interacting Kane-Mele systems [96,[161][162][163]. In the former case, superconductivity arises from the condensation of lowlying charge-2e skyrmions of the QSH order parameter. A recent microscopic treatment of this mechanism demonstrated pairing within a purely repulsive model, which was proposed as a theory for superconductivity in TBG [56,57,164]. The skyrmion mechanism was studied numerically both in coupled opposite-eld Landau levels using iDMRG [57] and in models with effective local interactions amenable to QMC simulation [165,166], as well as within effective eld theory approaches [58].

An interesting future direction would be to study the relationship between a skyrmion superconductor and that which we have presented here. In particular, while we have argued that our superconductor is "BCS"-like and is phenomenologically consistent with the weak pairing of holes at a Fermi surface, a skyrmion superconductor is, in its simplest manifestation, formed from tightly bound, bosonic skyrmion "molecules." Interpolating between these limits, for instance through spin polarons or baby skyrmions, is an active area of research [167,168]. In the case where both superconductors have p-wave pairing symmetry, such a transition could constitute a novel realization of the p-wave BCS-BEC phase transition [169,170].

Given that topology plays an essential role in a skyrmion superconductor by imbuing individual skyrmion textures with electric charge, it is important to understand the extent to which the weak-pairing superconductor presented here depends on the topology of its parent state. For instance, one can imagine modifying the single-particle problem so that the low-energy bands are exponentially Wannierizable but the spread of these Wannier functions is suf ciently broad to favor a spontaneous-but topologically trivial-LAF insulating state [97]. One could determine whether superconductivity still manifests in this context and whether it is similarly robust.

Finally, we further remark on routes to experimentally realizing the double Hofstadter model, both in solid-state and cold-atom settings. In the main text, we highlighted a precise mapping to the degrees of freedom of chiral TBG, but also commented on how these systems may differ energetically. Given these differences, an interesting direction of future research would be to identify solid-state materials that realize double Hofstadter systems more faithfully, for instance in the family of moiré TMDs. In the cold atom context, we presented a direct implementation of the double Hofstadter model using alkaline-earth atoms in an optical lattice. While our study primarily sought to understand zerotemperature behaviors, directly accessing ground state physics in cold atom quantum simulators remains an outstanding challenge. Namely, experiments typically either remove entropy by coupling to nite-temperature thermodynamic reservoirs [38,171] or indirectly probe the ground state through quasiadiabatic state preparation [41]. These obstacles motivate a systematic exploration of both the nite-temperature equilibrium phenomenology of the double Hofstadter model (e.g., applying the insights developed in Ref. [172]), as well as nonequilibrium dynamical state preparation in this system, which has the potential to drastically in uence the observed order [173].