Superlattice structures in twisted bilayers of folded graphene

The electronic properties of bilayer graphene strongly depend on relative orientation of the two atomic lattices. Whereas Bernal-stacked graphene is most commonly studied, a rotational mismatch between layers opens up a whole new field of rich physics, especially at small interlayer twist. Here we report on magnetotransport measurements on twisted graphene bilayers, prepared by folding of single layers. These reveal a strong dependence on the twist angle, which can be estimated by means of sample geometry. At small rotation, superlattices with a wavelength in the order of 10 nm arise and are observed by friction atomic force microscopy. Magnetotransport measurements in this small-angle regime show the formation of satellite Landau fans. These are attributed to additional Dirac singularities in the band structure and discussed with respect to the wide range of interlayer coupling models.

T wo-dimensional monolayer graphene exhibits outstanding intrinsic electronic properties 1 , paired with considerable possibilities of manipulation because of high robustness and accessibility at the surface. Furthermore, its atomically thin nature makes it extremely sensitive to ordered as well as disordered potential fluctuations in its vicinity. To this end, several groups have recently observed Hofstadter's butterfly 2 in the spatially modulated potential landscape of graphene on hexagonal boron nitride (hBN) [3][4][5] . Equally fascinating effects can be expected of twist induced moiré superstructures (Fig. 1a) in rotationally mismatched graphene double-layers [6][7][8][9] . Apart from the common Bernal-stacked variety 10 , these twisted bilayers (TBGs) constitute a whole new field of their own: Layers of large rotational mismatch effectively decouple [11][12][13] , exhibiting reduced Fermi velocities for decreasing interlayer twist in many cases 11,[14][15][16][17] . At the smallest angles, totally different electronic structures are expected 7,9,15,18 . In recent years, TBGs of various angles have been grown 19 and optical studies 20,21 as well as scanning tunnelling spectroscopy 11,[22][23][24] were performed on samples of different interlayer twist, revealing, for example, low-energy van Hove singularities (VHSs) and charge density waves. However, there is few systematic work on electronic transport so far, focusing only on large 17,25,26 or disordered small-angle systems 27 . Here we present a study on high-quality folded graphene monolayers of different twist angles y, the smallest of which lead to novel transport features in the form of satellite Landau fans, caused by twist-induced long-wavelength superlattices.

