Repulsion-Driven Topological Superconductivity in Twisted MoTe$_2$

The potential emergence of superconductivity from purely repulsive interactions has been extensively studied in non-topological settings, most prominently the intermediate-strength Hubbard model [1]. Surprisingly, numerically demonstrating that systems such as the (extended) Hubbard model exhibit all-electronic superconductivity remains a challenge. Two dimensional materials offer a distinct route to all-electronic superconductivity, as their flat bands feature non-trivial band geometry, which, together with orbital ferromagnetism, can result in Chern bands. In twisted MoTe$_2$ [2, 3] and rhombohedral graphene [4], fractional fillings of these Chern bands gives rise to fractional quantum Hall effects at zero magnetic field. Can the interplay of repulsive interactions, band geometry, and fractionalization give rise to all-electronic superconductivity? Recent experiments in twisted bilayer MoTe$_2$ (tMoTe$_2$) [5] report superconductivity near an $\nu=-2/3$ fractional Chern insulating state, bringing this question into a concrete microscopic setting.

The experimental phase diagram of twisted bilayer MoTe$_2$ (Figure 1a) in the doping $\nu$ and displacement field $D$ plane is replete with strongly correlated and topological phases. Doping holes into the two narrow (“flat”) Chern bands near $\theta=3.89^\circ$ (Figure 1b) produces a ferromagnetic Chern insulator at $\nu=-1$ with $\sigma_{xy}=-e^2/h$ [2, 3, 6, 7], part of a broader ferromagnetic regime at $\nu\sim-1.3$. Fractional Chern insulators (FCIs) [4, 8, 9, 10, 11, 12, 13, 14] are observed at $\nu=-2/3$ and other fractions, the first realization of the phase at zero magnetic field [2, 3, 6, 7]. In certain samples, doping additional holes on top of the FCI leads to, in sequence: a re-entrant integer quantum anomalous Hall (RIQAH) phase, a resistance peak at $\nu\sim-0.7$ and — finally and most strikingly — a superconducting dome centered at $\nu\sim-0.73$ [5]. (Thus far, superconductivity is only observed in a few samples, all near 3.9°.) The superconductor is unconventional, as ferromagnetism prohibits $s$-wave pairing, while the nearby FCI suggests strong interactions are essential.

The top valence band in MoTe$_2$ [15, 16, 17] is a Chern band close to the “ideal” limit [18, 19, 20] yet—unlike the lowest Landau level—its finite bandwidth allows for mobile charge carriers. Theorists have proposed a range of superconducting mechanisms in this setting, ranging from weak-coupling Kohn–Luttinger pairing enhanced by Chern band form factors and quantum geometry [21, 22, 23, 24, 25, 26, 27, 28, 29], to strong-coupling anyon superconductivity motivated by the nearby FCI [30, 31, 32, 33, 34, 36, 37, 38, 39, 40, 41, 42] [35]. Recent numerics have examined idealized models motivated by tMoTe$_2$: Ref. [36] used the density matrix renormalization group (DMRG) to study a model with a periodic potential but uniform quantum geometry, while Ref. [28] used exact diagonalization to study a model with non-uniform quantum geometry but no periodic potential. Chiral $f$-wave superconductivity was found in both cases. However, microscopic models of tMoTe$_2$ are considerably more complex, due to their periodic potential and their non-uniform quantum geometry. Whether repulsive interactions produce superconductivity in such microscopic models of tMoTe$_2$ is unknown. Here we study a microscopic model of tMoTe$_2$ with gate-screened Coulomb interactions using DMRG on an infinite cylinder. Our main findings are fivefold: (I). We find superconducting correlations at all accessible cylinder circumferences, consistent with a robust superconducting phase in the 2d limit. (II). The superconductor has chiral $f-if$ pairing and an odd number of chiral Majorana edge modes, making it a topological superconductor (TSC). (III). Its single-electron gap is $\Delta_{1e}<0.5$ meV, with zero-temperature superfluid stiffness $\rho_s\sim30$ µeV, giving an upper bound $T_c\lesssim500$ mK. (IV). The phase diagram contains: a superconducting dome at $-0.78\lesssim\nu\lesssim-0.72$ at zero displacement field; a metal; competition between FCIs and a topological CDW at $\nu=-2/3$; and an intervening charge ordered (CO) phase. (V). Superconductivity is favored by shorter-ranged interactions but persists up to effective screening lengths of $d_s\sim5$ nm.

Model & Methods

We study tMoTe$_2$ at twist angle $\theta=3.89^\circ$ and hole dopings $-0.8\le\nu\le-0.63$. Experiments find strong valley polarization throughout this doping range [2, 3, 5, 6, 7, 43], so we retain only the $K$ valley (spin is locked to valley). The single-particle Hamiltonian follows a standard continuum model for twisted transition-metal dichalcogenide homobilayers [15] that acts in the layer basis $\ell=(t,b)$ by

Here $h_\ell$ contains the monolayer valence-band dispersion and moiré-periodic intralayer potential, and $T(\mathbf{r})$ describes the moiré interlayer tunneling. App. A gives their precise forms and all other model details. The opposite layer potentials $\pm u_D$ are produced by a perpendicular displacement field $D$. We use the first-harmonic continuum parameters obtained from first-principles calculations in Ref. [44]. At $u_D=0$, the first-harmonic model has $C_{3z}$ and $C_{2y}\mathcal{T}$ symmetry, plus a pseudoinversion symmetry $\mathcal{I}$ that flips layers and sends $\mathbf{k}\to-\mathbf{k}$ within each valley [45, 46]. A finite displacement field preserves $C_{3z}$ but breaks both $C_{2y}\mathcal{T}$ and $\mathcal{I}$.

Diagonalizing $\hat{h}_K$ gives the dispersion $\epsilon_{\mathbf{k}}$ and Bloch eigenstates $|\psi_{\mathbf{k}}\rangle=e^{i\mathbf{k}\cdot\mathbf{r}}|u_{\mathbf{k}}\rangle$ of the top valence band, which has Chern number $C=+1$ in the $K$ valley. Let $\hat{c}_{\mathbf{k}}^\dagger$ create the Bloch state $|\psi_{\mathbf{k}}\rangle$. Projecting screened Coulomb interactions into this active band gives

where $A$ is the sample area and $\hat{\rho}_{\mathbf{q}}=\sum_{\mathbf{k}}\hat{c}_{\mathbf{k}}^\dagger\langle u_{\mathbf{k}}|u_{\mathbf{k}+\mathbf{q}}\rangle\hat{c}_{\mathbf{k}+\mathbf{q}}$ is the projected density operator at momentum $\mathbf{q}$. The interaction is normal ordered with respect to charge neutrality. We take dual gate screened Coulomb interactions $V_{\mathbf{q}}=e^2\tanh(qd_s)/(2\epsilon_0\epsilon_rq)$ [47], where $\epsilon_r$ is a relative dielectric constant and $d_s$ is an effective gate distance that sets the screening length¹.

