Moiréless correlations in ABCA graphene.

Atomically thin van der Waals materials stacked with an interlayer twist have proven to be an excellent platform toward achieving gate-tunable correlated phenomena linked to the formation of flat electronic bands. In this work we demonstrate the formation of emergent correlated phases in multilayer rhombohedral graphene--a simple material that also exhibits a flat electronic band edge but without the need of having a moiré superlattice induced by twisted van der Waals layers. We show that two layers of bilayer graphene that are twisted by an arbitrary tiny angle host large (micrometer-scale) regions of uniform rhombohedral four-layer (ABCA) graphene that can be independently studied. Scanning tunneling spectroscopy reveals that ABCA graphene hosts an unprecedentedly sharp van Hove singularity of 3-5-meV half-width. We demonstrate that when this van Hove singularity straddles the Fermi level, a correlated many-body gap emerges with peak-to-peak value of 9.5 meV at charge neutrality. Mean-field theoretical calculations for model with short-ranged interactions indicate that two primary candidates for the appearance of this broken symmetry state are a charge-transfer excitonic insulator and a ferrimagnet. Finally, we show that ABCA graphene hosts surface topological helical edge states at natural interfaces with ABAB graphene which can be turned on and off with gate voltage, implying that small-angle twisted double-bilayer graphene is an ideal programmable topological quantum material.

Two-dimensional (2D) van der Waals heterostructures with an interlayer twist have provided a new avenue for observing emergent tunable many-body electron phenomena. Recent experimental realizations include twisted bilayer graphene (tBG) near the so-called “magic angle” of 1.1° (1–3), twisted double-bilayer graphene (tDBG) (4–6), ABC trilayer graphene on near-perfectly aligned hexagonal boron nitride (hBN) (ABC-tLG/hBN) (7, 8) and transition-metal dichalcogenide heterostructures (9–12) [with predictions on a variety of other systems (13, 14)]. All of these systems host an interplay of two phenomena––the presence of one or more van Hove singularities (which we colloquially refer to as “flat bands” henceforth) at low energy where the density of states is sharply peaked, and the existence of a moiré pattern that creates a unit cell that is about a hundred times larger than the carbon–carbon nearest-neighbor distance in graphene. The large number of electrons with quenched kinetic energy make the flat bands conducive to interaction-driven phases (15). The enlarged moiré unit cell is thought to reduce both the flat-band bandwidth and the interaction energy scales, and also introduces easily accessible integer fillings that create Mott-like insulating states (1–12), the relation of which to nearby superconductivity is debated. A natural question that arises from all of these works is whether the moiré pattern is a necessary condition for the observation of correlated many-body phases, or whether it is simply sufficient to further reduce the flat-band bandwidth and hence the kinetic energy in the heterostructure.

In this regard, multilayer rhombohedral (ABC) graphene offers a different perspective toward achieving a flat-band edge without the use of a moiré potential (16). Indeed, in a seminal work (17), it was theoretically shown that the low-energy band structure of multilayer rhombohedral graphene has a sharply peaked density of state (DOS), with the band structure (where N is the number of layers) at low energy in the nearest-neighbor hopping approximation. This implies a peak in the DOS at charge neutrality in this material for , with an appreciable fraction of the entire band within this peak (18). Indeed, this physics is already at play in ABC-tLG/hBN (7, 8), where some of the flatness of the bands comes from the intrinsic band structure of ABC graphene, which is then further flattened and isolated by the moiré pattern from the hBN alignment. A facile alternative to flatten the bandwidth without introducing a moiré potential is to simply increase the number of layers of the rhombohedral stacked graphene. Unfortunately, isolating rhombohedral stacked graphene of any thickness is extremely difficult as it is less energetically favorable than the multilayer counterpart, Bernal stacked graphene. Since the difference between rhombohedral and Bernal graphene is simply a lattice shift, and the interlayer van der Waals forces are weak, it is well known that rhombohedral graphene reverts to the Bernal form when samples are processed with heat, pressure, or while performing lithography (19). In this work, we show that twisting two sheets of tDBG by a tiny (<0.1°) angle is a simple and robust method to create large area (up to micrometer-scale) rhombohedral graphene of four-layer thickness (ABCA graphene). We present gate-tunable scanning tunneling microscopy and spectroscopy (STM/STS) measurements at 5.7 K on these regions. We show that correlated phases can be achieved without the need for a moiré pattern and that rhombohedral graphene has unique topological properties.