Results
Folded graphene samples. Our folded layers are obtained by mechanical manipulation via atomic force microscope (AFM) or incidental flip-over during the exfoliation of natural graphite. Whereas the method of micromechanical cleavage is known to yield flakes of high crystalline order, a successive folding step induces little further contamination between layers, thus providing TBG of the highest quality. A further advantage of folded samples is that from sample geometry alone, the twist angle y can be estimated (see Fig. 1c): as graphene is most commonly terminated by armchair-or zigzag-edges, alternating in 30°steps 1,28 , an according set of straight edges provide a crystallographic reference direction. Interlayer twist then relates to the angle j between this reference and the folding edge by y ¼ 2 Á j. Note that this principle is not applicable without the common folding edge between layers. Owing to graphene's sixfold symmetry, y may be projected into the range of 0°oyr30°(see Methods section). In case of the shown optical image (Fig. 1b), the twist angle can thus be narrowed down to 1.5°±0.5°. Figure 1d shows AFM topography data of the same sample. Mono-and bilayer regions can clearly be distinguished, as analysed in the cross-sections in Fig. 1e. It is worth noting that we could not observe the formation of bubbles in the folded areas as can be found in samples fabricated by transfer methods 29 . This points towards little contamination between clean interfaces. Note also that at the folded edge of the TBG, a small but distinct elevation of 2.5 Å is present (grey circle, Fig. 1e), indicating a bended but unbroken interconnection between layers, which may contribute to electronic transport, as discussed later. The two rotated lattices can arrange in periodic superstructures reproducing the original honeycomb pattern on a twist-dependent length scale of a being the length of graphene's lattice vector 6,8,15 . The alternating interlayer registry (sublattice AB, BA and AA) modulates coupling and potential landscape in the TBG, which enables experimental visualization of moiré structures in graphene by scanning probe microscopy 22,24,30 . As recently ARTICLE shown 30 , friction AFM serves as a convenient tool for resolving large period superlattices in van der Waals heterostructures. Figure 1f shows a lateral force microscopy scan of the TBG area indicated by the black box in Fig. 1d. The dashed white star marks three distinct symmetry directions in the friction force plot. These are clearly confirmed by the prominent hexagonal pattern in the Fourier transform (Fig. 1h), which points to a trigonal lattice of period l ¼ 9 nm. Figure 1g shows a close up of the yellow box in Fig. 1f with accordingly added unit cells and lattice vector k.
Using equation 1, the resolved lattice matches the moiré structure to an interlayer twist of y ¼ 1.6°, which fits the geometrically estimated value of 1.5°±0.5°very well. Note that the large-scale corrugation in Fig. 1g stems from inherent roughness of the underlying SiO 2 substrate, which is amorphous and unlike hBN does not lead to any periodic superstructures.
Magnetotransport for different rotational mismatch. To investigate the electronic properties of graphene bilayers with different stacking, magnetotransport measurements in perpendicular fields up to B ¼ 13 T were performed. Figure 2 shows examples for AB stacking, small twist angle and larger angle. For each case, the resistance as function of gate voltage and magnetic field (middle row) as well as a cross-section at fixed charge carrier density over B À 1 (top) and the field effect curves at B ¼ 0 T (bottom) are shown. First, we focus on the limits of zero and large relative twist angle: in case of Bernal-AB-stacking (Fig. 2, left), Landau level spectrum and phase of the Shubnikov-de Haas oscillations in B À 1 show a Berry's phase of 2p (ref. 10). For larger rotational mismatch (y43°), the layers effectively decouple and linear monolayer spectra are recovered [12][13][14] , whereas interaction may be seen in a reduction of Fermi velocity, which is strongly angle dependent 11,14,15,17 . For the example in the right column of Fig. 2, the Fermi velocity is found to be 75±5% of the original value v F ¼ 10 6 m s À 1 (ref. 31), as obtained from temperaturedependent Shubnikov-de Haas measurements. These reduced values are used to estimate the rotational mismatch as y ¼ 3.25°±0.75°according to theoretical considerations 15 . In between those cases at approximately 0.3°oyo3°(ref. 9) lies a small range of twist angles, where rich physics like the Hofstadter butterfly are expected to be observable [6][7][8][9] . Indeed, our transport data are most complex for this small-angle regime: in the depicted example (Fig. 2e), a main Landau fan is found to be originating from around À 15 V backgate voltage. The Shubnikov-de Haas oscillations in Fig. 2b indicate charge carrier transport through a graphene monolayer (Berry's phase of p). In addition to this, a second fan arises at a backgate voltage of around 30 V. At zero magnetic field, it is accompanied by a small but notable shoulder in resistance and a dip in the otherwise linear conductance (Fig. 2h). This can be attributed to a large wavelength moiré pattern between the two rotated hexagonal lattices (Fig. 1a,f-h), which induces a periodic potential landscape on a twistdependent length scale. Bragg scattering by the correspondingly small superlattice Brillouin zone will lead to replica satellite Landau fans at higher energies 3-5 , whereas periodically alternating interlayer coupling should result in a more complex electronic spectrum 6,7,9,[32][33][34] .
Small twist angles. In the following, we focus on four TBG samples in the range of 0.3°oyo3°. Figure 3a,b shows the differential resistance dR/dB for two such small-angle devices.
Originating from the charge neutrality point, a monolayer Landau fan with minima at filling factors n is present (n is the charge carrier density, F 0 ¼ h/e the magnetic flux quantum and i an integer). Additional satellite fans can be identified at n ¼ À 1.5 Â 10 12 cm À 2 (Fig. 3a) and n ¼ À 3.5 Â 10 12 cm À 2 (Fig. 3b), respectively. Especially in Fig. 3b, somewhat less defined features arise in between at intermediate charge carrier concentrations. This could either hint to a more complex coupling or stem from limited long-range order, leading to deviations in the moiré period. Following the simple model of an additional superlattice potential, the shift in n is dependent on y and given by where the moiré unit cell (area A) of the fourfold degenerate system gets filled 4 . The expected Landau level picture is mapped out in Fig. 3c. From the transport data of different small angle samples, twist angles are determined by means of equations 1 and 2 and the relation between wavelength and unit cell area of the hexagonal moiré superstructure. The extracted twist angles are summarized in Fig. 3d for samples S1 to S4: The two values for sample S3 have been extracted from transport measurements before and after (S3*) annealing, which changed the overall doping level and notably increased the clarity of the secondary fan but did not affect inter-fan distance.
A more detailed analysis of sample S3* is given in Fig. 4. Figure 4a shows the measured magnetoresistance versus charge carrier concentration and magnetic field. Two Landau fans with origin at charge neutrality and about n ¼ 2.75 Â 10 12 cm À 2 , respectively, can be observed and are fit by a set of slopes similar to the schematic in Fig. 3c. The good match confirms a systematic nature of the observed features. Interestingly, in this sample, only every second slope of the main fan seems to be found in the satellite one as becomes evident at the intersections with the black horizontal lines. Although the main fan is described by fourfold degeneracy as extracted from the expected capacitive coupling constant, the Landau levels of the satellite one are thus eightfold degenerate. This can be verified on Fig. 4b, which depicts our resistance data against backgate voltage and inverse magnetic field: whereas the main fan crosses Landau levels of the second one at every equally spaced horizontal line, the satellite fan only intersects maxima of the first fan at either dashed or solid lines. This coexistence of different degeneracies in the presented example exceeds the simple model of a superlattice potential applicable to graphene on hBN 3-5 and might hint towards a more complex coupling in TBG. Figure 1c shows the differential resistance versus charge carrier density and magnetic field revealing traces of a further Landau fan at higher energies. Note that the origins of the observed Landau level sequences are roughly equally spaced on the axis of carrier density, which points to a systematic relation of the third one to the ones examined above.
Transport over the folded edge. Although coupling in the plane leads to the demonstrated novel electronic properties, also the edges play an important role, especially in case of folded samples featuring a bended interconnection as depicted in the inset of Fig. 5b. Figure 5 shows transport data of a device, in which two monolayer regions were separately contacted on both sides of the folding, thus forcing the charge carriers to either hop between layers or pass the curved folding edge. Such an edge (as observed in the AFM cross-sections in Fig. 1e) induces strong gauge fields 35 and also will exhibit different charge carrier densities. In addition to the usual charge neutrality point peak (filled circle) and features related to a satellite fan at positive gate voltages (white lines), another peak can be observed at around À 25 V. This maximum is independent of the magnetic field B applied perpendicular to the sample plane. As the average bending area will be parallel to B, it is tempting to attribute these feature to the folded edge and the predicted formation of snake states 35 .