We solve Eq. (2) using the infinite density matrix renormalization group (iDMRG) [48, 49, 50] on the infinite cylinder geometry using the TeNPy library [51]. We use a hybrid-Wannier basis and matrix product operator compression techniques developed by some of us [47, 52], and previously applied to continuum moiré Hamiltonians [53, 54, 55] including twisted MoTe$_2$ [20]. We use cylinders with circumference $N_yL_M$, where $N_y=4$–$7$ is the number of moiré unit cells around the cylinder and $L_M$ is the moiré lattice constant. Equivalently, the allowed transverse momenta form $L_y$ evenly spaced “wires” through the Brillouin zone at momenta $k_y=\frac{2\pi}{N_yL_M}(n_y+\Phi_y/2\pi)$, where $n_y$ is an integer and $\Phi_y$ the flux threaded through the cylinder. Periodic and anti-periodic boundary conditions (PBC, APBC) correspond to $\Phi_y=0$ and $\Phi_y=\pi$, respectively. We preserve electric charge $Q_e$, as well as translation symmetry in both the transverse and longitudinal directions. Enlarged MPS unit cells are used when required by the rational filling and/or charge order. We use DMRG bond dimensions up to $\chi=12288$.

Phase diagram

Figure 1c shows our central result: a superconducting dome with nearby metallic, charge-ordered, and fractional Chern insulating phases. Here we discuss the phases and their evolution, while their numerical identification is described in the End Matter and App. D. We first focus on strong, short-ranged interactions, $d_s=1.5\,\mathrm{nm}$ and $\epsilon_r=8$, where the competing phases are clearly resolved. Figure 4 shows that superconductivity persists over a broad range of screening lengths and dielectric constants.

Superconductivity appears over $-0.78\lesssim\nu\lesssim-0.72$ at $u_D=0$ and extends to finite $u_D$ in a dome-shaped region. Decreasing filling towards $\nu=-1$ drives a continuous transition to a metal. At $\nu=-0.8$, the metal persists throughout the $u_D$ range studied and has a circular Fermi surface centered at $\gamma$, consistent with the single-particle band minimum shown in Figure 1b. At still higher hole density, the metal connects to the valley-polarized Chern insulator at $\nu=-1$ found in previous calculations [45, 56, 57, 58, 59, 60, 61] and experiments [2, 3, 6, 7, 43].

At $\nu=-2/3$, we find an FCI, consistent with previous numerical studies [13, 20, 44, 45, 57, 62, 63, 64, 65]. Increasing $u_D$ drives the FCI into a topological CDW with ordering at $\mathbf{q}=\boldsymbol{\kappa}$ ($\sqrt{3}\times\sqrt{3}$ order) and quantized Hall conductivity $\sigma_{xy}=-e^2/h$ (Endmatter).

Doping the superconductor towards $\nu=-2/3$ at any displacement field drives a first-order transition into

a “charge ordered” (CO) regime whose density $\rho(\mathbf{q})$ is peaked at a wavevector $\mathbf{q}\approx\boldsymbol{\kappa}$, connecting to the CDW at $\nu=-2/3$. Endmatter and App. D show its order parameter and extent. The precise nature of the charge order is difficult to determine. At large $u_D$, DMRG on cylinders commensurate with $\sqrt{3}\times\sqrt{3}$ order finds $\rho(\mathbf{q})$ peaked at $\boldsymbol{\kappa}$ exactly. However, at lower $u_D$ or on incommensurate cylinders, the peak is detuned slightly from $\boldsymbol{\kappa}$, and the CO may be incommensurate. Away from $\nu=-2/3$, there is no sign of a two-dimensional many-body gap—indeed a phase which exists over a continuous filling range must be gapless—however, the DMRG proves hard to converge, and so a definitive identification of the gapless charge excitations remains a challenge.

Superconductivity

We first establish superconductivity at one parameter point, then characterize its stiffness, gap, pairing symmetry and topology below. To identify the ground state, we compute correlation functions of three operators: $\hat{c}^{\dagger}_{x,k_y}$, which creates an electron in a hybrid Wannier orbital on the ring centered at $x$ along the cylinder with transverse momentum $k_y$; the ring-averaged charge density operator $\hat{n}_x=N_y^{-1}\sum_{k_y}\hat{c}^{\dagger}_{x,k_y}\hat{c}_{x,k_y}$; and the momentum-resolved Cooper pair operator $\hat{\Delta}_{x,P_y}\sim c_{x,k_y}c_{x,-k_y}$, whose precise form is given in Appendix D.

On a quasi-one-dimensional cylinder, $U(1)$ charge symmetry cannot break spontaneously, so superconductivity manifests as a Luther-Emery liquid (LEL) with quasi-long-range pair correlations [66, 67, 68]. At $\nu=-0.73$, $u_D=0$, and $N_y=6$, we find that the $P_y=0$ pair correlations decay algebraically,

where $\eta_{\mathrm{SC}}$ is the superconducting exponent and $\hat{\Delta}_{x} \equiv\hat{\Delta}_{x,P_y=0}$. Figure 2(a) shows that, as the bond dimension increases, the pair correlations approach a power law. A spinless LEL has power-law correlations in the $(Q_e,k_y)=(2e,0)$ and $(0e,0)$ sectors, while correlations in the single-electron $(1e,k_y)$ and finite-momentum neutral $(0e,k_y\ne0)$ sectors decay exponentially [36, 66, 67].

To confirm this structure, we perform a scaling analysis of the MPS transfer-matrix [50] correlation lengths $\xi_Q$ in each relevant sector. At finite bond dimension $\chi$, algebraic correlations are cut off by the finite MPS correlation length [69, 70]. Following the scaling theory of infinite MPS, we extrapolate $1/\xi_Q$ against the transfer-matrix eigenvalue level spacing $\delta_Q$ [71, 72, 73] (details in App. B). Power-law correlations are signaled by $1/\xi_Q\to0$ as $\delta_Q\to0$, whereas exponentially decaying correlations extrapolate to a finite $1/\xi_Q$. Figure 2(b) shows that $1/\xi_{(0e,0)}$ and $1/\xi_{(2e,0)}$ extrapolate to zero, while $1/\xi_{(1e,k_y)}$ and $1/\xi_{(0e,k_y\ne0)}$ remain finite. Together, these results identify the ground state as a spinless Luther-Emery liquid at fixed circumference.

We next ask whether pairing strengthens toward the two-dimensional limit. For a two-dimensional superconductor compactified onto a cylinder,

where $\rho_s$ is the two-dimensional superfluid stiffness and $v$ is the phase-mode velocity [36, 68, 74]. For $\eta_{\mathrm{SC}}<2$, the superconducting susceptibility diverges at low temperatures as $\chi_{\mathrm{SC}}(T)\sim T^{\eta_{\mathrm{SC}}-2}$ [75]. We find $\eta_{\mathrm{SC}}=1.70$ at $N_y=5$ and 1.29 at $N_y=6$ (see App. § 4.4). The stronger condition $\eta_{\mathrm{SC}}<1$ marks the regime where pairing correlations dominate over competing $2k_F$ density correlations, (a 1d artifact with no 2d limit) [36, 66, 67]. The decrease from $N_y=5$ to 6, shown in the inset of Figure 2(a), is consistent with the $1/N_y$ scaling in Eq. (4), suggesting $\eta_{\mathrm{SC}}<1$ at modestly larger circumference and $\eta_{\mathrm{SC}}\to0$ in the two-dimensional limit.