Results

In the limit of small twist angles, it is well known that two twisted layers of graphene will prefer to form extended regions where the lattice is commensurate and Bernal stacked (called AB or BA stacking) (20), separated by sharp domain walls where the stacking order changes from AB to BA (21, 22). Where the domain walls intersect, the stacking order is locally AA, which is also where the low-energy electrons are localized (15). We begin by considering the atomic stacking nature of small-angle tDBG, as shown in Fig. 1 . The easiest way to think of tDBG is to consider the middle two layers which form the same structure as tBG. The rest of the structure of tDBG can then be visualized by the fact that the two upper and two lower layers are both Bernal stacked. Using this visualization technique, it becomes clear that the equivalent of the AA stacking location in tBG is a ABBC stacking in tDBG. One can then displace the top bilayer by one atom in either direction to visualize the other commensurate stacking sites in a manner exactly analogous to tBG. The consequence of this as shown in Fig. 1 is the formation of ABAB (Bernal) and ABCA (rhombohedral) stacked four-layer graphene in the commensurate regions of the sample. Because it is energetically favorable for commensurate regions to maximize in area, small-angle tDBG has large, uniform extended regions of rhombohedral and Bernal graphene that push any effect of the tDBG moiré superlattice to the interfaces of the commensurate domains. At small angles, tDBG thus forms large, uniform rhombohedral domains that are decoupled from effects of the large-wavelength small-angle tDBG moiré. The final stacking arrangement that arises is then shown in Fig. 1. The rhombohedral graphene is protected by the global symmetry enforced by the twist angle, making it stable and easy to produce. While in this work we show this case for small-angle tDBG, the atomic stacking arguments at hand will also extend to any angle including the magic angle of tDBG where the ABCA sites will likely play a role in electron localization along with the ABBC sites. These arguments also apply to monolayer–bilayer twisted trilayer graphene where the commensurate stackings become ABA and ABC [° tDBG is another interesting case where commensurate sites will be ABCB and BCBA].

Fig. 1.

Rhombohedral domains in small-angle tDBG. (A) Side view and top view of atomic-stacking considerations in a tDBG moiré pattern. Due to the enforced Bernal stacking of the top and bottom layer, the AA-site stacking equivalent in tDBG is ABBC. One can translate the bilayer–bilayer atomic alignment by one atom in either direction to visualize commensurate stackings around the ABBC site (as shown by arrows in B for a moiré unit-cell diagram). While the monolayer–monolayer case (which can be seen by the middle two layers) becomes AB/BA, when the top and bottom tDBG layers are considered the commensurate stacking sites become ABAB (Bernal) and ABCA (rhombohedral). (B) Cartoon moiré of small-angle tDBG with two example directions of atomic-stacking shift considered in A. (C) A nanoinfrared image of a large area of a small-angle tDBG sample. (D) An STM topographic image of small-angle tDBG (setpoints of 300 mV and 100 pA). (E) STS LDOS of ABCA and ABAB graphene domains in small-angle tDBG taken at the center of one of each domain (setpoints of 300 mV and 150 pA with a lock-in oscillation of 2.5 mV). (F) STS LDOS map at the energy of the ABCA flat-band edge (setpoints of 400 mV and 150 pA with a 2.5-mV oscillation).

