Spectroscopic characterization of charge carrier anisotropic motion in twisted few-layer graphene

Graphene, a layer of carbon atoms in a honeycomb lattice, captures enormous interest as probably the most promising component of future electronics thanks to its mechanical robustness, flexibility, and unique charge carrier quasiparticles propagating like massless high energy Dirac fermions. If several graphene layers form a stack, the interaction between them is, on the one hand, weak, allowing realization of various registries between the layers and, on the other hand, strong enough for a wide range tuning of the electronic properties. Here we grow few layer graphene with various number of layers and twist configurations and address the electronic properties of individual atomic layers in single microscopic domains using angle-resolved photoelectron spectromicroscopy. The dependence of the interlayer coupling on the twist angle is analyzed and, in the domains with tri-layers and more, if different rotations are present, the electrons in weaker coupled adjacent layers are shown to have different properties manifested by coexisting van Hove singularities, moiré superlattices with corresponding superlattice Dirac points, and charge carrier group velocity renormalizations. Moreover, pronounced anisotropy in the charge carrier motion, opening a possibility to transform strongly coupled graphene bilayers into quasi one-dimensional conductors, is observed.

Moreover, as a distinctive feature of the chirality, the propagation of Dirac fermions under external periodic potential is expected to be anisotropic, manifesting anisotropic group velocities 19 , anisotropic gaps in electronic structure at minizone boundary (MB) 19 and the formation of superlattice Dirac points 20 .
Experimental view on the anisotropy in the electronic structure of twisted nLG is particularly important for small θ, for which the above mentioned intriguing phenomena are predicted by tight binding models and at the same time first principle calculations become impractical as the supercell contains big number of atoms. Moreover, the results of ab initio calculations for rather large twist angles, suggesting effective decoupling of layers 21 , are in agreement with angle resolved photoelectron spectroscopy (ARPES) from multi-domain samples of nLG grown on SiC 22 . However, scanning tunnelling spectroscopy (STS) 11,12,14 , and recent ARPES 23 studies as well as our results obtained from single domains confirm the presence of interlayer coupling. In the experiment by Ohta and coworkers 23 2LG samples had rather large θ (> 5 o ) and the domains of ~100 μ m size comparable to typical photon beam dimensions in ARPES. The measurements indicate anisotropic gap at MB but non-interacting Dirac cones. Evidence of the velocity renormalization and, hence the interaction of Dirac cones was obtained by Landau level STS in magnetic field 14 , however k-resolved density of states and, consequently anisotropy, cannot be addressed by this method.
Here, we use microscopic ARPES (μ -ARPES) of single few micron sized rotationally misoriented domains of twisted nLG grown on SiC substrate and illustrate in detail the effects of the interlayer coupling on the electronic structure of individual graphene layers. We show that, when several concomitant couplings break the equivalence between the layers (n ≥ 3), the charge carriers may exhibit different transport characteristics in different layers related to a sharp decrease of the interlayer coupling strength as θ increases. Besides the presence of coexisting van Hove singularities, moiré superlattice Dirac points, anisotropic gaps at MB and anisotropic V f renormalizations, we observe non-dispersing band structure regions in strongly coupled graphene bilayers with a small θ (< 2.7 o ) providing a possibility to transform the bilayers into quasi-one-dimensional conductors by fine tuning of E f position.