Superconducting transition temperature.— We now compute the two energy scales that characterize the superconducting state: the charge stiffness $\rho_s$, which controls phase coherence, and the single-electron gap $\Delta_{1e}$, which sets the pair-breaking scale. We first determine $\rho_s$ from the dependence of the ground-state energy density $\mathcal{E}$ on the boundary-condition twist $\delta\Phi_y$. The minimum is at $\Phi_y=\pi$, and expanding to second order gives

Figure 2(e) shows $\rho_s$ as a function of $u_D$, extracted from quadratic fits to Eq. (5) (App. § 4.4). Figure 2(d,e) shows $\rho_s$

vanishes either as $v$ approaches $-2/3$ or as $u_D$ increases. Within the superconducting phase, $\rho_s$ is typically a few tens of $\mu\mathrm{eV}$.

We next estimate the superconducting gap $\Delta_{1e}$ using the plane-wave quasiparticle ansatz for MPS [50, 76, 77]. The quasiparticle gap is obtained from the MPS tangent space at fixed length of five rings. Figure 2(c) shows that the single-electron gap remains finite, with $\Delta_{1e}\lesssim0.5\,\mathrm{meV}$, while the charge-2e gap $\Delta_{2e}$ is consistent with zero, as expected for a superconductor with gapped single electrons and gapless Cooper pairs.

These energy scales constrain the superconducting transition in two distinct ways. Phase coherence is lost at the BKT transition, where $k_BT_{\mathrm{BKT}}=\pi\rho_s(T_{\mathrm{BKT}})/2$ [78], while the single-electron gap $\Delta_{1e}$ sets the scale at which pairing itself is destroyed. The superconducting transition is limited by the minimum of the phase-ordering and pair-breaking scales. For $d_s=1.5\,\mathrm{nm}$, $\Delta_{1e}\sim0.5\,\mathrm{meV}$ and $\rho_s\sim30\,\mu\mathrm{eV}$ correspond to several kelvins and $\sim500\,\mathrm{mK}$, respectively. Therefore the transition temperature is limited to $T_c\lesssim500\,\mathrm{mK}$ at this point.

$f$-wave topological superconductivity.— We next determine the pairing symmetry and topology of the superconductor. In a two-dimensional superconductor, the pairing symmetry is encoded in the pair wavefunction $\Delta(\mathbf{k})=\langle\hat{c}_{\mathbf{k}}\hat{c}_{-\mathbf{k}}\rangle$, which vanishes identically because our iDMRG preserves charge $U(1)$. We instead measure a proxy pair wavefunction $\widetilde{\Delta}(\mathbf{k})$ from the charge-2e operator with the slowest asymptotic decay along the cylinder (see App. C) [36, 79] [35]. Figure 3(a) shows that the phase of $\widetilde{\Delta}(\mathbf{k})$ winds by $-6\pi$ around the Fermi contour, corresponding to angular momentum $\ell=-3$ in our convention. This establishes chiral $f$-wave pairing. If the superconductor were in the weak-pairing limit, its topology would follow from the pair winding. For pairing in a Chern band [36, 80, 81, 82],

where $C$ is the normal-band Chern number and $\ell$ is the pair winding around the Fermi surface. In this limit, the number of chiral Majorana edge modes is $2c_{-} = C_{\mathrm{BdG}}$, giving $2c_{-} = -1$ for $C = 1$ and $\ell= -3$. Here $c_{-}$ is the chiral central charge quantized in half-integer units [80, 83]. Beyond weak pairing, however, pairing symmetry alone does not determine the topology [80, 84].

We therefore diagnose the topology directly within number-conserving DMRG [84, 85]. An odd $2c_{-}$ implies that a superconducting vortex binds an unpaired Majorana zero mode and has non-Abelian statistics, whereas an even $2c_{-}$ gives Abelian vortices [80, 83]. On a quasi-one-dimensional cylinder, this even–odd distinction appears as a gapless analogue of one dimensional symmetry protected topological order and can be detected by a nonlocal string operator [84, 86],

