Graphene moiré superlattices show an abundance of correlated insulating, topological, and superconducting phases. Whereas the origins of strong correlations and nontrivial topology can be directly linked to flat bands, the nature of superconductivity remains enigmatic. Here we demonstrate that magic-angle devices made of twisted tri-, quadri-, and pentalayer graphene placed on monolayer tungsten diselenide exhibit flavor polarization and superconductivity. We also observe insulating states in the trilayer and quadrilayer arising at finite electric displacement fields. As the number of layers increases, superconductivity emerges over an enhanced filling-factor range and in pentalayer extends well beyond the filling of four electrons per moiré unit cell. Our results highlight the importance of the interplay between flat and more dispersive bands in extending superconducting regions in graphene moiré superlattices.
Rich phase diagrams of quantum electronic phases have been realized in many graphene superlattice structures, but robust superconductivity is so far exclusive to twisted bilayer graphene (TBG) [1,2] and twisted trilayer graphene (TTG) [3,4]. Strikingly, TTG exhibits greater electric-field tunability [3,4], Pauli limit violation [5], and more strongly coupled superconductivity [3,4] in comparison to TBG. While these observed differences may serve as clues in identifying the superconducting origin of these systems, our ability to identify the truly universal features is ultimately limited by the dearth of robust superconducting moiré materials, suggesting that further progress lies in the discovery of other moiré superconducting systems.
We investigate twisted graphene multilayers where each successive layer is twisted by an angle ±𝜃 relative to the previous one in an alternating sequence (Fig. 1A). For an even number 𝑛 of layers, the spectrum at zero displacement field 𝐷 is expected to separate into 𝑛/2 independent TBG-like bands, each characterized by a different effective twist angle. When the number of layers 𝑛 is odd, in addition to (𝑛 -1)/2 TBG-like bands, one monolayer-graphenelike (MLG-like) band (essentially a Dirac cone) is expected [6] (see left column of Fig. 1B for examples when 𝑛 is 3, 4 and 5). The system may be conveniently modified through the application of a displacement field 𝐷, which controllably hybridizes the different bands (Fig. 1B, right column). Experimentally, we explore the properties of alternating twisted trilayer, quadrilayer, and pentalayer graphene (TTG, TQG, TPG) structures with 𝜃 = 1.52 ± 0.02 (device D1, trilayer), 𝜃 = 1.80 ± 0.04 (D2, quadrilayer), and 𝜃 = 1.82 ± 0.05 (D3, pentalayer), respectively [7]. These twist angles all lie close to the theoretically predicted "magic"
values needed to obtain one set of flat TBG-like bands (𝜃 = √2𝜃 ≈ 1.53, 𝜃 = (√5 + 1)𝜃 /2 ≈ 1.75 , and 𝜃 = √3𝜃 ≈ 1.87 assuming an effective TBG twist angle 𝜃 = 1.08 [6]; see section 4.1 in [7]). We find that TTG, TQG, and TPG all exhibit hallmark signatures of strong correlations (Fig. 1, C to E), including robust superconductivity and flavor symmetry breaking as revealed by pronounced resistance peaks around certain integer filling factors 𝜈 (number of electrons per moiré site; see section 1 in [7] for assignment of 𝜈).
In addition to the symmetry-breaking transitions previously reported in TTG [3,4,8], our TTG structure (coupled to a tungsten diselenide (WSe ) monolayer [9]) exhibits a previously unobserved correlated insulating state near 𝜈 = +2 at finite 𝐷 (Fig. 1C, inset; see also Fig. S4.
1B right column; see sections 3 and 4 in [7]) and is instead attributed to the interplay between an interaction-driven cascade transition and hybridization induced by the 𝐷 field (e.g., as captured by Refs. [10,11]). We have also detected an insulating state developing at finite 𝐷 fields in TQG near charge neutrality (Fig. 1D, inset, and Fig. 1G). However, in contrast to TTG, the TQG insulating state can be explained through the 𝐷-induced hybridization only. Importantly, the detection of insulating gaps in TTG and TQG implies a low level of disorder in our samples (see also Fig. S1).
The superconducting regions in all three structures extend over significantly larger filling factor ranges in comparison to TBG [2,12,9,13] where superconductivity is typically observed within 2 < |𝜈| < 3. Moreover, superconductivity emerges over successively broader regions of phase space, reaching 𝜈 ≈ +5 on the electron side for TPG (Fig. 1, C to E). Along with a zero longitudinal resistance 𝑅 observed in the characteristic 𝜈 vs. 𝑇 dome (Fig. 1, H to J), we also measure large critical currents (∼ 400 nA) and occasionally see oscillations in critical current caused by superconducting interference, substantiating the robustness of phase coherence [7].
Moreover, we measure high critical perpendicular magnetic fields 𝐵 (typically ∼ 0.8 T; see Fig. S3), indicating that the corresponding Ginzburg-Landau coherence lengths 𝜉 (approximately 10 -30 nm) are significantly smaller than those observed in TBG and deviate from the weak-coupling prediction, 𝜉 ≈ ℏ𝑣 /𝜋Δ with Δ ≈ 1.76𝑘 𝑇 ; this suggests a strongcoupling origin of superconductivity [3,4] (see section 2 in [7]). When combined with other recent experiments [5,14,15], these observations affirm the unconventional nature of superconductivity within the entire class of graphene moiré systems. Further, the measurements on three to five layers indicate that the addition of layers promotes superconductivity over a broader filling window despite the coexisting dispersive bands as well as the ostensibly increased vulnerability to disorder-both from the additional twist angles as well as from the sensitivity to the relative displacement between layers (Fig. S14).
In addition to the pronounced 𝜈-dependence, the observed superconducting pockets are highly tunable with electric displacement field 𝐷 (Fig. 2). A comparison of the three structures reveals, however, that TQG and TPG are more tunable than TTG. This is apparent both in the 𝐷dependent evolution of the filling range where superconductivity is measured (Fig. 2, A to C) as well as in the critical temperature 𝑇 (Fig. 2, D to F). Notably, superconductivity in TQG and TPG is fully quenched for all fillings at 𝐷/𝜖 = 0.75 Vnm and 𝐷/𝜖 = 0.6 Vnm , respectively.
In the case of TTG, however, superconductivity is present up to the maximum accessible electric field 𝐷/𝜖 = 1 Vnm . Nevertheless, 𝑅 versus 𝐷 and temperature measurements do show that superconductivity is suppressed for large 𝐷 fields at optimal doping (the doping at which superconducting 𝑇 is the highest) in all three heterostructures; further, they reveal that 𝑇 forms a 𝐷 symmetric dome maximized at small finite 𝐷 (Fig. 2, D to F, for electron-side data showing similar behavior see Fig. S5). We also note that TTG, TQG, and TPG all exhibit a similar variation of 𝑇 when viewed as a function of the potential difference 𝑈 between the top and bottom layers (Fig. S5, D andE; see also section 3 in [7] for the energy conversion from 𝐷 to 𝑈).
This layer-number invariance is consistent with non-interacting continuum-model calculations tracking the evolution of the inverse of the flat-band bandwidth with 𝑈 (Fig. 2G, bottom). The dependence of 𝑇 on 𝐷 in all devices qualitatively matches the predictions of Ref. [16] for TTG with one marked exception: the observed vanishing of superconductivity and the decay of 𝑇 appears to be linear in 𝐷 (Fig. 2, E and F, and Fig. S5), in line with predictions for multilayer graphene with rhombohedral stacking [17] and in contrast to the exponential 'tail' typically expected from the weak-coupling theory (and seen in the model of Ref. [16]).
Comparing the location of the superconducting regions with the evolution of the Hall density as a function of 𝐷 and 𝜈 in TTG, TQG, and TPG provides further insight into the intricate relationship between the superconducting phase and the correlation-modified Fermi surface (Fig. 3). As in previous TBG and TTG measurements, we observe symmetry-breaking electronic transitions (a 'cascade' of transitions) that are signalled by sudden drops in the Hall density magnitude (a 'reset') without a change in sign. These resets (see dashed lines in Fig. 3 [20,21,22] and naturally disfavor Cooper pairing of time-reversed partners.
At high 𝐷 fields, signatures of the cascade vanish and instead van Hove singularities (vHs) become more prominent, reflecting qualitative changes in the band structure (see yellow lines in Fig. 3, A to D, and Fig. S6 that track the vHs). Consistent with previous TTG measurements [3,4], the vHs in our TTG sample (as well as in TPG, see Fig. 3, E and F) crudely bound the superconducting regions. By contrast, the vHs in TQG cross well into the superconducting pockets-in fact, for electron doping, 𝑇 reaches its maximum exactly at the position of the vHs (Fig. 3D, orange dot, and Fig. S7,D to F). The interplay between the vHs and superconducting boundaries, as revealed in Hall density measurements, is complex: 𝑇 can be both enhanced and suppressed at the vHs depending on the layer number and possibly other details such as the precise twist angle (see Fig. S11 for data from the second TQG device).
Pentalayer measurements provide additional signatures that point towards a close relation between superconductivity and flavor symmetry-breaking cascades (Fig. 3, E and F). In contrast to TTG, in TPG we can access 𝐷 fields that are large enough to stifle superconductivity-which occurs simultaneously with the onset of the vHs and the apparent suppression of the cascade transitions (see red and light blue lines in Fig. 3F that mark the superconducting boundaries and the cascade transitions, respectively). For example, at low 𝐷 fields (|𝐷|/𝜖 < 0.6 Vnm ) around 𝜈 = +2, the Hall density resets close to zero, in line with a nearly complete flavor symmetry-breaking polarization. However, at higher 𝐷 fields (|𝐷|/𝜖 > 0.6 Vnm ), the Hall density is dominated by a vHs around 𝜈 = +2, and the cascade signatures are diminished. Superconductivity accordingly also vanishes. For hole doping, the disappearance of superconductivity similarly coincides with the weakening of the cascade. This on/off correspondence between the two phenomena suggests that they either share a common origin, such as a large DOS, or that the cascade serves as a prerequisite for robust superconductivity in graphene moiré superlattices.
As mentioned above, for low 𝐷 fields in TPG, the superconducting pockets are extraordinarily large, spanning -4 ≲ 𝜈 < -2 for hole doping and +2 ≲ 𝜈 ≲ +5 for electron doping (Fig. 1E, Fig. 2C, and Fig. 4). In particular, the electron-side range corresponds roughly to a density window of 6 × 10 cm , which is the largest filling range so far reported in a graphene-based superconductor. The observed superconductivity exhibits similar values of 𝑇 and 𝐵 (Fig. S3) as the trilayer and quadrilayer samples and is likewise accompanied by a weak oscillating pattern of critical current (Fig. 4C, inset), confirming superconducting phase coherence. We emphasize that the unprecedented persistence of superconductivity across a large filling factor range in TPG (and also TQG in comparison to TTG or TBG) cannot be explained in a minimal framework of alternating twisted graphene multilayers [6,23] In scenario (𝑖) , 𝜈 corresponds to 𝜈 ≈ +3 , the flat-band filling at which superconductivity is typically suppressed in TBG, suggesting that the superconducting phase space is largely the same for different multilayer magic-angle structures when considering just the flattest TBG-like bands. In scenario (𝑖𝑖), 𝜈 coincides with 𝜈 ≈ +4, precluding any simple analogy with TBG, although superconductivity can still be attributed to the flat bands.
Finally, scenario (𝑖𝑖𝑖) assumes full filling of the flat bands before superconductivity is suppressed at 𝜈 . This scenario includes the possibility that the distinction between the different TBG-and MLG-like bands breaks down even at 𝐷 = 0 owing to hybridization as well as potential multi-band superconductivity [24,25,26,27]. In this case, superconductivity in TPG is a more general phenomenon than in TBG because it occurs in either mixed bands or new, more dispersive bands (see section 5 in [7]).
From the perspective of the non-interacting band structure, the three scenarios all seem implausible, therefore interactions must play a crucial role. In particular, although the presence of the dispersive bands implies that |𝜈| -|𝜈 | > 0, this effect is much smaller than needed for either scenario (𝑖) or (𝑖𝑖). Coulomb interactions can significantly enhance |𝜈| -|𝜈 |, either by evening out the spatial charge distribution [28,29,30,31,32] or symmetry breaking. A simple model for TPG incorporating these mechanisms (section 4.5 in [7]) suggests a minimal flat-band occupation 𝜈 ≳ +3.8 at 𝜈 ≈ +5, diminishing the plausibility of scenario (𝑖) for electrondoped TPG which has 𝜈 ≈ +5. The relevance of this scenario is further undermined by the observation of vHs at 𝜈 ≈ +6 (Fig. S10D): under the reasonable assumption that the noninteracting band structure remains valid for the dispersive TBG-like bands (apart from a Hartree shift), scenario (𝑖) would instead place the observed vHs near 𝜈 ≈ +5. Taken together, these arguments effectively rule out scenario (𝑖). Note, however, that the presented line of reasoning is not straightforward for the other superconducting pockets (see section 5 in [7]).
Both scenarios (𝑖𝑖) and (𝑖𝑖𝑖) are indicative of the non-trivial role of additional bands in stabilizing superconductivity. Assuming well-defined flat and dispersive bands, in scenario (𝑖𝑖𝑖)
the former bands are completely filled, and superconductivity is supported fully by the latter nonflat bands. This assertion is at odds with the large dispersion of the remaining TBG-and MLG-like bands. However, although the exact mechanism underlying scenario (𝑖𝑖𝑖) is difficult to pin down, it is not without experimental support. For instance, a natural interpretation of the Hall density minimum around 𝜈 ≈ +4 for |𝐷| ≲ 0.4 Vnm is that it marks the complete filling of the flat bands, 𝜈 ≈ +4 (Fig. 4E, Fig. S9 and section 5 in [7]).
One possible realization of scenario (𝑖𝑖𝑖) consistent with the experimental observations is that the division of the electronic states into simple TBG-and MLG-like bands fails-obviating our very definition of 𝜈 and potentially allowing flavor polarization, and accompanying superconductivity, to persist well beyond 𝜈 = +4. Whereas such hybridization is expected for finite 𝐷 fields, mixing between flat, dispersive TBG-and MLG-like bands for |𝜈| < |𝜈 | may occur even at 𝐷 = 0 thanks to, for example, proximity to WSe , layer-to-layer charge inhomogeneity or distant-layer coupling (see sections 4.4 and 4.5.3 in [7]).
Our measurements demonstrate the increasing predominance of superconductivity in twisted graphene multilayer structures as the number of layers is increased from three to five and highlight the close relationship between the flavor symmetry-breaking transitions and superconductivity. Moreover, our findings suggest a scenario in which the symmetry-broken 𝜈 = ±2 state strongly favors the formation of the superconducting state whereas the cascade corresponding to 𝜈 = ±3 suppresses it. Interestingly, this scenario is consistent not only with previous TBG [2,12,13,9,33,34] and TTG [3,4,8] observations but also in part with the recently investigated ABC trilayers [35] and Bernal bilayers without [36] and with a WSe substrate [37] where superconductivity is observed near or within phases in which two out of four flavors are predominantly filled. These commonalities suggest that symmetry-broken states with similar types of polarization underlie superconductivity in all these graphene-based superconductors. In this context, the discovery of superconductivity in TQG and TPG together with recent work on untwisted bi-and trilayers dramatically expands the scope of graphene-based superconductors.
This expansion holds promise for resolving important questions related to the nature of the pairing mechanism in these systems and provides guidance for developing novel graphene-based superconductors and their applications.
Note added: After submission of this report, we have become aware of related work [38].
Device fabrication: All devices were fabricated using a 'cut and stack' method, in which graphene flakes were separated into pieces using a sharp tip (made out Platinum-Iridium); this approach prevents unwanted twisting and strain during tearing while allowing more control over the flake size and shape. After cutting, stacking procedure was as follows: first, a thin hBN flake (10 -30 nm) is picked up using a propylene carbonate (PC) film previously placed on a polydimethylsiloxane (PDMS) stamp. Then the hBN flake is used to pick up an exfoliated monolayer of WSe 2 (commercial source, HQ graphene) before approaching the graphene. After picking up the first piece of the graphene flake, the following layers are twisted by an angle ±θ relative to the previous one in an alternating sequence. Transfer stage rotation θ overshoots the target angle by 0.1 -0.2°to construct the measured angles. Care was taken to approach and pick up each stacking step slowly. In the last step, a thicker hBN (30 -70 nm) is picked up, and the whole stack is dropped on a predefined local gold back gate at 150°C while the PC is released at 170°C. The PC is then cleaned off with N-Methyl-2-Pyrrolidinone (NMP). The final geometry is defined by dry etching with a CHF 3 /O 2 plasma and deposition of ohmic edge contacts (Ti/Au, 5 nm/100 nm) and top gate.
Measurements: All measurements were performed in a dilution refrigerator (Oxford Triton) with a base temperature of ∼ 25 mK, using standard low-frequency lock-in amplifier techniques. Unless otherwise specified, measurements are taken at the base temperature. Frequencies of the lock-in amplifiers (Stanford Research, models 830 and 865a) were kept in the range of 7 -20 Hz in order to measure the device's DC properties and the AC excitation was kept < 5 nA (most measurements were taken at 0.5 -1 nA to preserve the linearity of the system and avoid disturbing the fragile states at low temperatures). Each of the DC fridge lines pass through cold filters, including 4 Pi filters that filter out a range from ∼ 80 MHz to > 10 GHz, as well as a two-pole RC low-pass filter.
1 Device Uniformity, Effect of WSe 2 and Twist Angle Assignment Device homogeneity and effect of WSe 2 : All devices investigated here show a high degree of twist angle homogeneity as characterized by four-point measurements between different pairs of contacts. Fig. S1 shows R xx versus carrier density with fixed top-gate voltage (V tg = 0 V), revealing that almost every pair of contacts shows superconductivity. More importantly, superconducting pockets from different pairs significantly overlap in the filling range, and resistance peaks at |ν| = 4 appear at the same density. Moreover, all findings related to the extent of the superconducting phase and the occurrence of the symmetry-breaking transitions in the ν-D phase diagram are highly reproducible. This also includes the observation of a gapped correlated insulator at ν = +2 in TTG, which has not been reported previously. In this context, we note that any significant twist-angle disorder would create conducting percolation pathways that quickly suppress insulating behavior.
We attribute the low level of disorder to the use of monolayer WSe 2 during device stacking, presumably originating from the increased lateral friction between WSe 2 and graphene (compared to the friction at the hBN-graphene interface). We note that this additional layer does not change the magic-angle condition (9,40), and the induced spin-orbit interaction (SOI) energy scale is ∼ 1 meV in twisted bilayers (9). Therefore, SOI is likely too small to significantly affect the overall band structure and directly impact the cascade physics at the magic angle (though may play a more important role for stabilizing superconductivity far away from the magic angle ( 9)).
Finally, we note that, in general, SOI is expected to manifest differently when the sign of D field is reversed, a feature that has not been observed in the experiment. The absence of D-field asymmetry is probably due to the small energy scale of SOI compared to the interactions and the weak tendency to polarize wavefunctions with D field in magic-angle graphene multilayers.
Twist angle assignment in multilayers: Twist angles were determined from high B field data and corresponding Landau-fan diagrams in a similar way as in TBG. From the slope of the Landau fan at charge neutrality (which is directly proportional to the gate-sample capacitance) and the voltage difference between charge-neutrality point (CNP) and |ν| = 4 filling, the corresponding |ν| = 4 electron density is obtained. We used two separate criteria for the assignment of |ν| = 4. First, at high D fields, resistive features (peaks) emerge (Fig. 2, A to C). We interpret these peaks presumably as the opening of the hybridization 'gaps' and corresponding full filling (|ν| = 4) of the 'gapped' bands. Second, at high B fields, quantum Hall insulating states develop around |ν| = 4, which typically cover a broader filling range where Hall conductance is quantized in accordance with the expectations from the dispersive bands (Fig. S2, also see discussion below). Electron density of |ν| = 4 directly determines the twist angle in the lowangle approximation θ 2 ≈ √ 3a 2 n |ν|=4 /8, where a = 0.246 nm is the graphene lattice constant. While twist angle uncertainty is introduced by the relatively broad experimental features in Landau fan diagram, this uncertainty can not explain the observed superconductivity reaching ν = +5 in the case of TPG. We note that signatures of the dispersive bands are also observed in Landau-fan diagrams and vHs (Fig. S2 and Fig. S10). For example, emerging Landau levels from the dispersive bands are typically observed through oscillations at low magnetic field. Since at low energies, the dispersive bands (Fig. 1B left) can be effectively treated as decoupled MLG-like bands when considering the Landau level spectrum (41), and the corresponding Hall conductance around |ν| = 4 will be quantized in a way that depends on the number of layers (see Fig. S2, G to M, for Hall conductance line cuts). Furthermore, both TQG and TPG devices show vHs signatures at |ν| ≈ 6, consistent with the expectation for vHs from roughly half filling of the dispersive TBG-like bands (see Fig. S10).
The above discussions implicitly assume that during fabrication process, errors in twist angles of successive layers are not large. In the following, we argue that this assumption holds in our devices to a large degree. When multiple, distinct twist angles are present (42), the resulting band structure is expected to become very complex (notably, no sense of single moiré periodic-ity remains) and deviates significantly from the observed experimental results. Our calculations have suggested similar complexity with varying stacking arrangements. Importantly, the resulting band structures noticeably deviate from the ones described in Fig. 1B, where for n layers at zero D field, there are n/2 (for even n) or (n -1)/2 (for odd n) sets of overlapping TBG-like bands (one flat) and 1 Dirac cone if n is odd. Although the parameter space is too large to provide a thorough analysis, we provide some references with varying geometries that do not agree with our experimental observations (43)(44)(45). Finally it is important to note that, under the condition that the twist angles don't vary much between successive stacking steps, previous TTG studies have shown that the lattice arrangement is energetically favorable in the A-twist-A stacking configuration (where the A sublattice from the first layer aligns with the third layer instead of being offset) and uniform twist angle between top and bottom layers (15,46).
T c and the coherence length: T c is determined by the following procedures. First, the high temperature R xx data is fitted using a linear function R(T ) = AT +B. Then, T c is defined by the value where R xx (T ) is a certain fraction (typically 10% as in Fig. 2) of R(T ). Ginzburg-Landau coherence lengths ξ GL are obtained from the B dependence of T c , by fitting the Ginzburg-Landau relation
is the superconducting flux quantum and T c0 is the critical temperature at zero magnetic field. We get ξ GL from the T c vs. B linear fit, where the intercept at the B axis is equal to Φ 0 /(2πξ 2 GL ). Following Ref.
(3), we use T c defined by 40% of the normal state resistance to evaluate the coherence length data in Fig. S3E (corresponding error bars are evaluated by using T c defined by 30% and 50% of the normal state resistance). As mentioned in the main text, ξ GL (B c ) is much smaller (higher) in the twisted graphene multilayers compared to TBG. One possibility for the reduction of ξ GL is the relative decrease of the characteristic moiré wavelength (see Fig. S3F).
Hall density analysis: Hall density shown in Fig. 3 is obtained by converting the anti-symmetric part of the R xy data, i.e., by subtracting data measured at positive and negative magnetic fields. We used either |B| = 0.9 T or 1.5 T in order to fully suppress superconductivity. Previously, it was found that in TTG (3), at high D fields superconductivity is bounded by regions corresponding to vHs, i.e., when Hall density changes sign. We approximately find a similar trend in our TTG and TPG structures, although vHs occasionally intrude superconducting pockets slightly. We note that the exact positions of vHs depend on the precise magnetic field used in the measurements (for example, see Fig. S9, A andC); however, this effect is relatively small relative to the observed intrusions. More importantly, TQG behavior is qualitatively different, as we find that positions of vHs and boundaries of superconducting pockets are independent. We also note that resets associated with the flavor polarization do not move in the B fields (B ≈ 1T) used to extract Hall density evolution. The occasional shift of these resets from integer ν values, may be attributed to either effects of interactions (i.e. Hartree correction, see section 4.5) or the details of cascade physics (18) at finite temperatures (47,48).
3 Insulating Behavior in TTG and TQG ν = +2 correlated insulators in TTG: In our TTG device, we observed a previously unidentified ν = +2 correlated insulator state (Fig. 1F and Fig. 1C inset). This assignment is based on the activating behavior reflected in exponential suppression of conductance with 1/T , which is not expected from a Dirac semimetal; see Fig. S4 for detailed D and ν dependence. The two-temperature scale in conductance is in line with the behavior of TBG correlated insulator where transport at higher temperature and lower temperature is governed by activation and variable-range hopping, respectively (see blue dashed line in Fig. 1C inset). Also, the insulating behavior is suppressed by an out-of-plane B field of B ⊥ > 0.4 T but is mostly insensitive to the in-plane B field (the weak suppression by the in-plane B field could possibly originate from a non-ideal in-plane-field alignment of the sample). These experimental observations are highly indicative of a gap that originates from strong interactions in TTG. We note, however, that formation of the fully gapped states in TTG requires a mechanism that additionally gaps out the MLG-like band, which may explain the presence of the gap only at finite D fields. Moreover, suppression of the gap with an out-of-plane magnetic field is at odds with the C 2 breaking scenario (49,50) and is more in line with incommensurate Kekulé spiral (51) or intervalley-coherent (10,11,52,53) orders in the flat bands. The insensitivity to in-plane field is suggestive of a spin-polarized insulator or otherwise insensitive to in-plane magnetic field. Finally, we can not rule out that the gap originates from induced SOI, since it is still possible that SOI promotes instabilities that favor the formation for certain ν = +2 insulating states in TTG. Future work will address the nature of this state in more detail.
Charge-neutrality gaps in TQG and conversion between U and D: Fig. S4, F andG, shows the charge-neutrality gap of TQG as a function of D field or potential difference U (between the top and the bottom graphene layer). From the continuum model, a gap in TQG is expected when finite D field is applied. However, the details of the gap evolution depend on the precise twist angle. When the twist angle is below the magic-angle value, a charge-neutrality gap opens as soon as a finite D field is applied. On the other hand, when the twist angle is above the magicangle value, a gap opens only at much higher D fields. The gap opening at D/ϵ 0 ≈ 1.1 V nm -1 in our TQG structure is consistent with the device being slightly above the magic angle. Note that the charge-neutrality gap is a good reference for matching the experimental D field with the potential difference U used in calculations since the interaction-driven Hartree correction vanishes at CNP. For a direct comparison, we enforce a more realistic flat-band bandwidth of ∼ 20 meV in the continuum model by slightly tuning away from the magic angle, and get a Udependent gap size (Fig. S4G). A good match between the experimental and the calculated gap is found when converting D into U with an empirical factor: U = 0.1×(n-1)×0.33 nm×eD, where e is the electron charge and n -1 is the number of graphene interfaces. This conversion is used for the other parts of the paper, for example, T c versus U in Fig. S5. We note, however, that relative comparison (i.e. scaling) between TTG, TQG, and TPG (in the context of T c ) does not rely on the precise D to U conversion.
In this section, we describe the non-interacting continuum model for multilayers and how symmetry considerations and various interaction terms affect the band structure of TTG, TQG, and TPG.
Band structure calculations are performed using a straightforward generalization of the TBG continuum model (54,55) extended to multilayer graphene systems (6,(56)(57)(58). As discussed above, we consider graphene multilayer systems with n layer = 3, 4, and 5 layers (n in the main text) in which the graphene sheets are twisted by alternating angles. In particular, we can envision grouping the layers into even and odd sets and then rigidly twisting these two groups by the twist angle θ; equivalently, each layer ℓ = 1, . . . , n layer is twisted by an angle θ ℓ = (-1) ℓ θ/2. For the moment, we focus on the case where the layers are all AA stack (i.e. stacked directly on top of one another) prior to twisting (see below, section 4.2).
It is appropriate to approximate the dispersion of the underlying graphene monolayers with the two Dirac cones about the valleys at K and K ′ . Note that because of the twist, the Dirac cones are located at slightly different momenta depending on whether the layer ℓ is even or odd, and we have denoted the Dirac cones' momenta here as K ℓ and K ′ ℓ . We thus define the spinors ψ ℓ,K (′) in terms of the microscopic graphene operators via f ℓ (r) = e iK ℓ •r ψ ℓ,K (r) + e iK ′ ℓ •r ψ ℓ,K ′ (r). Equivalently, in momentum space, we can write
ℓ . In our definition of ψ ℓ (and f ℓ ) both an A/B sublattice index and a spin index have been suppressed. Importantly, the small twist angle only mediates a very small momentum exchange between the neighboring layers so that states originating proximate to one valley do not mix with those originating proximate to the other. We thus focus for the moment on valley K and suppress the valley subscript until mentioned otherwise,
The band structure model can be separated into a sum of two parts: H cont = H D + H tun . The first term, H D is the intralayer Dirac term:
Here, v 0 ∼ 10 6 m/s is the Fermi velocity of the Dirac cones and σ x,y,z are Pauli matrices acting on the A/B sublattice indices of the spinors ψ ℓ . In our simulations, we assume that the Fermi velocity of the graphene monolayers does not differ layer to layer. Note that this assumption, specifically does not take into account effects such as graphene velocity renormalization that can occur in the top layer due to tunnelling between the graphene monolayer and the WSe 2 substrate, as we do not expect these effects to be large enough to have appreciable impact on the resulting band structure.
We assume that tunnelling only occurs between adjacent layers and that it takes the form
where
e -(-1) ℓ iq j •r t j ,
Here, R(ϕ) = e -iϕσ y is a 2 × 2 rotation matrix acting on vector indices, L M = a/[2 sin(θ/2)] is the moiré lattice constant, and σ ± = (σ x ± iσ y )/2 act on the sublattice indices. The parameters w ′ and w set the interlayer tunnelling strength; we discuss their values below. It will be convenient to define the dimensionless ratios
where k θ = 4π/(3L M ) = 2 sin(θ/2) • 4π/(3a). The total Hamiltonian may be written in matrix form as
As currently written, the diagonal Dirac terms, h D,ℓ (r), as well as the off-diagonal tunnelling terms, T ℓ,ℓ ′ (r), depend only on whether ℓ is even or odd. We can thus simplify the above expression by writing the Dirac terms as h D,2ℓ-1 (r) = h D,1 (r), h D,2ℓ (r) = h D,2 (r) and the tunnelling terms as T 2ℓ-1,2ℓ (r) = T (r), T 2ℓ,2ℓ+1 (r) = T † (r).
It has been shown ( 6) that a block diagonal form exists for Hamiltonians of the form Eq. (S5). We provide the specific transformations used for three, four, and five layers below.
Since the spectrum of the twisted multilayers breaks into independent sets of TBG-and MLGlike bands, we first briefly review the Hamiltonian of TBG. Thus, we start with
Provided that inversion and time reversal symmetries are preserved, the Dirac cones described by h D,ℓ (r) at K ℓ are preserved even when the interlayer tunnelling is added. Nevertheless, this tunnelling term breaks the (effective) continuous translation symmetry of h D,ℓ . Consequently, the set of conserved momenta are confined to reduced moiré Brillouin zone (BZ). Like the original monolayer graphene BZ, the moiré BZ forms a hexagon with the Dirac cones located at the corners. Here, we define K 1 = κ and K 2 = κ ′ (for the other valley,
For small twist angles, the intralayer Dirac terms are nearly identical, h D,ℓ (r) ≈ h D (r) = v 0 (i∂ x σ x + i∂ y σ y )-namely, the rotation in Eq. (S1) may be neglected to first order. In this case, the spectrum depends solely on the ratios η = w ′ /w and α = w/(ℏv 0 k θ ), where k θ = 4π/(3L M ) is the distance separating κ and κ ′ . Further, as shown in Ref. (55), the spectrum close to Dirac points at κ and κ ′ can be approximated using a simple perturbative scheme. In particular, in momentum space one finds
The magic angle is defined (55) by the condition v α,η = 0, which we see here should occur for α ≈ 1/ √ 3.
The Hamiltonian for the three layer system is
It maybe be transformed into a block diagonal form as
First, we note that the TTG spectrum separated into independent sets of bands-a TBG-like set described by the two-layer Hamiltonian h √ 2α,η,θ (an 8 × 8 object when sublattice and spin are included) and an MLG-like Dirac cone described by h D,1 . Using the reasoning above, we expect a set of flat TBG-like bands to occur when √ 2α = 1/ √ 3. Equivalently, we can assign an effective TBG twist angle describing these bands, θ eff TBG = θ/ √ 2 where θ is the physical twist angle of the system. If θ eff TBG is expected to yield flat bands in TBG, then we would similarly expect √ 2θ eff TBG to yield a set of flat bands in TTG.
For four layers, we start with
With the appropriate change of basis, we obtain
where φ = (1 + √ 5)/2 is the golden ratio. In this case, we therefore expect the TQG spectrum to possess two sets of TBG-like bands characterized by effective TBG twist angles θ/φ and θ/φ -1 .
The final system considered is the twisted pentalayer graphene. In the original layer basis, the Hamiltonian is
Once more, independent, co-existing TBG-and MLG-like subsystems are revealed with the appropriate change of basis:
There are now two independent TBG-like bands characterized by effective twist angles θ/ √ 3 and θ in addition to a MLG-like Dirac cone.
As indicated in Eq. (S7), the magic-angle value is essentially determined by the velocity of monolayer graphene v 0 and the interlayer tunnelling amplitude w. We fix v 0 for all considered configurations. The magnitude of the interlayer tunnelling amplitude is typically estimated to be around ∼ 100 meV. In case of TQG, a gap is expected to open at charge neutrality when finite D field is applied. However, it onsets for any |D| > 0 when the physical angle θ is less than φθ magic TBG , where θ magic TBG is the magic angle for TBG (as determined by v 0 and w). When θ is larger than φθ magic TBG , a gap still opens, but only above certain finite D fields. As Fig. S4 and Fig. S7J show, the latter scenario is observed in the TQG, device D2 (twist angle 1.8°), leading us to select w = 108 meV near the value used in Ref. (55). In particular, in the left panel of Fig. S7J, the R xx is plotted as a function of D field, displaying non-monotonic behavior-a resistance dip around D/ϵ 0 ∼ 0.5 V nm -1 followed by a steep increase at higher D, signalling the development of an insulating gap. Analogous trends are repeated on the right of Fig. S7J, which shows an increase in ν = 0 DOS (corresponding to the resistance dip) followed by a decrease to zero DOS. Similar reasoning can be applied to the TTG and TPG samples, although it is slightly more nebulous since a non-interacting gap is not expected to open in TTG and TPG for any D value at the CNP. Instead, when θ > √ 2θ magic TBG for TTG and θ > √ 3θ magic TBG for TPG, the system should become metallic with increasing D, whereas in the converse situation, the D field should immediately gap out all states except for a dispersive MLG-like Dirac cone. Following this line of reasoning, the results of Fig. S7 suggest that the twist angle in TTG is below the magic angle, whereas the one in the TPG sample is above the magic angle. Accordingly, we select w = 110 meV for TTG and w = 102 meV for TPG modeling. Notably, the resistance behavior and theoretical DOS shown in Fig. S7K for TPG are very similar to the results in Fig. S7J with the primary distinction being that the high-displacement field state does not display the activated transport of an insulator. Similarly, although not obvious from the DOS plot itself, the band structure of TPG at large D is semimetallic (as opposed to insulating).
The value of the interlayer AA hopping w ′ is expected to be less than the interlayer AB hopping, as a result of lattice relaxation (see next section). We chose w ′ = 60 meV for all three multilayers considered, which is in agreement with the estimates of Ref. (23) and similar to values used previously for TBG/WSe 2 structures (9). We note that, while consistent with experiment in the fashion outlined about, other factors could be also at play, modifying the behavior at CNP in ways not captured by our analysis. Ultimately, however, the choices made here are not expected to greatly influence any of the results in this section.
An important distinction between TBG and graphene moiré heterostructures containing additional layers is the band structure dependence of the relative layer displacement. Not only must the graphene sheets be stacked with alternating angles, as discussed in the main text and in the previous section, but moreover, the emergence of independent TBG-and Dirac-like bands only occurs when all odd (even) layers are AA stacked, i.e., stacked directly on top of one another. As stated above, we envision grouping the layers into odd and even sets, each stacked rigidly atop one another. The two sets are then twisted relative to one another by the twist angle θ. We have assumed that this stacking was realized in the previous presentation and now argue for the feasibility of this assumption.
In TTG, it has been numerically shown that this situation is energetically preferable: the system naturally relaxes into the odd/even aligned stacking configuration (56). This result is further experimentally verified in transport (3) and local probe (15) measurements. A simple heuristic supports these results and permits a generalization to additional layers. Starting from a bilayer system, the moiré superlattice is manifest on the microscopic lattice scale as the periodic variation of the relative interlayer stacking: one has AA regions at the moiré hexagon centres, while AB and BA stacking regions represent the moiré hexagon vertices. The AA regions have a relatively high energy compared to the Bernal-like region and the lattice accordingly responds by relaxing to minimize their area. We now consider adding a third layer with the same relative twist angle as the first layer, but for the moment arbitrarily displaced from that layer. A moiré superlattice is of course also generated between the new layer and the second layer, and the system once again seeks to minimize (maximize) the area of the AA (AB/BA) regions. Crucially, if the first and third layers are misaligned, the AA regions between the first and second layers are misaligned from the AA regions between the second and third layers, frustrating the ability of the lattice to relax. Only when the first and third layers are aligned will the AA region occur at the same locations and only then can the system optimize its energy through relaxation. These arguments clearly generalize to quadrilayer and pentalayer systems-we need only consider the moiré pattern generated by each adjacent pair of graphene sheets to conclude that relaxation is optimized by an odd/even aligned configuration. (A complementary explanation is that TTG is necessarily an intermediate step in the construction of the TQG and TPG devices, and thus the Here, ℓ, ℓ ′ again indicate the system layers and σ x acts on the sublattice indices. The fact that this symmetry is not onsite means that we cannot generically partition the system into even and odd parity sectors like we did for TTG and TPG. For TQG, hybridization between subsystems is therefore not prohibited by symmetry.
In obtaining the independent TBG-and MLG-like bands (subsystems) listed above, a number of assumptions were made and one may be concerned about the relative robustness of these results. For TTG, at least, this question may be dismissed so long as mirror symmetry is present; above, we showed that this mirror symmetry protects the block diagonal subsystem form obtained for TTG. Similarly, mirror symmetry disallows mixing in TPG between certain (but not all) subsystems. However, the mirror symmetry is explicitly broken by the application of a displacement field as well as by the WSe 2 substrate used in the experiment. Below, we show that these modifications induce mixing between all subsystems. We additionally consider other mirror-preserving effects that may result in subsystem mixing in TQG and TPG.
Besides the displacement field, we find that the subsystem-mixing energy scales discussed below are relatively small compared to the input parameters of the continuum model, i.e., compared to an interlayer tunnelling of w and w ′ . More importantly, they are also smaller than the observed bandwidth of TBG, which spectroscopic measurements indicate is ∼ 40 meV for samples close to the magic angle (59)(60)(61)(62). The subleading magnitude of the effects we explore below thus bolsters our use of the alternating-angle continuum model, at least as a starting point. We note that the relatively small subsystem hybridization discussed here could be significantly magnified by interactions.
In the main text, we allude to the fact that a finite displacement field mixes the TBG-and MLGlike subsystems obtained in the previous sections. This effect is included in the Hamiltonian through the addition of
Specifically, we have h
disp = U diag(1/2, 1/6, -1/6, -1/2), and h
(5) disp = U diag(1/2, 1/4, 0, -1/4, -1/2). Here, the scale U is defined as outlined in section 3.
One of the exterior layers of the samples considered here is placed adjacent to WSe 2 . This type of construction was first shown to induce spin-orbit coupling in twisted bilayer graphene in Ref. (9). The presence of WSe 2 breaks not only the spin symmetry, but also the mirror symmetry in systems considered here and possibly induces subsystem mixing. Similarly, the invariance of the spectrum under U → -U (D → -D) will also no longer strictly hold. Nevertheless, the magnitude of the induced spin-orbit scale has been measured to be approximately 1 -5 meV in TBG, smaller than the other scales of the theory (e.g., the bandwidth). In effect, in rotating to the subsystem basis, the spin-orbit terms are "spread" across an increasing number of layers by the unitary transformations V TTG , V TQG , V TPG .
The chemical potentials of the different layers may also differ in a way that is symmetric under onsite mirror actions U (4/5) mirror of Eq. (S15). In particular, we may have h (4) µ-var = diag(δµ 1 , δµ 2 , δµ 2 , δµ 1 ), which takes the subsystem-basis form
where δµ = δµ 1 = -δµ 2 . Similarly, for TPG, a term like h
µ-var = diag(δµ 1 , δµ 2 , δµ 3 , δµ 2 , δµ 1 ) also preserves the mirror operator U (5) mirror but can be shown to induce inter-subsystem mixing within the even parity sector. As we demonstrate in section 4.5, such a term is naturally generated by the Coulomb interaction. We specify to TPG in that section, but the reasoning is analogous for TQG (and for TTG, although this term will not induce mixing between sectors because of the mirror symmetry).
Although generically present, the Coulomb interaction-generated terms of this form are relatively small compared to the other terms present. The calculations presented below estimate that values of |δµ ℓ | < 5 -10 meV for TPG are generated as one dopes away from charge neutrality. We expect the results for TQG to follow the same trend.
Our Hamiltonian so far only includes tunnelling between neighboring layers. Generically, however, hopping between next-nearest neighboring layers occurs as well. For TQG, we could therefore consider hopping between layers 1 (2) and 4 (3):
Assuming for simplicity that T nnl. = T † nnl , in the subsystem basis, this term takes the form
Note that because we assume next-nearest layers are stacked AA relative to one another, to first order, no spatial dependence is expected in T nnl . The subsystems are similarly mixid with the five-layer analogue h
nnl . Reference ( 56) computed the values of T nnl expected in TTG (where it will not induce subsystem mixing) and found that a typical scale |T nnl,ij | ∼ 5 -10 meV, which translates to |T nnl,ij |/ √ 5 ∼ 3 -5 meV (i and j are sublattice indices).
As mentioned in section 4.2, the moiré lattice relaxes in order to minimize AA regions and maximize AB/BA regions. This relaxation effect ultimately depresses the value of w ′ (interlayer AA tunnelling) relative to w (interlayer AB/BA tunnelling) as a result of out-of-plane corrugation. For interior layers, which neighbor more than a single sheet, the effects of relaxation are naturally stronger than for exterior layers. Consequently, the value of η = w ′ /w appropriate for tunnelling to and from interior layers is reduced. Our assumption below Eq. (S5) that T ℓ,ℓ+1 (r) depended only on whether ℓ was even or odd is no longer valid. Unsurprisingly, this effect once again mixes the subsystems in TQG and TPG. Reference ( 23) estimated the magnitude of this effect and determined that it should be in the range 5 -10 meV for the twist angles considered here.
The presence of flat-band subsystem in the low-energy theory of the multilayer graphene structures necessitates the consideration of interaction-driven band structure corrections. In the following, we focus specifically on the case of TPG as its phase diagram demonstrates the strongest deviation from the minimal paradigm that a multilayer structure maybe thought of as a TBGlike Hamiltonian with spectating additional bands. Rather, as we argue, the dispersive TBG-and MLG-like subsystems play a crucial role in extending the filling range of the superconducting pocket in accordance with scenario (ii) and (iii). Here, we consider three types of interaction corrections: (a) an in-plane Hartree correction; (b) a two-parameter effective model mimicking generic Hartree-Fock modifications of band structure; (c) an out-of-plane Hartree correction allowing for inhomogeneous charge distribution between the layers. We demonstrate that these effects generically lead to two consequences for the electronic spectrum: promoting charge redistribution to the non-flat bands and also leading to possible symmetry breaking between the non-flat and flat bands.
We begin with an in-plane Hartree effect. As demonstrated experimentally in previous work on TBG (32) and TTG ( 15), filling-dependent interaction effects, specifically Hartree and Fock corrections, drastically alter the electronic dispersion. Here we incorporate only a Hartree mechanism (28)(29)(30)(31) in the analysis, arguing that its key effect, a relative shift of flat bands up in energy with respect to the non-flat bands, is the simplest mechanism through which the size of the superconducting pocket in TPG is extended. We then supplement this discussion with a phenomenological Hartree-Fock-like theory. Before proceeding, we stress that the main purpose of the analysis in this section is to demonstrate that scenario (i) wherein flat bands are filled only to ν flat ≈ +3 at ν ≈ +5 is highly unlikely, thus highlighting the non-trivial role played by the dispersive TBG-and MLG-like bands.
The foundations of the Hartree calculation in TPG described below are identical to the analysis in Refs. (32) and (15). We reproduce them here for the convenience of the reader. We introduce the Coulomb interaction into the system through
In section 4.1, we introduced creation and annihilation operators, ψ † (r) and ψ(r), that correspond to the non-interacting eigenstates given by the Hamiltonian of Eq. (S5). Here and in what follows, we suppress the layer, valley, sublattice and spin subscripts. In Eq. (S24), V (r) = e 2 /(ϵ|r|) is the Coulomb potential and δρ(r) = ψ † (r)ψ(r) -ρ CN (r), where ρ CN (r) = ⟨ψ † (r)ψ(r)⟩ CN is the expectation value of the density at the charge-neutrality point. The use of δρ(r) instead of ρ(r) in the interaction is motivated by the expectation that the input parameters of the model H Tn layer G = H cont already include the effect of interactions at the charge-neutrality point. The dielectric constant ϵ in the definition of V (r) is used as a fitting parameter; see discussion below for details.
We study the effect of the interacting continuum model of magic-angle TPG through a selfconsistent Hartree mean-field calculation. Instead of solving the many-body problem, we obtain the quadratic Hamiltonian that best approximates the full model when only the symmetric contributions of H C are included, i.e., the Fock term is neglected. Thus instead of H cont + H C , we study the Hamiltonian
H is the Hartree term at filling ν, H
and the last term in Eq. (S25) simply ensures there is no double counting when one calculates the total energy. In the above equation, V (q) = 2πe 2 /(ϵ|q|) is the Fourier transform of the Coulomb interaction V (r) in Eq. (S24), and the expectation value ⟨ Ô⟩ ν = ⟨ Ô⟩ occ -⟨ Ô⟩ CN only includes states that are filled up to ν relative to charge neutrality, as defined by diagonalizing the Hamiltonian H
MF . Typically, for a "jellium"-like model, the expectation value in Eq. (S26) vanishes save for q = 0, which is subsequently cancelled by the background charge-allowing one to set V (q = 0) = 0 and completely ignore the Hartree interaction. However, because the moiré pattern breaks continuous translation symmetry, momentum is only conserved modulo a reciprocal lattice vector. We therefore obtain
where the prime above the summation over the moiré reciprocal lattice vectors indicates that G = 0 is excluded. The self-consistent procedure begins by assuming some initial value of H
H and diagonalizing the corresponding mean-field Hamiltonian H
MF to obtain the Bloch wavefunctions and energy eigenvalues. These quantities are then used to re-compute the expectation values that define H (ν) H and thus H (ν) MF subject to the cascade treatment described above. This process is repeated until one obtains the quadratic Hamiltonian H (ν) MF that yields the correlation functions ⟨•⟩ ν used in its definition.
Due to the accumulation of electronic density at the AA sites of the stacking sequence, the Hartree potential is dominated by the first 'star' of moiré reciprocal lattice vectors (28,63), which in our conventions corresponds to
(1, 0) T for n = 1, . . . , 6, with R(ϕ) a rotation matrix. The restriction to the G n 's paired with the 2π/6 rotation symmetry of the continuum model greatly simplifies the calculation of the Hartree term. Notably, V (G) k ′ ⟨ψ † (k ′ + G)ψ(k ′ )⟩ ν must be the same for all G n , and, instead of Eq. (S27), we use
The self-consistent procedure in this case is identical to that described in the previous paragraph, but due to the reduced number of reciprocal lattice vectors that are included in the summation, the calculation is computationally easier. Convergence is typically reached within ∼ 6 iterations.
We now proceed to discuss the precise effect of the Hartree correction. Since the Hartree correction couples bare graphene states at momenta k and k + G, its effect is most pronounced for subsystems of the Hamiltonian whose eigenstates require multiple bare graphene states originating from multiple moiré BZs, e.g. k + G states with G extending beyond the second BZ. As such, Hartree affects the flat-band subsystem most severely since its eigenstates deviate the most from the bare graphene states, while the MLG-like subsystem is affected the least. As a result, this correction gives rise to an energy offset that shifts the flat bands upwards in energy with respect to the rest of the energy spectrum (technically the dispersive TBG-like subsystem is also shifted slightly with respect to the MLG-like subsystem). This effect has been seen both theoretically and experimentally in TTG (15,64) Thus we expect it to be present in TPG, as is confirmed through our simulations; see Fig. S12, A andB. Physically, this effect arises simply because the charge distribution from the non-flat subsystems is more homogeneous in the unit cell and, therefore, it contributes less to the potential of Eq. (S26).
We now discuss what happens when one starts from charge neutrality and electron dopes the system. Due to the shift of the flat band upwards in energy relative to the non-flat bands, more charge can enter the non-flat bands upon doping (increasing ν) than a naïve non-interacting model predicts. As a result, the filling range of the flat TBG-like band superconducting pocket may be extended since the filling of the flat bands ν flat can continue to lie in the range amenable to superconductivity, whilst the total filling ν increases by adding charge to the non-flat bands. This is the central idea behind the scenarios (i) and (ii) discussed in section 5.
In the simulations for flat-band filling ν flat > +2, we consider two ways to fill the otherwise 4-fold degenerate bands: an uncascaded model where all 4-fold degenerate bands are filled equally, and a simple cascade model where two of the flat-band flavors (say spin ↑ for K, K ′ ) are shifted down in energy such that the highest energy of the shifted bands falls on the Dirac point of the unshifted flat bands, c.f. Fig. S12C. In the absence of Hartree-induced band inversion (as in fact we will consider in the following section), the shifted bands (↑ bands) are fully filled at ν flat = +2 and the unshifted bands (↓) and the dispersive bands contain the remaining ν -2 charge. With the Hartree-induced gamma point inversion, the two sets of the flat bands (shifted -↑, unshifted -↓) become partially filled near ν flat = +2. The shifted (↑) band is mostly filled and the unshifted (↓) is mostly empty. This simple approach qualitatively reproduces the effect of a cascade at |ν flat | ≈ 2 under the assumption that the specific nature of the cascade state (i.e. spin or valley polarized) is irrelevant for the consideration of the total filling. We caution, however, that a cascade transition is an effect originating from an interplay between Hartree and Fock corrections (see Ref. (32) for further discussion), and that Fock corrections, which we neglected so far, can give rise to many effects.
The most crucial ones, including bandwidth broadening (63,(65)(66)(67) and gap opening in the flatband subsystem (10,11,52,53,60,65,68), may actually affect the charge distribution across the different subsystems. These effects are neglected in the current analysis-an approximation we will justify in the following section.
The cascade, as shown in Fig. S12G, allows charge to enter the unfilled flat bands more easily compared to the uncascaded ground state. This behavior is expected since a cascade minimizes the contribution the Hartree term by redistributing charge away from parts of the flat bands which overlap more strongly with the Hartree potential-in particular, for the parameters considered here and within the relevant range of filling factors, the cascade is the ground state solution. Note, however, that while including only the Hartree correction is sufficient to initiate cascade, Fock must be included in order for it to persist over the experimentally observed range (see Ref. (32) for further discussion on the interplay of Hartree and Fock corrections and the onset of cascade).
To quantitatively estimate the fillings of the different subsystems, it is necessary to parametrize the strength of the Coulomb interaction, e.g., the dielectric constant ϵ that enters into Eq. (S26). Although, in principle, the dielectric constant is fixed primarily by the substrate and any interaction corrections can be accounted for via a self-consistent treatment, in practice (28,32,63,66), it can be treated as a fitting parameter. If a bare value of the interaction is used, then the resulting interaction corrections are too large and lead to unobserved predictions (15,32). We use the cascade near ν flat ≈ +2 to constrain ϵ. We identify ν flat = +2 with the experimental onset of the cascade transition, which occurs near ν ≈ +2.15; see Fig. 4E. By choosing ϵ so that ν flat reaches +2 at the same point the total filling ν reaches +2.15 (see Fig. S12E), we find ϵ ≈ 11.15. We emphasize that although this is an approximate fitting relying on the particular model of a cascade, the Hartree-induced flat-band energy shift is a robust and important effect (Fig. S12F). Using the value of ϵ ≈ 11.15, we find that at ν ≈ +5, the flat bands are filled to approximate ν flat ≈ +3.8 (see Fig. S12G), further demonstrating the implausibility of scenario (i).
We note that due to the hybridization of different bands under a finite displacement field, the assignment of flat TBG-, dispersive TBG-and MLG-like subsystem becomes, to some degree, arbitrary as bands start to hybridize. For the purpose of qualitative discussion, however, we can evaluate the spectral overlap of each finite-field eigenstate with the zero-field basis and assign a label of "flat/dispersive TBG-like/MLG-like" based on the largest overlap. Within this convention, we find that the displacement field enables easier charge accumulation in the "flat" subsystem as opposed to the "non-flat" subsystems, thus suppressing the total filling range over which superconductivity can reach (see Fig. S12H).
Preferential filling of the dispersive TBG-and MLG-like subsystems is also enabled by other interaction effects, for example, gap opening due to the Fock correction. This term, as mentioned previously, also plays a key part in the symmetry-breaking cascade as well as band broadening. While a careful microscopic analysis of Hartree and Fock effects in multilayer devices is necessary, here we introduce a simple phenomenological model intended to capture the qualitative effects of Hartree and Fock corrections on the filling of the non-flat subsystems. We hope that the simple parametrization of this model can be used as a benchmark for its validity against a more rigorous analysis.
To mimic the effects of Hartree and Fock, we add two additional ingredients to the noninteracting model of Eq. (S13): a constant energy shift of the flat-band subsystem ∆ H and an intralayer sublattice potential
A schematic of the effect of the two added terms on the band structure is shown in Fig. S12D. We further mimic a cascade by shifting two copies of the flat-band subsystem below the Dirac points of the band structure and closing the 'Fock' gap in the cascaded bands.
We present the result of this analysis in Fig. S12I, where the filling factor of the flat-band subsystem versus ∆ H and ∆ F is plotted for a fixed total filling ν = +5 (corresponding to the edge of the superconducting dome in TPG). We find that in order for scenario (i) to apply, i.e. ν flat ≈ +3 at ν = +5 in TPG, the corresponding parameters are unrealistic. Especially, a ∆ F ≈ 20 meV would yield an insulating gap of 40 meV, which far exceeds a typical correlated insulating gap of few meV experimentally observed in the context of TBG (59)(60)(61)(62). Notably, for a ∆ F below ∼ 5 meV, there is negligible effect of the Fock gap on the flat-band filling redistribution, consistent with our earlier Hartree-only approximation.
In the above analysis of Hartree-induced band shifting and deformations, we only focused on the effect of the in-plane Hartree potential, assuming a homogeneous charge distribution across the layers. It is known, however, that in multilayer (untwisted) graphene, interlayer screening due to inhomogeneous charge distribution over different layers can be important (69,70). To check the importance of this effect, we follow the self-consistent treatment in Ref. (70) and determine the screening-induced potentials on each layer. Particularly, in the language of the Hamiltonian in Eq. (S5), we show that even at zero external displacement field, the induced potentials are not uniform on all layers and can lead to hybridization between the flat TBG-like and the MLG-like subsystems as presented above in section 4.4.3.
We model the graphene layers as parallel plates with zero thickness and respective electron charge densities en i . The local displacement field between the ith and (i + 1)th layer is then given by F i,i+1 = e( i j=1 n j -N j=i+1 n j )/(2ϵ), where N is the total number of layers, and we take ϵ = 11.15, the same as the in-plane dielectric constant determined above. Here, we use the in-plane dielectric constant value because electron densities are delocalized over multiple layers and cannot be simply treated as some classical charge on a particular layer. Thus, the effective dielectric constant should be much larger than the vacuum's value. The above local displacement field produces the local potential difference between the ith and the (i + 1)th
, where ν j is the filling fraction projected to the jth layer, d is the interlayer distance, and L M is the moiré lattice constant.
The role of this mechanism is shown in Fig. S13. We find that the self-consistently generated potential differences shift the flat bands upwards in energy. Similar to the in-plane Hartree correction, it enables further charge filling of the dispersive TBG-like bands, which is in line with scenario (ii) for the extended TPG superconducting pocket. We stress that unlike the inplane Hartree correction, the out-of-plane Hartree term leads to the hybridization of the different subsystems. This hybridization, in addition to other effects described in the previous section, may facilitate symmetry breaking or a breakdown of an approximate assignment of flat and dispersive TBG-like bands (see the discussion above concerning band mixing, section 4.4.3), in line with the condition for scenario (iii).
Here we present several scenarios that can result in the superconductivity of TPG extending to ν ≈ +5, and discuss these scenarios in the context of experimental observations. We note that in the discussion below, ν denotes the total number of electrons per moiré site, and ν flat denotes the number of electrons per moiré site added to the flat TBG-like bands. For TBG and TTG, the strongest superconducting pockets normally start from |ν| = 2 and end around |ν| = 3. Therefore, a conventional scenario would suggest that TPG could behave in a similar way, i.e., flat TBG-like bands are filled to ν flat = +3 when superconductivity is diminished at ν = +5. This scenario implies that the additional two electrons per moiré site are distributed in the dispersive TBG-and MLG-like bands due to the interaction effects discussed in section 4.5, with a large portion of the charge carriers being hosted by the dispersive TBGlike bands. Since vHs of the dispersive TBG-like bands are normally found around half filling, the corresponding Hall density signatures are expected to occur at the same filling, i.e., ν = +5 in this scenario. However, in the experiment we observe vHs signatures originating from the dispersive TBG-like bands near ν ≈ +6 instead (see Fig. S10). This line of reasoning allows us to completely rule out scenario (i), therefore, we conclude that superconductivity exceeds flat-band filling ν flat = +3 for electron-doped TPG. As a result of interactions, a fraction of electrons are preferentially distributed in the dispersive TBG-and MLG-like bands. It is therefore possible that for total filling of ν ≈ +5, the flat TBG-like bands are filled close to ν flat ≈ +4, with the extra one electron per moiré site being distributed in the other bands. We explored this possibility in more detail in section 4.5, which shows the filling correspondence between ν flat and ν for various interaction terms and dielectric constants (see Fig. S12). In this scenario, the modeling suggests that the filling of the flat bands is nearly four (ν flat > +3.8), which is well outside typical TBG behavior. The last scenario suggests either that the flat TBG-like bands are fully filled before the suppression of superconductivity, in which case superconductivity would exist in the more dispersive bands, or that the distinction between the different TBG-and MLG-like bands breaks down due to hybridization (i.e. mixing), even at D = 0. As discussed in previous sections, such mixing between flat, dispersive TBG-and MLG-like bands can happen when mirror symmetry is broken. Moreover, layer-to-layer charge inhomogeneity (see Fig. S13 and section 4.5.3) or distant-layer coupling (see section 4.4) allow for band hybridization even in the presence of mirror symmetry.
In the context of scenario (iii), we speculate on the polarized or un-polarized nature of the active bands in the regime +4 < ν < +5. One simple possibility is that the dispersive TBGlike bands spontaneously break the flavor symmetries on its own, with the flat bands playing relatively little role. Another possibility is that hybridization obviates the distinction between flat and dispersive bands such that flvor polarization is allowed to persist far beyond. Further experiments will be needed to help ascertain the extent of flavor polarization that persists to ν = +5.
Experimentally, starting from low D fields, we observe a drop in Hall density at ν ≈ +4 which surprisingly does not affect superconductivity in any abrupt way (superconductivity continuously evolves and is present until ν ≈ +5). As the D field is increased, this Hall density drop is gradually replaced by a transition where Hall density changes sign (Fig. S9). The high Dfield transition can be interpreted as a 'gap' feature emerging in the band structure similar to TTG (3). Further measurements of R xx show that the corresponding ν ≈ +4 feature does not shift with temperature (Fig. 4A) and is significantly broadened at high B fields, resembling the feature associated with the flat-band gap in TTG (Fig. S2, A andE). These observations indicate that the ν = +4 feature is naturally explained as either marking the end of the flat bands or resulting from band details due to hybridization, which is in line with the scenario (iii). In this context, the alternative possibility that ν = +4 corresponds to a flavor-polarization reset at ν flat = +3 is highly unlikely. Finally, we note that this line of argument cannot fully rule out scenario (ii) due to the potential presence of small dispersive pockets in the flat bands that may remain unfilled near ν = +4.