We first use nanoinfrared imaging to visualize small-angle tDBG on a large scale. Fig. 1 shows a nanoinfrared image of a large-area region of small-angle tDBG. We see that there are two distinct sets of domains (bright and dark in the image) with quasihexagonal order. The average size of the domains in this sample is of order 200 nm on a side, though several larger domains of over a micrometer in size are seen, presumably due to local strain or pinning of particular stacking sites. We have observed such domain structures of hundreds of nanometers in width in every small-angle tDBG sample. The contrast between the two domains clearly shows that the two domains have distinct electronic responses, as we might expect for Bernal and rhombohedral graphene. Having seen the optical signatures of small-angle tDBG structure on a large length scale, we turn to STM/STS to study the atomic-scale spectroscopic properties. In Fig. 1 we show an STM topographic image of a homogeneous small-angle tDBG region with triangular domains that are about 220 nm per side. This can directly be contrasted with a small-angle tBG STM topography in . In tBG, the STM topographic contrast is dominated by the AB/BA domain walls and the AA site centers, while the two commensurate AB and BA domains have identical contrast. In the case of tDBG, the contrast differs between the two sets of triangular domains, ABCA and ABAB.

We next turn to the electronic structure of the domains. In Fig. 1 we show STS local DOS (LDOS) on the ABCA and ABAB domains elucidating their starkly different electronic properties. While the ABAB domains show a gradually increasing density of states away from the Fermi level (consistent with measurements on multilayer graphite), the ABCA domains have an extremely sharp enhancement in the DOS near the Fermi level. In Fig. 1 we show an STS LDOS image at the energy of the LDOS peak. This image shows that the contrast between the ABCA and ABAB domains is uniformly observed in the low-energy electronic structure, and thus the peak in the LDOS in Fig. 1 is simply a property of the ABCA graphene itself and not due to confinement or impurity effects. Furthermore, the uniform spatial electronic structure proves that these ABCA domains have spatially maximized their commensurate area, pushing the unnatural atomic stackings of the spatially evolving tDBG moiré to the domain walls and decoupling the ABCA graphene from the effects of the superlattice potential. Thus, although we use a tDBG moiré to achieve the ABCA graphene regions, they are no longer affected by the much larger moiré, allowing the study of ABCA graphene without any superlattice potential affecting its low-energy electronic properties. From the topography and STS map in Fig. 1 , it can also be seen that the ABAB Bernal domains are convex whereas the ABCA domains are concave. This domain shape can simply be understood as a competition between the desire on the part of the ABAB domains to maximize their area (being more energetically favorable) and the desire to minimize the length of the domain wall (similar to surface tension). This behavior is also different from tBG, where the two AB/BA domains are equivalent, leading to simple straight domain walls as seen in and previous works (22, 23).

Having clearly identified the large regions of uniform ABCA graphene in our sample, we now focus on the low-energy gate-dependent spectroscopic structure. It is important to note that in STM measurements, the spectroscopic signals are dominated by the first conducting layer of the sample, and hence by the top layer of ABCA graphene. To begin, we compare the surface-measured experimental LDOS (Fig. 1) to layer-dependent theoretical LDOS as calculated by density-functional theory (DFT) and shown in Fig. 2. The low-energy density of states is concentrated on the outer layers (i.e., layers 1 and 4) and we find that when examined over a large energy range, the LDOS predicted by DFT for the top layer matches the observed spectroscopy in Fig. 1 strikingly well, including the high-energy bands appearing as steps in the spectrum.

Fig. 2.

Displacement-field-dependent LDOS in ABCA graphene. (A) Layer-dependent LDOS as calculated by self-consistent DFT. (B) STS LDOS curve at a gate voltage of 50 V and top-layer DFT-calculated DOS at 0.8 V/nm. (C) Cartoon of the STM/STS measurement geometry under an applied back-gate voltage. The sample is biased relative to the grounded tip. The back gate dopes the sample and inherently also induces a displacement field. (D) Displacement-field-dependent low-energy band structure of ABCA graphene as calculated by DFT. (E) A series of experimental STS LDOS (solid) and DFT top-layer LDOS (dotted) curves at various displacement-field magnitudes. (F) Extracted displacement-field gap sizes for 251 STS LDOS curves equally spaced between our displacement-field extremes. DFT-calculated displacement-field gap sizes are overlaid. STS measurements in B and E were at setpoints of 400 mV and 150 pA with a lock-in oscillation of 2.5 mV.