where $\hat{N}_{x'}$ is the total electron number on the transverse ring at $x'$. For zero-momentum pairing, bosonization gives [84]

with $\eta_b = \eta_{\mathrm{SC}}/4$. Figure 3(b) shows algebraic decay with $\eta_b \approx0.41$, while the Cooper-pair exponent $\eta_{\mathrm{SC}} \approx1.65$ satisfies $\eta_{\mathrm{SC}}/4 = \eta_b$ within fitting error. This establishes an odd number of chiral Majorana edge modes and hence non-Abelian vortices.

Screening-controlled superconductivity.— Finally, we examine how the range and strength of the Coulomb interaction controls superconductivity. We parameterize the interaction strength at hole density $n = |\nu|/A_M$ as $U \equiv nV_{q=0}(d_s,\epsilon_r)/2$ where, as above, $d_s$ is an effective screening length related to, but not necessarily equal to, the geometrical gate spacing. As shown in Figure 4(a), weak interactions favor a metal, while stronger repulsion generally favors the charge order, with the topological superconductor in between. Lines of constant $\epsilon_r$ indicate representative parameters for moiré devices. Figure 4 establishes that superconductivity can reach $d_s \approx5$ nm. However, numerical convergence becomes more difficult at larger $d_s$.

Discussion

We have established repulsion-driven chiral $f$-wave topological superconductivity with non-Abelian vortices in a microscopic continuum model of MoTe$_2$ with gate-screened Coulomb interactions. The superconductor appears at $-0.78 \lesssim\nu\lesssim-0.72$ at zero displacement field and survives at finite $u_D$, showing that it does not rely on inversion symmetry. Whether the pairing originates from a Kohn–Luttinger instability, strong-coupling anyon superconductivity, or another mechanism remains an open question. However, the increase of $\rho_s$ with doping away from $\nu=-2/3$ (Figure 2(d)) is consistent with doped-FCI scenarios.

Our numerics make several concrete experimental predictions. Increasing hole doping away from $\nu=-2/3$ enhances the superfluid stiffness, while increasing $u_D$ strongly suppresses it even as the single-electron gap remains comparatively robust. Microwave measurements of $\rho_s$ [87, 88] together with tunneling spectroscopy [89] can directly test these trends. Phase-sensitive measurements that account for the spin-valley, finite-momentum, and Bloch-orbital structure of the intravalley pairs can test the predicted chiral $f$-wave pairing [90, 91], while vortex spectroscopy and thermal transport can potentially probe its topology [92, 93, 94]. Electrostatic screening is another experimentally tunable knob, driving the system from charge order through the TSC to the metal (Figure 4). It can be modulated experimentally by varying the gate distance or dielectric environment [95, 96].

The charge-ordered regime provides a possible microscopic origin of the reentrant Hall state. At $\nu=-2/3$, increasing $u_D$ drives a transition from the FCI to a topological CDW. In experiments, loss of valley polarization (which is beyond the scope of our single-valley calculation) may preempt this transition [43]. Upon doping away from $\nu=-2/3$, we continue to find near-$\sqrt{3}\times\sqrt{3}$ charge order, although DMRG convergence becomes more challenging and its detailed nature is less certain. In this doped regime, disorder may Anderson localize the residual carriers, exposing the integer Hall response of the topological CDW, thus providing a route to the reentrant integer quantum anomalous Hall (RIQAH) regime [5, 97]. This charge order regime is distinct from the translation-symmetric $\mathbb{Z}_3$ orthogonal metal [98, 99, 100], or coexistence of FCI, charge-ordered, and superconducting domains [101] and thermally stabilized non-Fermi liquid metals [40, 41, 42, 102, 103]. These scenarios can be distinguished experimentally through imaging and transport. The predicted near-$\sqrt{3} \times\sqrt{3}$ modulation can be directly detected in microscopy [104, 105], while the $\mathbb{Z}_3$ orthogonal metal, domain-coexistence, and charge-order scenarios give distinct trajectories in the $(\rho_{xy}, \rho_{xx})$ plane [99, 100, 101].

Note Added.— During the preparation of this work, Refs. [106, 107, 108] appeared. Ref. [106] finds $p$-wave rather than $f$-wave superconductivity at similar fillings in a related tMoTe$_2$ model, while Refs. [107, 108] find broader families of chiral superconducting states in minimal Chern-band models with different pairing symmetries.

Acknowledgments

We thank T. Senthil and Zijian Wang for helpful discussions. Our computations were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. Tianle Wang acknowledges support from the Schwinger Foundation Postdoctoral Fellowship (JSF-25-09-0004). N.Y.Y. acknowledges support from a Brown Investigator Award. M.P.Z. and Tianle Wang acknowledge support from the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division under Contract No. DE-AC02-05CH11231 (Theory of Materials program KC2301). D.E.P. acknowledges support from NSF CAREER Award No. DMR-2542485. Taige Wang and A.V. acknowledge support from the Simons Collaboration on Ultra-Quantum Matter, a grant from the Simons Foundation (651440, A.V.). Taige Wang acknowledges additional support from the Harvard Quantum Initiative Fellowship. The Flatiron Institute is a division of the Simons Foundation.

End Matter

Topological CDW and its transition to the FCI.— To study competition between charge order and topological order at $\nu=-2/3$, we performed DMRG at $N_y=6$ in a larger variational space whose translation symmetry was reduced to permit enlarged unit cells with $\sqrt{3} \times\sqrt{3}$ order. At $u_D=0$ we find the well-studied Laughlin FCI, while at high displacement field we find a state with exponentially decaying correlations in all channels, signaling a many-body gap. Figure 5(a) shows this state’s real-space density $\rho(\boldsymbol{r})$, where a clear $\sqrt{3} \times\sqrt{3}$ modulation is directly

visible. Related commensurate charge-ordered states have been found previously in exact-diagonalization studies of twisted MoTe$_2$ [17, 63, 109, 110], including a study of their real-space density in Ref. [109].

The density forms well-localized “fidget spinner” motifs, akin to the Wannier functions in the topological bands of twisted bilayer graphene [111], suggesting that the charge order may retain nontrivial topology. To confirm this, we thread $2\pi$ magnetic flux through the cylinder, and measured the pumped charge $\Delta Q$, which reflects the Hall conductivity $\sigma_{xy}=-(e^2/h)\Delta Q$. Figure 5(b, inset) shows the expected $\sigma_{xy}=-\frac{2}{3}e^2/h$ for the FCI at $u_D=0$ and $\sigma_{xy}=-e^2/h$ in the CDW, confirming it as a topological CDW (TCDW). Charge pumping across a range of displacement fields, Figure 5(b), places the FCI/TCDW transition near $u_D\approx0.45\,\mathrm{meV}$.

Charge Order Parameter.— Away from $\nu=-2/3$, charge order persists over a finite range of fillings, manifesting through peaks in the Fourier space charge density $\rho(\boldsymbol{q})$ away from the moiré reciprocal lattice. As shown in Figure 13 (App. § 4), $\rho(\boldsymbol{q})$ develops pronounced peaks near the six $\boldsymbol{\kappa}_j$ points at the corners of moiré Brillouin zone, showing that the order remains close to $\sqrt{3}\times\sqrt{3}$ even at incommensurate filling. For all dopings, we observe non-trivial peaks in $\rho(\boldsymbol{k})$ only in the vicinity of $\boldsymbol{\kappa}$.

We quantify this charge order by

The allowed momenta $\boldsymbol{q}$ depend on the filling-dependent DMRG unit cell on the cylinder, and do not always include $\boldsymbol{\kappa}_j$. On $\boldsymbol{\kappa}$-commensurate cylinders and at larger displacement fields, the peak lies exactly at $\boldsymbol{\kappa}_j$, while at lower displacement field or on incommensurate cylinders it can shift away. As noted above, the state should be gapless away from $\nu= -2/3$, but DMRG with expanded unit cells is expensive and suffers from convergence issues, preventing a clear determination of the nature of any gapless charge excitations. We leave a definitive characterization of this phase to future work.

Phase identification

Fig. 6 summarizes the diagnostics used to distinguish the TSC from its competing metallic and charge-ordered phases in Figure 1(c) and Figure 4. In both cases, the TSC is sandwiched between a metal on one side and charge order on the other.

We identify charge order using $\rho(\kappa)$ from Eq. (9), classifying points with $\rho(\kappa) \ge0.075$ at $\chi= 2048$ as charge ordered. To distinguish the metal from the TSC, we use the single-electron transfer-matrix correlation length. Following the finite-entanglement scaling procedure used in the main text and described in detail in App. § 4, we define $(\widetilde{\xi}_Q)^{-1} \equiv\lim_{\chi\to\infty} \xi_Q(\chi)^{-1}$. Gapless charge-$Q$ excitations give $(\widetilde{\xi}_Q)^{-1} = 0$, while a finite value signals exponentially decaying correlations in that sector. Among points without charge order, we identify the metal by $(\widetilde{\xi}_{1e})^{-1} \le0.045$.

The two criteria above are necessary but not sufficient to identify a superconductor. For converged states that are non-metallic and non-charge ordered, we verify superconductivity using the diagnostics established in Figure 2, including algebraic Cooper-pair correlations, a gapped single-electron sector, and gapless charge-$2e$ correlations. We identify this region with the TSC characterized in Figure 2 and Figure 3. Figure 6 shows $\rho(\kappa)$ grows upon approaching charge order, while $(\tilde{\xi}_{1e})^{-1}$ decreases toward zero upon approaching the metal.

Superconducting correlations.— Figure 7 shows the inverse superconducting exponent $\eta_{\mathrm{SC}}$ across both phase diagrams, defined via $\langle\hat{\Delta}^{\dagger}_{x}\hat{\Delta}_{0}\rangle\sim r^{-\eta_{\mathrm{SC}}}$ and computed via scaling collapse fit (App. § 4.4). This is expected to be non-zero in a metal, where its value is set by the gapless electron-electron correlations. As discussed in the main text, $\eta_{\mathrm{SC}} < 2$ gives a divergent superconducting susceptibility, while $\eta_{\mathrm{SC}} < 1$ marks the regime where Cooper-pair correlations dominate over competing $2k_{F}$ density correlations. The estimates of $\eta_{\mathrm{SC}}$ across the phase diagram are less precise than the detailed finite-entanglement analysis in Figure 2, so we do not use them to determine the phase boundaries, but they are sufficient to track the evolution of pairing correlations.

Throughout the TSC region, the Cooper-pair correlations remain strong. On entering the charge-ordered region, the fitted $\eta_{\mathrm{SC}}$ increases rapidly and the pair correlations become short ranged. Toward the metal, $\eta_{\mathrm{SC}}$ instead increases more gradually as the Cooper-pair correlations weaken. Figure 7(c,Figure 7d) shows these trends along two cuts through panel (a), at fixed $\nu= -0.73$ and fixed $u_{D} = 0$, respectively. The behavior of $\eta_{\mathrm{SC}}$ is consistent with the superconducting region identified independently from charge order and the single-electron correlation length.

References

Microscopic Model

We consider a simple microscopic model for tMoTe$_2$: a continuum model for the valence bands with screened Coulomb interactions. We now provide a self-contained description of the model.

Lattice Geometry

Our conventions for the TMD geometry follow [15]. We take monolayer TMD lattice vectors

where $a_0=3.52$ Å is the monolayer TMD lattice constant [44]. We focus on $\theta=3.89^\circ$.

Two monolayers with relative twists of $\theta/2$ and $-\theta/2$ for the top and bottom layers form a moiré pattern with a superlattice generated by

where $L_M$ is called the moiré length. The minimal length reciprocal vectors are

where $R_\theta$ is the counterclockwise rotation matrix. We consider the reciprocal lattice generated by

The high-symmetry points of the moiré lattice, shown in Figure 8(b), are at

Brillouin Zone Geometry

Our DMRG conventions follow [47] and [20]. For DMRG it is convenient to use a rectangular Brillouin zone. We use $\mathcal{B}=[-G_x/2,G_x/2]\times[-G_y/2,G_y/2]$ where

so that $\mathbf{a}_1=2\pi\left(\frac{1}{G_x},\frac{1}{2G_y}\right)$.

Instead of a Monkhorst-Pack grid, we instead take a rectangular grid of size $L_x\times L_y$, where we will take $L_x\to\infty$. Define grid basis vectors

We then take a $k$-grid

where $0\leq\Phi_{x,y}\leq2\pi$ and

For example, for $N_y=5$, $n_y=-2,-1,0,1,2$ and for $N_y=6$, $n_y=-3,-2,-1,0,1,2$. When $L_x\to\infty$, this gives $L_y$ evenly spaced “wires” in the Brillouin zone, shown in Figure 8(c).

The phase $\Phi_y$ controls the boundary conditions. If $\Phi_y=0$, then one wire will go through $\Gamma$, giving periodic boundary conditions (PBC); $\Phi_y=\pi$ corresponds to anti-periodic boundary conditions (APBC). We use APBC unless otherwise stated.

To understand the real-space geometry of the whole system, we define

which satisfy $L_a\cdot\delta k_b=2\pi\delta_{ab}$. The whole system is then a torus spanned by $\mathbf{L}_x,\mathbf{L}_y$, making it rectangular. We require that the unit cell torus spanned by $\mathbf{a}_1,\mathbf{a}_2$ tiles the whole system, which imposes the requirement

We may then take $L_x\to\infty$, turning the torus into an infinite cylinder.

Single-Particle Hamiltonian

We use a standard continuum model [15] for the valence bands of tMoTe$_2$, which we reiterate here for the reader’s convenience. The model is identical to the one from previous DMRG work [20]. The kinetic term for the $K$-valley acts in the bottom/top layer basis $(\ell=t,b)$ as

where $u_D$ is the vertical displacement field,

is the quadratic band maxima of the TMD valence bands,

is the on-layer moiré potential, and

We take parameters from first principles calculations of MoTe$_2$ [44]: $(a_0, m^*, V, \psi, w) = (3.52\ \text{Å}, 0.6m_e, 20.8\ \text{meV}, -107.7^\circ, -23.8\ \text{meV})$ at $\theta= 3.89^\circ$. At these parameters, the top valence band for the $K$ valley has Chern number $C = +1$, and the second-to-top band has Chern number $C = -1$. Within a valley, the model has $C_{3z}$ and $C_{2y}T$ symmetry, as well as moiré translation symmetry.

We solve the kinetic Hamiltonian in the Bloch basis

where $b$ labels bands and $\hat{T}_{\mathbf{a}_j}$ applies a (passive) translation by $\mathbf{a}_j$. In practice, we use a plane-wave basis with six hexagonal shells of planewaves, as needed for convergence. Figure 9 shows the bandstructure, dispersion, and real-space charge density of the top band. Additionally, it computes the Berry curvature and trace condition violation across the Brillouin zone finding, as in prior work [20], that the band is perturbatively close to ideal at these parameters.

Many-Body Hamiltonian

We define electron creation operators $\hat{c}_{k b}^{\dagger}$ that create the Bloch states $\lvert\psi_{k b}\rangle$. We consider an interacting Hamiltonian

where

• $:\hat{O}:$ denotes normal ordering with respect to charge neutrality (where all the valence bands are fully filled)

• the kinetic term is

• in terms of the system area $A$ and relative momentum $\mathbf{q}$ (not restricted to the Brillouin zone), the screened Coulomb interactions are

where $\epsilon_r$ is the relative dielectric constant and $d_s$ is both the sample-to-top-gate and sample-to-bottom-gate distance [47]

• the density operators are defined as

in terms of the form factors

For numerics, we project this Hamiltonian into the top valence band in the $K$ valley. We have chosen conventions so that this projection is exactly the restriction of the band sums to a single index. In terms of the dispersion $\epsilon_k$, form factors $\lambda_q(k)$, and creation operators for the top band, the projected Hamiltonian is

Computational Basis

Following previous DMRG work on moiré systems [20, 47, 112], we take a computational basis of hybrid Wannier functions that are maximally localized along $x$ and planewaves along $y$. In terms of the Bloch wavefunctions $|\psi_k\rangle$ for the top band, these are defined as

The gauge $\zeta(\mathbf{k})$ is chosen such that [47, 112]

where $\hat{x}$ is the position operator and $\hat{P}=\sum_k|\psi_k\rangle\langle\psi_k|$ is the (single-particle) projector to the top band. Note that for wires with flux $0\leq\Phi_y\leq2\pi$, we have

interpolating between periodic and anti-periodic boundary conditions. Here the polarizations are the centers of the Wannier orbitals along the $x$-axis: [47]

(Here we have set the electric charge to be $e=1$.) The modern theory of polarization [113, 114, 115] shows that the Chern number of the band is given by the winding of the polarization through the Brillouin zone as $C=\int\frac{dP_x}{dk_y}dk_y$, which here is $C=1$. Figure 10 shows the winding the polarization, the exponential decay of the Wannier functions along $x$, and their real-space profiles.

At the many-body level we define Wannier creation operators $\hat{c}^{\dagger}_{n,k_y}$ via (note minus signs)

In the main text and below we write the computational basis as

where $x=x_n=\mathbf{e}_x\cdot(\mathbf{a}_1n)$ is the $x$-position of the center of the $n$th ring. Note that the exponential factor $e^{+i\mathbf{k}\cdot(n\mathbf{a}_1)}\neq e^{ik_xx}$, but instead has an extra phase. The many-body Hamiltonian Eq. (35) can be readily expressed in this basis. This puts the tMoTe$_2$ Hamiltonian onto a cylinder geometry with a circumference of $L_y$ unit cells that is infinite along $x$ but composed of maximally-localized hybrid Wannier functions. This choice ensures that the kinetic hoppings decay quickly along $x$. However, the quasi-long range screened Coulomb interactions still decay slowly, requiring special techniques to perform DMRG in practice.

MPO Compression

To perform infinite DMRG, one must express the Hamiltonian as an infinite matrix product operator (MPO), which for fermions is a tensor network of size $2\times2\times D\times D$ where $D$ is the virtual bond dimension. Naive constructions of MPOs truncate the interaction

range in the $\hat{d}^{\dagger}_{n,k_y}$ basis at some cutoff length $M$. Not only does this cutoff cause errors, potentially biasing the Hamiltonian, but it requires huge bond dimensions $D\sim2000$ for only $M=2$ at $L_y=6$ [47]. To avoid these issues, we apply the technique of iMPO compression developed by some of us [47, 52], which faithfully represents the full Hamiltonian including its quasi-long range interactions. In practice we check that each matrix element is represented with error $\epsilon<1\times10^{-3}$ meV out to separations of 200 sites along the 1d chain. The resulting iMPOs have bond dimension $D\sim120$ for $L_y=6$, and grow modestly with $L_y$.

Practical Application of The Scaling Theory of iMPS

This Appendix briefly reviews the scaling theory of matrix product states and how it may be applied in practice to extrapolate correlation lengths to bond dimension $\chi\rightarrow\infty$ in a principled manner. The essential idea is that at finite $\chi$ all MPS correlation functions decay exponentially $C_\chi(x)\sim e^{-x/\xi}$—making it impossible to encode a power law prefactor $C_\eta(x)\sim x^{-\eta}e^{-x/\xi}$ exactly, including the case $\xi= \infty$. However, the limit

is well-controlled in the right variable $\delta(\chi)$. The classic references on this topic, from which we draw extensively below, are [71, 72, 116, 117].

Suppose we have an injective iMPS of bond dimension $\chi$ with transfer matrix

with $(L_i|R_j) = \delta_{ij}$ and $1 = \lambda_0 > |\lambda_1| \ge|\lambda_2| \ge\cdots$. Here and below rounded kets $|v)$ denote vectors in the doubled virtual space that the transfer matrix acts on. Figure 11(a) shows an example transfer matrix spectrum.

