The scaling laws of edge vs. bulk interlayer conduction in mesoscale twisted graphitic interfaces

The unusual electronic properties of edges in graphene-based systems originate from the pseudospinorial character of their electronic wavefunctions associated with their non-trivial topological structure. This is manifested by the appearance of pronounced zero-energy electronic states localized at the material zigzag edges that are expected to have a significant contribution to the interlayer transport in such systems. In this work, we utilize a unique experimental setup and electronic transport calculations to quantitatively distinguish between edge and bulk transport, showing that their relative contribution strongly depends on the angular stacking configuration and interlayer potential. Furthermore, we find that, despite of the strong localization of edge state around the circumference of the contact, edge transport in incommensurate interfaces can dominate up to contact diameters of the order of 2 μm, even in the presence of edge disorder. The intricate interplay between edge and bulk transport contributions revealed in the present study may have profound consequences on practical applications of nanoscale twisted graphene-based electronics.


Lateral Force Measurements
The lateral shear force was measured during the sliding process in order to verify that sliding is performed under superlubric conditions, thus ensuring the existence of an angular mismatch at the bilayer graphene interface 1 . The shear force was evaluated using the relation = 2 , applicable for small shear distances 1 , where is the mesa radius and = 0.227 [J‧m -2 ] is the adhesion energy of graphite 1  We note that when shear forces are applied to bare graphitic surfaces, peeling and rippling effects may appear.
Nevertheless, the situation in our experiments is quite different, as the sliding interface is buried deep inside the pillar and is supported by two thick graphitic slabs. This induces very strict constraints on the out-of-plane motion of the carbon atoms. Furthermore, as discussed above, our sliding interfaces are incommensurate, thus the shear forces are expected to be very small and hence rippling and buckling are highly unlikely to occur. One may think that some rippling may occur at the exposed surfaces during the sliding process. Nevertheless, we do not see any evidence for such rippling in the frictional behavior of the system, that demonstrates power law scaling of the friction forces with the contact area with an exponent of 0.3 as expected for incommensurate flat circular contacts 1,2 . This is also supported by the fact that there is no experimental evidence of wear in our contacts (Supplementary Figure 1b). Furthermore, the bottom panel of Fig. 1A (see main text) shows that the measured current is symmetric with respect to the fullyeclipsed configuration. Namely, for a given absolute shift value the current is the same when the top mesa shifts towards the center or away from it. If, in the former case, puckering of the lower surface would occur in front of the sliding surface 3 , the current would not be symmetric. Another experimental support for these claims comes from the fact that all current vs. shift distance scaling laws derived in this paper rely on the circular geometry of the interface.
Random rippling and buckling effect would not obey such straight-forward scaling laws with the shift distance. Based on all the above, we may exclude surface rippling effects and rely on the rigid sliding interface model.

Fitting Procedure to Extract Interfacial Bulk and Edge Resistance
A numerical fitting procedure was employed to obtain the current vs. sliding distance profiles, ( ), and to extract the corresponding resistance of the sheared interface, , and its Edge (E) and Bulk (B) contributions. ( ) was calculated based on the equivalent electrical circuit depicted in Fig. 1d

Experimental Data of Current vs. Voltage on a Single Graphitic Contact
In the main text, we presented surface plots of the current vs. voltage and interfacial shearing distance (Fig. 2a). To construct these diagrams using the same contact, we performed the experiment by shearing the interface in steps