Discussion
We have studied folded double-layer graphene devices of various twist angles. At small rotational mismatch, the formation of superlattice structures is observed by AFM, whereas systematic satellite Landau fans emerge in transport measurements. These indicate the generation of secondary Dirac singularities in the band structure at experimentally accessible energies comparable to observations on heterostructures of closely aligned graphene on hBN [3][4][5] . Underlining the important role of interlayer twist, our observations also point to the requirement of long-range order and a well-defined moiré period over a large sample area, as discussed later on. A first intuitive measure for the visibility of superlattice effects in transport measurements is the mean free path l m in relation to the moiré period l, which is proportional to the charge carrier mobility m, readily derived from the field effect measurements at zero magnetic field. The examined samples exhibit typical values of mB2,500 V s cm À 2 , corresponding to a mean free path of l m ¼ 46 nm at a charge carrier concentration of 2.5 Â 10 12 cm À 2 . As l m exceeds the superlattice period, which is around 10 nm in the examined samples, the requirements for an observation in transport should be given. The superlattice effects in graphene on hBN [3][4][5] are significantly more pronounced than in our systems, which can be explained by the high charge carrier mobility in devices on a smooth hBN substrate, exceeding the one in our systems by at least an order of magnitude 4 . Although the mean free path may determine visibility and sharpness of the superlattice features, another factor seems to be responsible for clarity and definition of the observed satellite Landau fans. This ARTICLE becomes evident for sample S3, which has been measured before (S3) and after (S3*) an additional outside annealing step. The satellite features before had a frayed appearance, similar to the examples in Fig. 3a,b. Afterwards, a well-defined Landau fan emerges as depicted in Fig. 4, whereas the mobility remains virtually unchanged. A possible explanation for this would be an improvement in global interlayer registry: as local strain can result in slight variations of interlayer twist over the sample area, causing comparatively large deviations in moiré unit cell (equation 1) and electronic structure, the annealing process and subsequent cooling may have led to a more homogenous interlayer registry, thus improving the clarity of superlattice transport signatures. In a similar vein, as such signatures have not been reported on epitaxial grown bilayers (which often exhibit small twist angles around 2°(refs 12,27)), the high crystal quality of exfoliated and folded samples is likely to play a crucial role in the observed superlattice physics.
Although the observed transport signatures as well as a large wavelength moiré superstructure between two hexagonal lattices constitute striking similarities to the hBN case, the physics in a graphene bilayer are somewhat more complex: in addition to the common superlattice potential, there exists a periodically alternating coupling between the two electrically active graphene layers. In the absence of coupling, the description of a TBG in momentum space is given by two rotated sets of Dirac cones displacing the K-points of two layers by DK ¼ 2 K sin(y/2) (K being the magnitude of the vector K to the Dirac point) 15,32 . Accounting for interlayer coupling, theoretical models predict an angle-dependent reduction of Fermi velocity v F and merging of the displaced Dirac cones in low-energy VHSs [14][15][16] . Landau quantization in such a system would yield a single layer sequence of doubled degeneracy below and a Bernal bilayer sequence above the VHS 36,37 . The energy E VHS at the singularity and the transition in quantization can be estimated as E VHS ¼ v F ' DK/2 À t y , v F being the Fermi velocity and t y the interlayer hopping amplitude 14,22,24 . Although the latter may vary greatly, depending on the system 22 , a value of t y D0.1 eV seems to be a good estimation 14,24 . Using the single-layer relation between charge carrier concentration n and energy E, n ¼ E 2 /(v F 2 ' 2 p), the singularity should occur around n ¼ 2 E VHS 2 /(v F 2 ' 2 p) (the additional factor 2 accounts for the extra twofold degeneracy below the VHS). For the device exhibiting a larger twist angle, as presented in the right column of Fig. 2, the experimentally probed range lies well below this transition, which would be around V bg ¼ 107 V above or below charge neutrality. This is in agreement with the observation of single-layer behaviour in the measurements, whereas stronger influence of the backgate in the lower sheet leads to asymmetric charge distributions and separation of two distinguishable Landau sequences. Together with an observable reduction in Fermi velocity, these measurements are in accordance with the above model. However, looking at the systematic formation of multiple Landau fans in the small-angle regime, another explanation is needed. Chu et al. 32 point out the possibility of the coexistence of interlayer coupling phenomena and superlattice Dirac points generated by periodic potential fluctuations like in graphene on hBN. Mele et al. describe a coupling mechanism that suppresses the renormalization of Fermi velocity as well as the formation of a VHS, producing second-generation Dirac singularities at the crossing points instead 33,34 . The latter model is of special interest for the case of very small angles, as the scenario of a continuous reduction in v F is no longer valid below a certain angle (depending on the model 15,16 ) leading, for example, to flat bands and a localization of charge carriers 15 . Both models 32,33 could explain the interesting observation of multiple Landau fans and the coexistence of different degeneracies in one spectrum as suggested in Fig. 4. Our findings bear further evidence to the rich and complex physics in twisted graphene bilayers especially at the smallest angles, encouraging more theoretical and experimental work on the topic. Moreover, the presented transport signatures suggest that an emergence of electronic superstructure effects in double-layer systems is not limited to ultraclean graphene-boron nitride samples and could also be expected to appear in other upcoming van der Waals heterostructures 38,39 .

Methods
Sample preparation. Our samples were prepared by exfoliation of natural graphite and placed onto Si/SiO 2 substrate. Monolayers and flipped over TBGs were selected via optical microscopy. A multimode AFM with diamond-coated tip of high spring constant (about 40 N m À 1 ) was utilized to fold further samples by programmed tip movement over the sample edges.
Devices were processed using standard electron beam lithography techniques to define and evaporate chromium-gold contacts.
Depending on geometry of the folded sample, contacts were placed in varying positions to the folded edge, which led to different contributions of the bend (the strongest being discussed in Fig. 5).
Outside annealing under argon atmosphere at 400°, before the measurements as well as in situ annealing under vacuum inside the cryostat has been carried out.
Atomic force microscopy. The resolution of the moiré superlattice as shown in Fig. 1f was acquired by friction AFM. A Multimode II with a J-type scanner was used in ambient conditions at a constant, regulated temperature. Several hours of pre-scanning were performed to reduce thermal drift. Best results were obtained with triangular Pyrex-Nitride probes of spring constant k ¼ 0.3 N m À 1 and tip radius rr10 nm. A low set point in the attractive regime was chosen at zero proportional gain and close to zero integral gain. Piezo xy-offset was held at zero. A scan rate of 10 Hz proved to be high enough to reduce thermal noise and low enough to keep piezo xy-oscillations small. The depicted example in Fig. 1f shows the lateral deflection signal in retrace direction at a scanning angle of 50°.
Measurement setup. The transport measurements were performed in a 4 He bath cryostat using DC and low frequency AC setup. A source current was driven through the sample and voltage measured in two-and four-terminal setups. The measurement shown in Fig. 5 was acquired at 2 K, all other shown data were acquired at 1.5 K.
Measurements on the Bernal bilayer presented in the left column of Fig. 2 as well as the small-angle sample S4 presented in Fig. 3b were conducted in two-terminal configuration. All other measurements were performed in a multi-terminal setup. Whereas for sample S1, S2, S4 at least one of the contacts lies on the double-layer system, the device geometry of S3 is in such a way, that two contacts lie on each side of the folded monolayer, allowing four-terminal measurements across the folded edge.
Analysis. The capacitive coupling constant a between sample and backgate is estimated via the parallel-plate capacitor model. Thickness (d ¼ 330 nm) and dielectric constant (e r ¼ 3.9) of the oxidized silicon top layer on the used substrates yield a value of a ¼ 6.53 Â 10 10 V À 1 cm À 2 . Fitting the observed Landau level sequences to i Á g ¼ a DU/[B/(F 0 )] with g as factor of degeneracy, g is readily determined, as applied to Fig. 4a.
Owing to the sixfold symmetry of graphene, a twist angle y'' A[0°, 360°) can always be projected into an interval of y' A [0°, 60°] by y' ¼ y'' mod 60°. Although, depending on the axis of rotation, a 60°twist can mean the difference between AB-and AA-stacking, there is small difference in the overall moiré pattern between y' and 60°-y' for 0°aya60°in the low-angle regimen 15 , allowing for a further projection to y A (0°, 30°] by y ¼ 60°-y' for y'430°.