In terms of this, the connected correlation function

takes the form

where $w_j$ are form factors determined by $\hat{O}$, $|R_j)$ and $|L_j)$. We decompose the eigenvalues into their magnitude and phase as

where $a$ is the lattice spacing, $\epsilon_j \ge0$, and $\xi_j = 1/\epsilon_j \ge0$ is called the correlation length. Such a correlation function decays exponentially: $C_O(x) \xrightarrow{x\gg1} w_1\lambda_1^x = w_1e^{-i\phi_1x}e^{-x/\xi_1}$. We note the phase $\phi$ contains a wealth of physical information. For instance, for $(1e,k_y)$, it is related to the $k_x$ values where the $k_y$ wire intersects the Brillouin zone.

Correlation function often takes the form $C_\eta(x) \sim e^{-x/\xi}x^{-\eta}$, which cannot be represented as any finite sum of exponentials, and is therefore beyond an exact MPS representation at finite $\chi$. However, they can be represented by a continuous integral over decaying exponentials. Concretely,

which follows from the definition of the $\Gamma$ function. As $\chi\to\infty$, the points $\lambda_j$ must therefore accumulate into lines in the unit disk. We refer to these lines as “branches”, as they produce the branch cuts in the Stieltjes transform representation of power laws like Eq. (49).