We now turn to the effect of a gate voltage on our measurements. In Fig. 2 we show a curve at our maximum positive gate voltage of 50 V (limited by gate leakage). We see that the shape of the measured spectrum has undergone significant changes. The low-energy peak and high-energy band edges are still present in the spectrum (although they have shifted due to doping), but in addition a clear gap has emerged below the Fermi level. To better understand the origin of this gap, in Fig. 2 we show a model for our geometry and basic electrostatics. In our measurement, the sample is biased relative to the tip, while the gate voltage is applied relative to the sample. The application of a back-gate voltage results in both electrostatic doping of the sample as a whole, as well as an electric field applied across the layers of the sample. In trilayer ABC layer graphene (24), it has indeed been shown that electric fields can cause a gap to form. To better understand this effect in four-layer ABCA graphene, we calculate the electric-field dependence of the band structure by means of self-consistent DFT in Fig. 2. Upon application of a displacement field to the sample, a gap emerges in the ABCA graphene band structure as seen in Fig. 2, which shows the DFT-calculated top-layer DOS on top. While the band structure of the top layer is electron–hole asymmetric in an applied displacement field, the bottom layer has a band structure that is approximately electron–hole mirror-symmetric with respect to the top layer. We also note that in the presence of a displacement field, the low-energy bands develop a topological character, with each of the valleys in the conduction and valence bands having a Chern number of 2, but with opposing signs (25, 26) ().

For the STS curve at +50 VG in Fig. 2, we estimate our displacement field to be 0.85 V/nm (we assume the tip to be one electrode of a parallel plate capacitor held at 0 V and account for intrinsic sample doping which displaces charge neutrality and the zero-field condition). The DFT prediction for the top-layer DOS with an electric field of 0.8 V/nm (solid black curve in Fig. 2) matches the experimental spectrum well. We note that in the DFT calculation, the electric-field value is fixed at infinity, and the charge on the layers is allowed to redistribute itself self-consistently. We have measured the experimental LDOS for a range of gate voltages, a few of which are displayed in Fig. 2. Overlaid on these curves are the DFT predictions at the corresponding electric-field values. As the field is increased, the displacement field gap grows in size. For electric fields pointing from the gate to the tip, the gap is located to the left (lower energy) of the main peak in the LDOS, while for fields pointing in the other direction the gap is located to the right (higher energy). In Fig. 2 we track the displacement-field-dependent gap size for 251 STS LDOS curves extending between our gate-voltage extremes of −75 and 50 V with DFT calculated gap size overlaid. For displacement-field magnitudes larger than about 0.3 V/nm, DFT and experiment are in excellent agreement. The curves shown in Fig. 2 are all in this category. While the match between experiment and DFT is excellent at high fields in Fig. 2, for low fields the single-particle DFT fails to capture the experimental results. In particular, DFT predicts that for electric fields between −0.1 and 0.1 V/nm, the gap produced by the displacement field is negligible. This is a consequence of screening that reduces the field inside the ABCA graphene. The experimental data are in stark contrast to this prediction, always displaying a nonzero gap in this region that is of order 5–10 meV. The fact that the gap matches quantitatively at high positive and negative field negates any possibility of extraneous fields (for example from the tip) affecting this conclusion. We have also checked the effect of the hBN substrate with DFT calculations. Our results show that even though a band gap can be opened by hBN substrate in certain perfectly alignment at zero field, such gap will be closed when we apply certain displacement field (see more discussion in ). This is also inconsistent with our experimental observation. We are thus forced to consider other explanations for the observed gap at charge neutrality in experiment.