Results and Discussion
Low Energy Electron Diffraction patterns (Fig. 1a inset) taken from the samples clearly show multiple and elongated diffraction spots reflecting various rotations of graphene domains, although it can be seen that there are also some preferential orientations, for example at 30 o with respect to SiC lattice vectors. This rotational misalignment can be directly visualized in the real space using scanning photoelectron microscopy imaging taken with the electron analyzer oriented towards a selected K point 24 . On Fig. 1a various rotational domains of regular shape can be observed. In the following, using angle resolved three-dimensional band structure mapping, we determine θ between the layers, the number of layers, velocities of the charge carriers and other important characteristics in individual nLG domains (see Methods).
Angle integrated spectra taken from different domains show saddle type VHS peaks shifting towards E f as θ decreases (Fig. 1b), similar to the observed by STS 11,12,[25][26][27] but only for the occupied electronic states contributing to the photoemission signal. Our results for the positions of the occupied VHS peaks summarized in Fig. 1c are in agreement with theoretical predictions and experimental results of ref. 11 and ref. 12, while some discrepancy can be found if our data are compared with those in ref. 25 and ref. 26. However, in the latter works there is noticeable asymmetry of VHS states with respect to E f . This may be related to graphene doping and also effects of the substrate or STM tip, whereas nLG domains in the present work are decoupled from the substrate and thus have charge neutrality level or Dirac K-points at E f . Note that in the spectrum with θ = 2.7 o ,−2.1 o the peak is double corresponding to two coexisting VHSs in a 3LG domain with two different rotations. This spectrum can be decomposed into two peaks after subtraction of a linear background with the positions of corresponding VHS peaks shown on Fig. 1b as two points with the smallest angles in the plot. We also report an angle integrated spectrum for a 4LG with one layer rotated by 7.2 o over 3LG with Bernal stacking showing three peaks discussed in detail in the last section.
Observation of nLG minizone and main features of electronic structure related to interlayer coupling. A large real space superlattice formed by a pair of twisted layers can be transformed into a minizone 14,16,[27][28][29] with K and K' points originating respectively from the main K and K' points of the two layers (Fig. 1d). We map k-space region around the coloured Dirac points in Fig. 1d, which, as shown in the following, can be easily identified by ARPES.
If n ≥ 3, and there are several periodicities when some of the layers experience two different concomitant couplings, the description in terms of a common minizone appears complicated. As we demonstrate below, in many cases the electrons in the superlattices corresponding to different periodicities can be treated independently because of big differences in the coupling strength. Constant energy photoelectron intensity distribution (k-map) at E = E f (I(k x , k y , E = E f )) in Fig. 2a shows the main K, K' and superlattice K' Dirac points from a tri-layer domain with the top layer (K) rotated by 5.   Fig. 2d. The intensity drop from K to K' due to attenuation of photoelectrons originating from the bottom layers is used here and below to discriminate the photoemission-related features of different layers, while further decrease from K' to K' can be rationalized considering that the latter Dirac points are formed as a result of a weak interlayer coupling. K points not clearly visible in Fig. 2a are also present but with much lower intensity concluding the formation of the minizone of the top and the middle layers. One can observe faint signatures of K points in Dirac spectrum where K point is expected (Fig. 2i). We speculate that such difference in photoemission between K and K' is due to a final state effect in photoemission process, in which the photoelectrons from lower layers experience stronger scattering by the top layer atoms. As a result, experimental k-maps follow triangular symmetry and theoretical band structure calculations, having hexagonal representations, are expected to have double number of features if compared to our experimental data.
First, we discuss in detail the electrons in the top layer of the above 3LG and show that all general features of the electronic structure predicted by tight binding calculations and related to the interlayer coupling and 'spinor' wave function of Dirac fermions can be clearly resolved: (I) The superlattice gap with the maximum opening of ~175 meV (at blue dots in Fig. 2e) is non-uniform along MB and its evolution can be followed in detail in Fig. 2f. The bands manifest 2/3π rotational symmetry around K as shown in Fig. 2e with zero gap along ΓK, while clearly open in KM (Fig. 2h); (II) The minizone has additional superlattice zero gap points below E f . For θ = 5.6 o such points are found in the middle between K and K', K and K' etc., (white dots in Fig. 2a) representing half of the theoretically predicted points for 'even' type superlattice 20 . We observe also 'odd' type superlattices, e.g. for θ = 4.3 o , where the points in the middle between K and K', K and K' correspond to the maximum MB gap opening (Fig. 3a), making characteristic difference between these two types of superlattices with the gap between K and K' present for the 'odd' (Fig. 3a) and absent for the 'even' (Fig. 3b) structures; (III) V f of the charge carriers in the top layer is renormalized from ~1.1 × 10 6 m/s of isolated graphene to 0.71 × 10 6 m s −1 (Fig. 4b) in quantitative agreement with tight binding model predictions 13 for a bilayer with θ = 5.6 o .
Charge carrier velocity renormalization vs twist angle. Dirac spectrum of a layer rotated by a large angle (Fig. 4a) is similar to the spectrum of isolated graphene with V f = 1.07 × 10 6 m s −1 as measured by ARPES 10,30 and Landau level STS in magnetic field by others 14,31 . For θ < ~10 o the velocity decrease becomes appreciable and can be followed comparing the dispersions in Fig. 4. The absence of the photoemission signal along KM for E < E f − 350 meV (Fig. 4b) and in both KM and ΓK for E < E f − 180 meV (Fig. 4c) is due to the gaps in corresponding MBs. V f ΓK in the top layer versus θ is summarized in Fig. 4d. Despite the data are obtained from nLG domains with different n, the results fit with good agreement to tight binding model prediction for a bilayer, with the best fit of the equation 2 (1), for t = 95 meV, comparable to calculated t = 110 meV, where t ≈ 0.4t and t is the interlayer coupling for Bernal stacking 13 . These results are in general agreement with V f deduced from STS measurements 14 and theoretical calculations 15 , however with slightly higher velocity renormalization. We tentatively explain this small discrepancy as an effect of electron-phonon coupling leading to observation additional velocity renormalization in photoemission 32   the layers near E f is manifested by the Fermi surface at K' (Fig. 2a) elongated in k x , making difficult to calculate θ using the procedure described in Methods. We estimated the twist angle by discerning the two shifted 'triangular' profiles of Dirac cones of the two bottom layers on k-map well below E f (dashed lines in Fig. 2d), where the interaction between Dirac cones of the layers is not pronounced. In accord with tight binding calculations of the electronic bands, predicting localization and hence flattening of the of electronic bands with the energy close to E f for a bilayer with small θ 15,16,28,33 , the bands at K' and, accordingly K', have weak dispersion (shown in the box of Fig. 2g).
As can be seen in other nLG domains, an equilateral triangle representing the Fermi surface of the underlying layers, surrounding the main Dirac point K of the top layer, is a common feature of all superlattices observed by us for different domains. The only exception is the discussed here bilayer, in which the supelattice Dirac points of the bottom layer (K" indicated as black circles on Fig. 2a) are absent at E f . We tentatively attribute this observation to Fermi surface instability related to charge localization or CDW observed by STS in a bilayer with small θ 11 . In fact, CDW instability is expected to remove spectral weight at the band regions around saddle point of VHS and to leave zero gaps, where the Fermi surface is well defined [34][35][36] . The minizone of the bilayer is too small compared to our k-resolution to deduce CDW Q-vector in detail and to make any temperature dependent systematic study, but the distortion of the bands near E f observed at the temperature of 110 K used in the experiment can be clearly observed in Fig. 2: the electronic structure of the bilayer gradually develops from non interacting bands (Fig. 2d) through gapped and zero-gap superlattice points at MB (Fig. 2c) into isotropic region (Fig. 2b) followed by elongated Fermi surface (Fig. 2a) indicating that CDW gap is ~20 meV.
Beside the discussed above phenomenon, the relief of electronic bands at a saddle type VHS is expected to result in the appearance of K-E plains with non-dispersing portions of the bands 17 . We observe these regions at E ~ E f − 40 meV within ~0.05 Å −1 in k-space, (Fig. 2i, between small dotted lines in K-E plain indicated in Fig. 2c), which can be more extended for greater twist angles as shown and discussed also in the following section (Fig. 5f). These non-dispersing band regions allow realization of peculiar transport phenomenon, when a localized hole-like perturbation, 40 meV below E f for the present bilayer with θ ~ 1 o , would propagate with a group velocity V = ћ −1 ∇ k E perpendicular to the flat band in the direction indicated by the arrow and also along its symmetrically equivalent directions. More interestingly, such flat bands can be adjusted to E f by uniform doping of the domain and cause a transition of the graphene bilayer from two-dimensional into a quasi one-dimensional conductor.
Coupling strength vs twist angle and interplay of several concomitant couplings for n ≥ 3. In  of slow or even localized electrons in the middle and bottom layers. Summarizing this result we can conclude that on the different 'floors' of twisted nLG domain the charge carriers can propagate with different group velocities and, if a given layer is found between the two layers with different twist angles, smaller rotation has bigger influence on the electronic structure of the layer. The difference in θ of the middle to the bottom and the middle to the top layers for the discussed above 3LG is big, however this rule holds when the difference in the twist angles is smaller. In the following example of another 3LG (Fig. 5), despite the absolute values of the rotation angles of the middle layer (K') with respect to the top layer (K, θ = 2.7 o ) and to the bottom layer (K", θ = − 2.1 o ), are similar, there is considerable difference in the corresponding VHS peak positions at E ~ E f − 150 meV and E ~ E f − 73 meV respectively (Fig. 1b). The group velocities of the charge carriers in the layers are also different with V f ≈ 0.2 × 10 6 m/s for the bottom and the middle layers, whereas V f ≈ 0.35 × 10 6 m/s for the top layer (Fig. 5d,e). The behaviour of the electrons in the middle layer is essentially more affected by the interaction with the bottom layer than with the top layer rotated by greater angle with respect to it. Indeed, the electronic band along K' -K" branWch is symmetric with similar V f in both the middle and the bottom layers (Fig. 5e). This result indicates that apparently, a minute change in θ dramatically affects the interlayer coupling. In fact, the interlayer tunnelling term t(Δ k) in tight binding calculations rapidly drops to zero with the increase of transfer momentum Δ K 16,[37][38][39] . In Fig. 5c we plot the ratio of the photoelectron intensity at E f for the superlattice Dirac point K' to the intensity at the main Dirac point K' I(k = K', E = E f )/I(k = K', E = E f ) that can be used as a quantitative measure of the relative coupling strength. Its quick decrease as Δ K increases illustrates the sharp dependence of the interlayer coupling on θ. However, the ratio drops for small angles as shown for Δ K = 0.032 Å −1 representing ~1 0 twist from the previous section. This drop in the ratio can be associated with the discussed above CDW instability. We obtained I(k = K', E = E f )/I(k = K', E = E f ) < 0.23 by fitting I(k, E = E f ) of Fig. 2a in the region around K' with four 2d gaussians representing two main (red and black dots in Fig. 2a) and two 'expected' (black circles in Fig. 2a) superlattice Dirac points situated in a similar fashion as K, K', K' of the big minizone. However, the actual ratio can be smaller or even zero since the experimental energy resolution is comparable to CDW gap.
Also in the minizone formed by the middle and the bottom layers of this domain we observe flat non dispersing bands (Fig. 5f). The group velocity of the excitation close to corresponding VHS peak, i.e. at E ~ E f − 60 meV, is directed along the red arrow perpendicularly to K' -K" and its equivalent directions (dotted triangle in Fig. 5b) and can be estimated ~1.5 × 10 5 m s −1 from the corresponding K-E dispersion in Fig. 5g. We could find clear evidence of the flat bands, which may be used to realize one-dimensional transport, in the minizones corresponding only to small θ such as 1.1 o and 2.1 o presented in the paper. In both cases E f shift needed for switching to quasi-one-dimensional conductivity is less than 100 meV and can be achieved by small doping level below 0.003 e − per unit cell, not changing bilayer electronic structure if each layer is doped by the same amount 8 .
Finally, we note that stronger coupling delocalizes the electronic bands between the layers and makes another twisted weakly coupled layer to interact with such system as a whole. To demonstrate this interaction we present an extreme case: a layer rotated by 7.2 o on top of 3LG with Bernal stacking (Fig. 6), representing the strongest interlayer coupling in nLG. Near E f , the electronic structure of the underlying tri-layer is composed of three common electronic bands, two reaching E f and one parabolic with a maximum at E = E f − 0.5 eV. Having odd number of layers, the carriers in one of the two bands at E f are massless 5 . These bands are described with good agreement by tight binding Hamiltonian and overlap parameters taken from ref. 40. Good correspondence suggests that electronic structure of the 3LG is not perturbed by the presence of the fourth rotated layer on top, while the top layer on the contrary has clear signatures of its interaction with the underlying 3LG domain. Indeed, there are three gaps in the minizone of the top layer (Fig. 6b,c) and three VHS peaks in the angle integrated spectrum (Fig. 1b) corresponding to each of the three electronic bands of the underlying tri-layer. Also V f at K is renormalized to 0.85 × 10 6 m s −1 . Interestingly, the data taken from the same top layer but from its part overlaying a 1LG with the same twist are presented in Fig. 3b. In this case there is only one VHS peak (Fig. 1b) and V f of the underlying layer in Fig. 3b is also renormalized, which is not the case in the massless carriers branch of the bottom tri-layer in Fig. 6. The superlattice potential created by 3LG is apparently different from a single layer showing no gap opening in K-K' direction characteristic for even type superlattice. The lower energy VHS position for the layer rotated by 7.2 o over 3LG is closer to E f than VHS position for the same rotation over 1LG by ~30 meV (Fig. 1b), in qualitative agreement with the results of comparison for a layer rotated over 1LG and 2LG with Bernal stacking 27 . However, the difference in VHS peak energies is too small to expect appreciable variation in V f in the top layer in these two cases.
In conclusion, the rich variety of electronic structure phenomena giving a possibility of a wide range tuning of the electronic structure of twisted few layer graphene is observed. In the presence of several interlayer couplings the electronic properties can be completely different in different layers of the same domain providing independent channels of electron transport. The group velocity of the charge carriers propagating in the moiré potential of strongly coupled (slightly twisted with θ < 2.5 o ) layers is not only renormalized, but can be anisotropic, giving unique possibility to transform the graphene from the two-dimensional into a quasi-one-dimensional conductor using two tuning parameters: E f (by doping) and θ. In the case of smaller twist angles (~1 o ), when the electron localization is expected, the electronic bands become weakly dispersing and hint on the appearance of CDW instabilities. We believe that the obtained results can further stimulate the device development, in which different layers may have separate electric contacts and the electronic properties of the graphene 'floors' can be engineered by their respective twists.

Methods
The experiment was performed using 27 eV photon beam from synchrotron radiation undulator beamline 3.2L of Elettra light source producing linearly polarized light (along k x direction) focused to ~0.6 μ m spot with Schwarzschild objective 24 and incident at 45 o with respect to the sample surface. ARPES data were acquired with a hemispherical electron energy analyzer registering the photoelectrons with total energy resolution of ~30 meV and angular resolutions of ± 0.15 o along the angular dispersion direction of the analyzer (k y direction) and ± 0.33 o across (k x direction). The sample was mounted on a scanning stage, which enables positioning and raster imaging with respect to the fixed photon beam. 3d photoelectron intensity distribution maps I(k x ,k y ,E) from the microscopic areas were taken by rotating the electron energy analyzer with respect to the sample, using a two-axis goniometer.
The samples, n-type (5 × 10 16 cm −3 ) C-face 6H-SiC single crystal (0, 0, 0, 1) oriented surfaces polished to < 1 Å RMS (Novasic) were first resistively heated in vacuum up to 1000 o C to clean the substrates and finally annealed in 1 bar of Ar at 1350-1500 o C for several minutes to obtain nLG flakes. nLG was prepared in a separate chamber directly connected to the measurement station and subsequently transferred for the measurements without vacuum loss. Whereas the graphene growth on Si-face SiC is epitaxial, for C-face in Ar atmosphere few micron sized domains of various rotational compositions between the layers are grown 41 . In accord with the scope of the present research, we were searching and comparing the results for the domains with various twists between the layers, while also the domains with Bernal stacking can be found as in ref. 42.
Dirac points, i.e. the crossing energy positions of Dirac conical spectra in graphene, are found at hexagon K vertexes of two-dimensional Brillouin zone readily accessible with ARPES 10 , and θ between the adjacent layers can be evaluated in the reciprocal space directly from Δ K = K-K' (Fig. 1d). |ΔK| = 2|ΓK|Sin(θ/2) and ΓK is the principal vector of graphene Brillouin zone (|ΓK| ≈ 1.7 Å −1 ). The number of layers in a particular domain is equal to the number of main non-equivalent K points (K, K', K" as in Figs 2, 3 and 5) or, for Bernal stacking, is inferred by the number of electronic π -bands 5,8,40,42 . For example, for the domain of 4LG in Fig. 6 there are two main K points and to one of them correspond three π -bands.
In order to select the rotational domains, the images (Fig. 1a) were taken by registering the photoelectrons at E f within ± 250 meV and ± 5 o angle windows 24 keeping the analyzer at 43 o with respect to the sample normal corresponding to photoemission from |ΓK|. In this case, the domains with the top layer oriented with K-point in the same direction as the analyzer give stronger contribution to the intensity, whereas the gradations of intensity on the image reflect the relative rotational orientation of the domains.
The constant energy photoelectron intensity distribution surfaces (k-maps) were obtained from I(k x ,k y ,E) within the energy window of ± 10 meV. V f of the charge carriers was estimated using linear dispersion relation for Dirac fermions E = ћV f K, where K is the absolute value of the wave vector measured with respect to Dirac point.