To quantify this, it is standard to define scaling variables with units of inverse length:

for particular $\lambda_i$ and $\lambda_j$. As $\chi\to\infty$, $\delta_{ij}(\chi) \to0$. This can be thought of as the small length scale at which the power laws is gives in to the exponential decay [116]. While appealing, this leaves open an obvious practical question: which $\lambda_i$ and $\lambda_j$ should be used?

To illustrate the difficulty, suppose for an operator $\hat{O}$ with charge $Q$, the true ground state correlation takes the form:

which has both multiple power laws (potentially with different phases $\phi_j$) and multiple exponential terms. Each of these functions will give rise to a different branch in the unit disk on which eigenvalues will accumulate. Our practical prescription is then as follows: isolate the eigenvalues belonging to each branch in the correlation function via a careful inspection of the transfer matrix eigenvalue spectrum and take $\delta_{ij}$ corresponding to separation between the largest ones.

We note several practical issues. First, $\lambda_i$ and $\lambda_j$ should necessarily be from the same symmetry sector. This is easy to enforce when the symmetry is conserved, but additional non-conserved symmetries (perhaps non-Abelian ones) can also act on the transfer matrix spectrum. Second, the phase of $\lambda_i$ is often helpful in identifying a branch, but multiple branches can have the same phase. Third, if there are “bystander” exponential correlations like $E_1\lambda_e^x = E_1e^{i\phi_1 x}e^{-x/\xi_1}$ with long correlation lengths, then the power law eigenvalues $\lambda_j$ will cross $\lambda_e$ as $\chi\to\infty$. The index $j$ for the scaling variable $\delta_{ij}$ changes at this point, and we must use continuity to assign $j$ correctly. Fourth, even within the eigenvalues $\lambda_j$ for a particular power law branch, there may be finite-$\chi$ eigenvalue crossings. For instance, with the fourth and fifth eigenvalues swapping places between $\chi= 2048$ and $\chi= 4096$. An avoided crossing will create deviations from the scaling behavior, while an unavoided crossing will create a labeling issue. Only once such issues are accounted for can a practical identification of $\delta_{ij}(\chi)$ be made for each $\chi$ and used for extrapolating to $\chi\to\infty$. This is shown in Figure 11(c) where, for instance, several eigenvalue crossings are visible in the higher correlation length values.

The fits used in practice are shown in Figure 11(d). Other points are fit using a similar procedure. As $\xi^{-1} \ge0$, the slightly negative extrapolated values like $-0.0023$ are an error incurred by the extrapolation procedure.