To highlight the gap at zero field, in Fig. 3, we show two LDOS curves at small displacement-field values. The curve shown in blue is exactly at charge neutrality (and zero field) with a gap of 9.5 meV centered precisely at the Fermi level. The curve shown in magenta is the smallest gap observed at any value of gate voltage, 4.5 meV in magnitude, measured at a field of about −0.1 V/nm. This gap in the LDOS is consistent throughout the entire ABCA domain as can be seen in Fig. 3 where we plot the LDOS at a sequence of points across a whole ABCA and ABAB domain at charge neutrality. Given the clear disagreement with single-particle DFT calculations, this gap is a consequence of electron correlations. We note that we can rule out strain effects since strain values are minimum, as can be seen in Figs. 1 and 4. Recently, in transport measurements, a large gap of 80 meV has been observed in four-layer ABCA graphene (27). The gap observed in our experiments is an order of magnitude smaller, and we see no evidence for the large transport gap in our microscopic measurements. Our measurements have been repeated on several samples showing ABCA domains, indicating the robustness of the gap values we measure. One question that arises is why ABCA graphene displays such a correlated gap, while ABC graphene has not displayed this phenomenon in previous spectroscopic measurements (28). To answer this question, we quantify the sharpness of the ABCA LDOS peak by extracting the half-width of the most prominent peak in the spectrum at each value of gate voltage. A plot of these widths as a function of gate voltage is shown in Fig. 3. For a range of gate voltages near charge neutrality (gray region in Fig. 3), the presence of the correlated gap in the peak makes it difficult to define a peak width cleanly. At larger gate voltages, we find that the peak width increases linearly with gate voltage, likely a quasiparticle lifetime effect. An extrapolation of the linear behavior to charge neutrality shows that without splitting, the peak would have a half-width of 3 meV when the sample is charge neutral. Previous measurements of ABC graphene have shown a peak width of 10–20 meV near charge neutrality (17) and about 50 meV on a conducting substrate (29). This sharpness of the peak in the DOS of ABCA graphene is far lower than trilayer ABC graphene (as expected in the low-energy approximation) or even the flat bands of magic-angle tBG (24, 30–32) which are about 10 meV in half-width (23). It is important to reiterate that at these small angles in tDBG, the commensurate ABCA regions are energetically maximized and have pushed any effects of the tDBG moiré to the domain walls, thus this finding is purely associated with ABCA graphene. ABCA graphene is thus a simple flat-band system, where no isolated flat bands emerge as the flatness comes from the band edge, where kinetic energy is sufficiently quenched such that a correlation-induced gap emerges without a relevant moiré periodicity.

Fig. 3.

Electronic correlations in ABCA graphene. (A) Two STS LDOS curves at small displacement fields (blue at neutrality and ∼0 V/nm and magenta at −0.1 V/nm) where in a single-particle picture there should be no gap. The van Hove singularity splits and shows a correlated gap of 9.5 and 4.5 meV, respectively. (B) STS LDOS linecut at charge neutrality across an ABCA and an ABAB domain. (C) ABCA van Hove singularity half-width as a function of gate voltage. The shaded region in gray is excluded as the van Hove singularity half-width cannot be cleanly extracted due to splitting. (Inset) An STS spectrum with the peak still far from the Fermi level, exhibiting a 4-meV half-width. STS measurements in A–C were taken at setpoints of 300 mV and 200 pA with a lock-in oscillation of 0.5 mV. (D) Tight-binding mean-field DOS for ferrimagnetic and charge-transfer excitonic insulator ordered states. (E and F) Visual representation of the two ordered states in ABCA graphene. The mostly irrelevant middle layers are shaded. Ferrimagnetic spin order (E) and charge-transfer insulating order (F) are shown by arrows pointing up or down as well as big or small spheres, respectively (shown on select sites only). (G) Asymmetry of peaks in mean-field tight-binding model as a function of U, intralayer, and U′, top–bottom layer interactions. The most likely ordered state in each region is labeled.

Fig. 4.

Gate-tunable topological edge states ABCA graphene. (A and B) STS LDOS map cuts at −20 meV (within the overlapping gap of ABAB and ABCA graphene) and −10 meV (at the flat-band edge peak). (C) Band structure of ABAB graphene at a displacement field of 0.8 V/nm as calculated by DFT with self-consistent charge screening. The size of the blue circles is proportional to the projection of the wavefunction on the top layer. (D and E) STS LDOS map cuts (from same map as A and B) at 37.5 and 22.5 meV. (F) Two-dimensional Fourier transform of D and E (Right) and a linecut through the primary electron scattering peaks (Left). (G) Tight-binding LDOS profiles at fixed spatial intervals crossing an ABAB–ABCA domain wall with their corresponding distances (with respect to the domain boundary) color coded and written for the beginning and end point as well as the dominant in-gap topological edge state. (H) STS LDOS profiles at fixed spatial intervals crossing the ABAB–ABCA domain wall at the represented and color-encoded distances (with respect to the domain wall as determined at higher energies) showing the dominant in gap topological edge state about 10 nm into the ABCA domain, consistent with tight binding in G. STS measurements in A, B, D, E, and H were taken at setpoints of 400 mV and 150 pA with lock-in oscillations of 2.5 mV (A, B, D, and E) and 1 mV (H).