Electronic Transport Calculations
To evaluate the interlayer transport behavior of the graphene junction (modeled herein as a bilayer graphene system composed of two finite hexagonal or circular graphene flakes) as a function of the layers relative position and bias voltage we adopt the approach presented in Ref. 4. The current-voltage characteristics are evaluated via the Landauer scattering formalism 5 under which the current, is related to the transmittance probability through the system, ( ). The transmittance probability ( ) appearing in Supplementary Eq. (1), is calculated using the nonequilibrium Green's function technique for elastic electronic transport: Here, / ( ) are the retarded ( ) and advanced ( ) Green's function matrix representations of the junction and / ( ) are the top and bottom broadening matrices. The formers are given by: and where is a unit matrix of dimensions of the system and is the matrix representation of the device's Hamiltonian given by the following block form: and / = ( / ) † . Finally, the broadening matrices are given in terms of the self-energies as: The broadening, ℏ , is chosen to be sufficiently large to obtain a smooth density of states of the top and bottom flakes to mimic their periodic counterparts. The results are tested to be insensitive to this choice (see Supplementary note 5).
For computational efficiency we transform the transmittance probability expression of Supplementary Eq.
(3) to the diagonal basis of the dressed Hamiltonian: To this end, we denote by the transformation matrix that transforms the complex symmetric matrix r to its diagonal representation ̃: By inserting −1 or its conjugate transpose between each pair of matrices in Supplementary Eq. (3) and using the cyclic property of the trace operation we obtain: Note that, in the new basis, the retarded and advanced Green's functions matrices still obey the required relation: and have diagonal representations: This allows us, at the expense of a single complex symmetric matrix diagonalization and a single evaluation of ̃ and ̃, to evaluate ̃( ) and ̃( ) at any value of while avoiding matrix inversion. Note also that, since r is complex symmetric, is complex orthogonal such that = −1 (or * = ( −1 ) † ) and † = ( −1 ) * [2]. Hence, we have: It is important to note that the four transmittance probability contributions presented above ( → , → , → , and → ) only specify where electrons enter and exit the bilayer structure regardless of the specific route that they take when crossing the interface itself. This latter information is embedded in the Green's functions appearing in the trace formula (Supplementary Eq. (3)) and cannot be rigorously separated into bulk and edge contributions when there is coupling between the two regions that leads to delocalization of the wave functions over the entire flake surface.
We comment that the experimental values are analyzed under the assumption that the edge and the bulk contributions can be treated as parallel resistors. In practice, this means that the coupling between the edge and bulk regions can be neglected or that the wave packet of the transmitting electrons is sufficiently localized to pass through a specific surface region. These assumptions become valid with increasing junction lateral dimensions.

Validity Test for the Value of the Broadening Factors
Within our approximate treatment of the electronic transport problem the coupling of the modeled interface to implicit particle reservoirs is done by introducing broadening factors, ℏ , which mimic the imaginary part of the implicit bath self-energy within the wide-band approximation. This introduces a finite lifetime to the atomic states of the graphene flakes to which the broadening factors are applied. In energy space this translates to broadening of the -function eigenstates of the isolated flake into Lorentzian functions of width ℏ . To avoid artefacts resulting from the discreteness of the spectrum of the finite model flake, the value of the broadening factors should be sufficiently large to result in a continuous density of states (DOS). Care should be taken also not to use too large broadenings to allow for the characteristic electronic structure of the flakes to be manifested in the calculation.
To validate that our choice of broadening factors fulfils these requirements we plot the density of states of one of the flakes constructing the interface while assigning each eigenvalue, , a Lorentzian function of width ℏ such that the total DOS is calculated as follows:

Convergence Test with respect to the Number of Modeled Layers
The results presented in the main text were obtained using a bilayer model system with broadening factors applied directly to the atoms of the contacting layers. This allowed us to perform calculations on large-scale junctions with reasonable computational burden. However, contacting the interfacing layers directly to the implicit leads has the effect of exaggerated broadening of the transmittance probability curves. Therefore, to confirm the validity of our qualitative conclusions, we repeated some of the calculations using a fourlayer model system.
The Hamiltonian of the device can be written in block form as follows: where

Illustration of the Wave Functions of Eight Eigenstates Within the Fermi Window of Circular Junctions
To obtain the molecular orbital plots of the edge states presented in Fig. 3a-c of the main text we diagonalized the tight-binding Hamiltonian 6 of the relevant junction models and plotted the absolute squared molecular orbital expansion coefficients over the different atomic sites. In Fig. 3d of the main text we presented the angularly averaged molecular orbital weights as a function of distance from the flake edge averaged over eight molecular orbitals residing within the Fermi window. For completeness, we present in Supplementary Fig. 9 the individual eight molecular orbital absolute squared weights of the 20 nm circular bilayer system all showing pronounced edge character. The corresponding angularly averaged molecular orbitals weights as a function of distance from the flake edge are presented to the right of each molecular orbital for the 10 nm (blue circles), 15 nm (red circles) and 20 nm (grey circles) junctions.