Proxy Pair Wavefunction

In this section we define the “proxy pair wavefunction” defined in the main text and argue that it reflects the pairing symmetry of the superconductor in the 2d limit. This section draws from, but also expands upon, previous work from some of us on this topic [36, 79] [35].

Consider the 2d Cooper pair operator

and let $\hat{\Delta}(\boldsymbol{r};\boldsymbol{k})$ be the Fourier transform in its relative coordinate. The ground state wavefunction of a 2d superconductor spontaneously breaks electromagnetic symmetry $U(1)_Q$, whose order parameter is the pair wavefunction

whose value is $\boldsymbol{r}$-independent for pairing at zero momentum, which we assume throughout. The transformation of the pair wavefunction $\Delta(\boldsymbol{k})$ under the space-group, particularly rotations, determine the superconducting pairing channel. Within DMRG, however, we only have access to cylinder wavefunctions $\Psi(L_y,\chi)$ that preserve $U(1)_Q$. How may $\Delta(\boldsymbol{k})$ be probed from that data?

We briefly recapitulate the argument from [79], eliding details to streamline the discussion. In the far-field regime $|\boldsymbol{r}| \gg|\delta\boldsymbol{r}|$, cluster decomposition (indirectly) implies

which is an expression of long-range order. (Here and below, $=$ is equality in the far-field regime.) If we fix a particular $\boldsymbol{r}$ in the far-field regime, we may therefore compute the pair wavefunction up to an overall constant as

If we now go to finite-but-large circumference $L_y$, long range order becomes quasi-long range order

For sufficiently large $L_y$, one expects the limit to be smooth, with

in the far-field regime $x \gg\delta x$ and, of course, $0 < y < L_y$. Finally, at finite but large bond dimension, one expects

where $\delta_{2e}(\chi)$ is the scaling variable from Sec. B that approximately sets the crossover from power law to exponential. That is, in the intermediate regime of far field by less than $1/\delta_{2e}$, the finite bond dimension correlation function matches the infinite bond dimension correlation function. We have thus argued

As $\eta(L_y)$ and $D(r;\delta r|L_y,\chi)$ are directly computable from DMRG data, one may use them to approximate $\Delta(0;\delta r)$ and, by Fourier transforming, $\Delta(k)$. We emphasize the above is an argument where error bounds are not made precise, though each step is physically reasonable.

We now move beyond previous references. By Eq. (47), the four-point function takes the form

where $\lambda_j$ are the transfer matrix eigenvalues in the $2e$ sector, which depend on $L_y$ and $\chi$, but are independent of $r=(x,y)$ and $\delta r$. In the power-law regime, the sum is dominated by the eigenvalues from the branch of the transfer matrix spectrum responsible for the power law:

One property is immediately obvious: as these terms limit to a single power law, each of the $w_j(y,\delta r)$’s must transform under the same irrep of the space group. A single $w_j$, in particular $w_1$, is therefore sufficient to determine the pairing symmetry.

We therefore define the “proxy pair wavefunction” to be

which is Fourier transformed in the relative coordinate to yield $\widetilde{\Delta}(k)$ in the main text. We have checked that using $w_2$ results in a small quantitative error, confirming the robustness of this proxy. Further, $\widetilde{\Delta}(\delta r)$ is significantly cheaper to compute than $D(r;\delta r|L_y,\chi)|_{x|^{-\eta(L_y)}}$, which was used in previous work [79].

DMRG Details & Phase Identification

In this section we describe the details of how each phase is identified. We first describe the DMRG computations that were performed, and the observables used for their analysis.

DMRG Computations

DMRG computations were performed at a variety of fillings, displacement fields, screening lengths, and dielectric constants, as shown in the figures. To permit a filling $\nu=p/q$, the iDMRG unit cell must obey $q|N_xN_y$ (i.e. $N_xN_y \pmod q=0$), necessitating different choices of $N_x$ for different fillings. Some phases also required special constraints. For example, Laughlin FCIs require $3|N_x$ to permit all three root states. In practice, our representative points use $N_x=5$ and $N_y=6$, giving 30 total sites within the iDMRG unit cell, which supports filling $\nu=-22/30=-11/15\approx-0.73$.

DMRG computations were performed at each parameter point as various bond dimension $\chi$, reaching at least $\chi=2048$ and sometimes as large as $\chi=12288$. (Used for example at largest $d_s$.) For most points, two sets of DMRG computations were performed: those conserving the maximal translation symmetry, and those preserving only a subgroup of translation symmetry so as to allow $\sqrt{3}\times\sqrt{3}$ order to develop. For instance, at $N_y=6$, these would have translation by one unit cell around the cylinder, conserving momentum 6-fold momentum, versus translation by three unit cells around the cylinder, conserving 2-fold momentum. The Hilbert space with reduced translation symmetry (RTS) is a strict superset of the maximal translation symmetric (MTS) space, so the variational energy is expected to be strictly smaller from DMRG in this space at the same bond dimension. This is indeed observed in practice. However, the reduced translation symmetry makes the computations substantially more expensive, even at the same bond dimension. We therefore use the RTS states to analyze the formation of charge order, but perform the analysis of non-charge ordered states with MTS states.

DMRG Observables

For each parameter and each RTS or MTS, we monitor the convergence of the state across a variety of difference convergence parameters and physical observables to ensure it is well-converged and meaningful. Phase assignment is then made computing a basket of physical observables, which are examined manually to ensure they are numerically and physically reasonable for each point. The most important observables are: the correlation lengths in the $0e$, $1e$, and $2e$ sectors, whose extrapolation was described in App. B above; the momentum-space occupation $n(k)$; the real-space charge density $\rho(r)$; its Fourier transform $\rho(q)$; the structure factor $S(q)$; and the Cooper-Cooper correlation function $D$.

Most of these are simple. For instance,

which is easily expressed in terms of the computational basis using Eq. (42). Similarly,

where, due to Poisson summation, the only allowed momentum points are at $q=G$ in the reciprocal lattice determined by the translation symmetry of the state. The structure factor is more complex, but was described in detail in Ref. [20].

To define the pair correlation function $D$, we first define a Cooper operator

where $X$ is the center of mass position, $x$ is the relative position between the pairs, and $k_y$ is their momentum. (We specialize to pairing at zero total momentum.) In terms of these, we define the Cooper-Cooper correlation function