Discussion

Our DFT result at zero displacement field is in consistent with previous calculations (33). In the previous work, Pamuk and coworkers showed that a magnetic gap can be opened in the ABCA multilayer graphene with the DFT calculations if one increases the electron–electron correlations by manually including certain portion of exact exchange in the DFT functionals. In order to understand the nature of the emergent insulating state, we theoretically investigate the ordered states that arise from adding correlation effects in this regime more systematically. To this end, we use a tight-binding model with a Hamiltonian of with the noninteracting part, and treat electron–electron interactions described by within a mean-field approximation (). Note that in doing so we treat both valley degrees of freedom. We first concentrate on short-range Hubbard-type interactions between the up and down spin electrons of the material within each layer by including where i runs over all lattice sites and is the density at site i. We treat the term in a self-consistent mean-field decoupling yielding the results shown by the magenta curve in Fig. 3 (). The correlated (mean-field) state in this case is a ferrimagnetic state within the topmost layer which spontaneously breaks the spatial unitary group [SU(2)] invariance of the system and an antiferromagnetic ordering across the topmost and bottommost layer (opposite ferrimagnetic state in the bottommost layer) as shown in Fig. 3. To have a more complete model, we now add longer-range interactions across topmost to bottommost layers on A and B sites, respectively. This is a valid choice of the next-order interactions as we know (from both tight binding and DFT) that the top and bottom layer dominate the low-energy physics. Toward this end and to keep the model as simple as possible, we add an additional contribution, to H and study the ordered states in the self-consistent mean-field approximation as we sweep U and U′, the intralayer interactions and the top–bottom layer interactions. When turning on U′, at values U′ ∼ U, the ferrimagnetic state discussed above (at U′ = 0) gives way to a charge-transfer (excitonic) insulating state which spontaneously breaks the inversion symmetry between bottom and top layer, spontaneously transferring charge across the A sites of the topmost and B sites of the bottommost layer as shown in Fig. 3 (). Due to the spontaneous charge transfer, the local density of states at the A sites of the topmost layer becomes electron–hole asymmetric in contrast to the LDOS in a ferrimagnetically ordered state (also shown in Fig. 3) that maintains electron–hole symmetry. We calculate a full phase diagram of the system as a function of U and U′, shown in Fig. 3. The color scale on this plot displays the electron–hole asymmetry, which is seen to be generically large in the charge-transfer excitonic insulator state. Based upon this criterion alone, our experimental spectra at most gate voltages match the expectation for the charge-transfer excitonic insulator state. The only exception is the region around where the gap is a minimum, where it is observed to be more electron–hole symmetric. While the electron–hole asymmetry gives us some insight into the nature of the correlated gap, future measurements at high magnetic field can definitively answer whether the state involves magnetic order. We note explicitly that our simplified model leaves room for future investigations, trying to include the effects of correlations on a more quantitative level. Although our model shows results consistent with the experimental observations, other gap-opening mechanisms should be the subject of future studies. Two of such potential candidates (strain and substrate effects) have been ruled out above, leaving correlations as a plausible choice (note that the effects of the substrate might strengthen the charge-transfer mechanism over the ferromagnetic one.

We next look at the consequence of the topological character of ABCA graphene in an applied displacement field. In Fig. 4, we show an STS LDOS image of several ABCA and ABAB domains at an electric field of 0.475 V/nm and −20 meV in energy. We clearly see the presence of edge states near ABCA–ABAB domain walls dominate the LDOS. Fig. 4 (at the flat-band edge peak, −10 meV) and D-E (above the flat-band edge) show more energy slices from the same LDOS map as Fig. 4. In these images, there is no signature of the edge states, which only appear at particular energies and gate voltages. Toward understanding these edge states, we consider the nature of the ABAB domains in comparison to the ABCA domains. In Fig. 4 we show the DFT calculated flat-band structure of ABAB graphene under a displacement field (with the projection on the top layer shown by the size of the blue circles). The band structure is highly layer polarized. While the low-energy band labeled 1 dominates the top layer, band 2 has nearly zero weight. Under the application of a displacement field, a gap is observed on the top layer (between bands 1 and 3), although the ABAB graphene four-layer structure as a whole remains metallic. This is also confirmed by STS measurements in the ABAB regions under a displacement field, Fig. 4. We further confirm this band structure with quasiparticle interference mapping. Fig. 4 show STS LDOS map cuts at two energies in the conduction band of ABAB graphene. The wavelike structure in the ABAB domains originates from backscattering from the ABCA domain walls. Shown in Fig. 4 are 2D Fourier transforms (right) of the LDOS maps in Fig. 4 revealing the dominant hexagonal peaks from the scattering pattern. A line cut through the peaks (Fig. 4 , Left) gives a direct visualization of the top-layer quasiparabolic band as well as a faint signature of the valence band (it dominates on the bottom layer), matching the theoretical band structure under field. Having confirmed the theoretically predicted band structure of ABAB graphene under a displacement field, we now consider its topological properties. As the domain boundary between ABCA and ABAB graphene varies smoothly across a scale of 5–10 nm, strong intervalley scattering is expected to be negligible. We can consider the valley Chern number as a good topology quantity across the domain-wall boundary. An analysis of band 1 yields a valley Chern number of 1, with opposing signs in the two valleys () in ABAB graphene. Thus, when considering the interface between ABCA and ABAB graphene seen in STM, both develop gaps in the presence of a displacement field with a difference in Chern number of 1 for each valley. This in turn implies that there are helical edge states of topological origin at the interface with one mode per valley [in contrast to the interface between AB and BA bilayer graphene where two modes per valley are expected (22, 34, 35)]. Although ABAB graphene remains metallic under displacement field, the system is strongly layer polarized. The topological edge states reside only in the top layer, while the other bulk states localized in other layers. Due to the layer polarization, scattering between the topological edge states in the top layer and other bulk states in ABAB graphene are significantly reduced although such scattering is not prevented by topology. The dominating edges in the LDOS map of Fig. 4 are a direct manifestation of these topological surface helical edge states.

For further insight into the observed interface properties of this topological edge state, we study a finite-size spatial region that includes a sharp boundary between ABCA and ABAB graphene at similar conditions as in experiment using tight binding (). In Fig. 4 we show the results of these calculations for the LDOS at a number of positions across the ABCA–ABAB interface. The calculations show that the topological edge mode is localized around 8 nm away from the interface inside the ABCA domain, characterized by an in-gap state within the overlapping displacement field gaps of ABCA and ABAB graphene. To compare to theory, we show the corresponding experimentally measured LDOS spectra across the interface in Fig. 4. We see remarkable agreement between theory and experiment for the spatial LDOS dependence. Within the overlapping displacement-field gaps of the LDOS, an in-gap state is present which in experiment is located 10 nm away from the domain boundary (36) inside the ABCA region, in close accord with theory. Edge states are composed of linear combinations of the bulk states within each domain. Since the DOS of the ABCA domain is much higher than that of the ABAB domain, the contribution of the bulk states from the ABCA domain is much larger, so the resulting edge states localize more on the ABCA side. These measurements establish that the interface of ABCA–ABAB domains (which can be easily realized by small-angle tDBG) feature surface helical topological edge states which can be turned on and off with displacement field and gate voltage.

Our experiments show that small-angle tDBG moiré patterns are an easy platform to create large, stable, four-layer rhombohedral graphene domains. We commonly observe rhombohedral patches that are hundreds of nanometers to a micrometer in size, making them suitable for any technique including transport where electrodes can be directly placed on the rhombohedral regions. One can additionally explore stacking more bilayers at small angles toward even thicker rhombohedral graphene. Spectroscopy on ABCA graphene reveals a van Hove singularity of 3–5-meV half-width and a correlated gap at neutrality, making it a unique correlated van der Waals system where there is no need for a moiré superlattice to modify the low-energy electronic properties. Further exploration of ABCA graphene can provide answers to certain fundamental questions, for instance regarding the importance of moiré potentials and integer fillings in generating emergent phases such as superconductivity in neighboring regions of the phase diagram. The chiral nature of the tunable topological edge states in rhombohedral graphene also makes it of great interest to observe gyrotropic effects and toward Floquet engineering.

Materials and Methods

Experiment.

Our fabrication of tDBG samples follows the established “tearing” method, using polypropylene carbonate (PPC) as a polymer to sequentially pick up hBN, half of a piece of graphene followed by the second half with a twist angle. This structure is flipped over and placed on a Si/SiO2 chip. PPC polymer was removed by vacuum annealing at 250 °C. Afterward, direct contact is made to the tDBG via µ-soldering with Field’s metal (37), keeping temperatures below 80 °C during the entire process to minimize the chance of layers fully rotating back to Bernal stacking which happens on annealing the structures.

STM measurements are taken in a home-built 5.7 K system at ultrahigh vacuum. STM tips are prepared and calibrated for atomic sharpness and electronic integrity on freshly prepared Au (111) crystals. Samples were measured with multiple tips to ensure consistency of results.

The size of the experimental displacement field gap in Fig. 2 is defined as the distance from the van Hove singularity to the gap edge. We fitted a Gaussian function to the van Hove singularity in order to obtain its energy position. The gap edge is determined by the energy position at which the LDOS intensity reaches a plateau.

Nanoinfrared imaging in this work was performed using a commercial (Neaspec) scattering-type scanning near-field optical microscope. In this measurement, midinfrared light (continuous-wave CO2 gas laser at a wavelength of 10.6 µm) was focused on the apex of a metallic tip. The scattered light was collected by a cryogenic HgCdTe detector (Kolmar Technologies). The AFM tip was excited at a frequency of about 75 kHz, with a tapping amplitude of about 60 nm. The far-field contribution to the signal can be eliminated from the signal by locking to a high harmonic (here we used the third harmonic of the tapping). The phase of the backscattered signal was extracted using an interferometric detection method, the pseudoheterodyne scheme, by interfering at the detector the scattered light with a modulated reference arm.

Tight Binding and DFT.

The band structures and layer-projected LDOS of multilayer graphene with ABCA and ABAB stacking are calculated with first-principle calculations based on DFT as implemented in the Vienna Ab initio Simulation Package (38). The pseudopotentials are constructed with the projector augmented-wave method (39). The exchange-correlation functionals are treated at the local-density approximations (40) level. Plane-wave basis sets are employed with an energy cutoff of 600 eV. A dense k-point grid of 150 × 150 × 1 is adopted in the self-consistent calculations to ensure the convergence of the band structures and the LDOS. A vacuum region larger than 15 Å is added along the z direction to eliminate artificial interactions between slab images. Dipole corrections are included in all the calculations with displacement fields.

To simulate the STS measurement across the ABCA–ABAB domain walls, we calculate the LDOS at different locations across the domain walls with tight-binding model calculations. We model the domain walls with a flake of four-layer graphene with area larger than 200 nm × 200 nm in the real space. In the region with x < −1 nm, the stacking sequence is ABCA. In the region between x = −1 nm and x = 1 nm, the x coordinates of all lattice sites in the top two layers are linearly expanded such that the stackings become ABAB in the region with x > 1 nm. We model the system with the following tight-binding Hamiltonian (41): , where is the hopping parameter between pz orbitals at the two lattice sites and and it has the following form: where a = 1.412 Å is the in plane C–C bond length, c = 3.36 Å is the interlayer separation, n is the direction cosine of along the out of plane axis (z axis), is the intralayer (interlayer) hopping parameter, and is the intralayer (interlayer) decay constant. For the interlayer interaction, only the hopping between atoms in adjacent layers are considered. This tight-binding model has been shown to reproduce the low-energy structure of graphene calculated by local DFT calculations with the following value for the parameters: , , and (41). For the calculation of LDOS, we employ the Lanczos recursive method (42).