which is defined for $X>(x+x')/2$ so the annihilation operator is strictly to the right of the creation operator. (Outside this regime, various minus signs are required, whose details we elide.)

In a superconductor (or metal), one expects that at least some components of this correlation function have power-law behavior in the large bond dimension limit: $D(X)\sim X^{-\eta_{\mathrm{SC}}}$. However, individual components do not necessarily behave in such a nice manner. For instance, they can have additional oscillations $e^{i\zeta X}X^{-\eta_{\mathrm{SC}}}$, or can be a mix of exponentially decaying and power law decaying correlations in different $(x,k_y)$ channels. Physically, however, one expects there to be a ”Cooper operator” which is some particular superposition $\Delta(X)=\sum_{x,k_y}u_{x,k_y}\widehat{\Delta}(X,x,k_y)$ that is a low-energy excitation in the superconductor whose structure should dominate the long-range correlations.

Based on such intuition, we define superindices $I = (x,k_y)$ and $I^{\prime} = (x^{\prime},k_y^{\prime})$, whereupon the Cooper-Cooper correlation function becomes a matrix $[D(X)]_{I,I^{\prime}} = D(X,x,k_y,x^{\prime},k_y^{\prime})$ at each $X$. We then perform a singular value decomposition

where $\sigma_j(X)$ are the singular values of $D(X)$. In practice we therefore plot

the largest singular value, as a function of the center of mass separation. The main text plots labelled

$\langle\hat{\Delta}(x)\hat{\Delta}(0)\rangle$ indeed plot $\sigma_0(x)$. We find that the corresponding singular vector $u_0(X)$ converges quite quickly with $X$ in practice, in line with the above intuition.

Phase Identification

We focus on three representative parameter points whose only difference is the relative dielectric constant $\epsilon_r$, while screening length $d_s$, filling, etc are held fixed. The parameters of the points are described in Table I.

Point	$\epsilon_r$	$d_s$ (nm)	$U$ (meV)	Phase Identification
A	5.04	1.4	39	Charge Order
B	7.56	1.4	26	Superconductor
C	15.12	1.4	13	Fermi Liquid

Figure 12 shows the key correlation functions and the extrapolation the correlation lengths in various sectors across the three parameter points. The Fermi liquid is clearly identified: its 1e and 2e correlation functions have clear power-law behavior, and the inverse correlation lengths $\xi_{0e}^{-1}$, $\xi_{1e}^{-1}$, $\xi_{2e}^{-1}$ all vanish upon extrapolation, as expected. The superconducting data shown here recapitulates the phenomenology of Figure 2.

Finally, the point at $\epsilon=5.04$ shows much shorter-range correlation in 1e and 2e, and only $\xi_{0e}^{-1}$ vanishes for RTS state. We interpret the state as forming a charge order that gaps out the 1e and 2e charge excitations. We note that the unit cell size here does not allow $\sqrt{3}\times\sqrt{3}$ order. Indeed, the characteristics of this state are consistent with a gapped one dimensional phase with a unit-cell-size dependent density wave, whose thermodynamic limit is ill-defined. This has the effect of gapping out charge excitations, potentially spuriously. Future work will examine states where $\sqrt{3}\times\sqrt{3}$ order is possible.

Figure 13 shows key observables in each state. For the Fermi liquid, the occupations change sharply from $n(k)=0$ to $n(k)=1$ across a circular Fermi surface centered at $\gamma$. The structure factor look roughly as expected for a Fermi liquid, and the charge density is peaked strongly at the MX and XM centers shown in Figure 9(c), with no additional translation breaking.

The superconductor looks qualitatively similar to the metal in these observables, which is to be expected in a quasi-1d system. The occupations $n(k)$ show the clearest distinction, with a relatively smooth crossover from empty to occupied near the Fermi surface. Here our phase identification therefore depends more strongly on the extrapolated correlation lengths.

The charge-ordered state, by contrast, is immediately apparent from $\rho(r)$. Due to the $(N_x,N_y)=(5,6)$ iDMRG unit cell, the $\kappa_{t,b}$ points are not accessible momenta. However, $\rho(q)$ is strongly peaked in the momenta closest to $\kappa$, signaling charge order. As the density is incommensurate with a $\sqrt{3}\times\sqrt{3}$ charge order, the peak in $\rho(q)$ is slightly offset from $\kappa$ to realize a charge order commensurate with $(N_x,N_y)$. At accessible bond dimensions, correlations in the $Q=1e$ sector remain gapped. (We note that going beyond $\chi=2048$ incurs significant expense when translations are broken to permit charge order.) The spatial patterns look somewhat similar to the fidget spinner patterns in Figure 5, except that many spinners now have a “dimer”-like center with two peaks or a “water molecule”-like center with three peaks. The strange and inconsistent features of this state will be the topic of future work.

Finally, Figure 14 provides a three-dimensional plot of the occupations across the Brillouin zone in each phase.

Plane Wave Excitations

To obtain the charge gap for 1e and 2e excitations, we computed the excitation spectrum $\omega_{\pm}(k)$ using the plane-wave quasiparticle ansatz MPS [50, 76], with $\pm$ denoting the addition and removal of respective particles. The momentum of the excitations $k=(k_x,k_y)$ is determined by the MPS charge number $k_y(n_y)$ and phase $k_x=(\theta+j)/N_{\mathrm{uc}}\ \mathrm{mod}\ 1$, $j=0,1,2,N_{\mathrm{uc}}-1$, where $\theta$ is the phase difference between sites of the plane wave superposition ansatz, and $N_{\mathrm{uc}}$ is the size of the enlarged unit cell. The charge gap is consequently defined as $\Delta=\min_k\omega_+(k)-\max_k\omega_-(k)$.

Curve Fitting Details

Stiffness

The superfluid stiffness follows from the second-order change in the energy density given in Eq. (5), which we reproduce here for convenience:

where $\delta\Phi_y = \Phi_y - \pi$. In practice, we compute the DMRG ground state at a range of $\delta\Phi_y$’s and fit them to a parabola. While large $\delta\Phi_y$ values always behave quadratically, very close to $\Phi_y = \pi$ the changes in energy can be more erratic. In some cases, this may be set by the norm error or convergence of the DMRG, but a physical origin cannot be excluded. We therefore perform a parabolic fit excluding the points closest to $\Phi_y = \pi$. At finite $u_D$, the coefficient of the linear term in $\delta\Phi_y$ becomes non-zero and is fitted as well. We estimate our fitting error between different choices of which points to include to be on the 10% level at most, and usually on the 1% level. A well-fit point is shown in Figure 15.

Power Laws

The correlation function $C_O(x;\chi) = \langle\hat{O}(x)^\dagger\hat{O}(x)\rangle$ for an operator $\hat{O}$ with quasi-long range order manifests in DMRG as a power law $C(x) \sim x^{-\eta}$ in the limit $\chi\to\infty$. However, at each finite $\chi$, $C(x)$ must decay exponentially at large distances. To fit $\eta$ from finite $\chi$ data, we take advantage of the scaling collapse form

for some scaling function $g(u)$. Following e.g. [79] we may then define a loss function for $\eta$ that is minimized when the curves from different $\chi$ are closest. To that end, define

and put

where $N_\chi$ is the number of bond dimensions, $\bar{y}(u;\eta) = \frac{1}{N_\chi}\sum_\chi y_\chi(u;\eta)$, $\Delta u = u_1-u_0$ is range over which all curves are defined. As each $y_\chi(u;\eta)$ has data at different $u$’s, the data is linearly interpolated to perform the integral as a Riemann sum. We fit using only $x > 10L_M$ to remove the non-universal “head” of the curve. Figure 16 shows typical fits.