Minibands in twisted bilayer graphene probed by magnetic focusing

Detection of magnetic focusing peaks allows probing of minibands in the energy spectrum of twisted bilayer graphene.


Section 1. Determining twist angle.
We used two independent methods to find the actual twist angle between the two graphene layers. First, the twist angles were found from xx ( ) and xy ( ) measurements using the position of the secondary neutrality points (NP): where is graphene's lattice constant, the carrier density, = 4/∆ the superlattice unit cell area, and ∆ the position of the secondary NP. For device D1, NP are shown in Figs 4A and 4B in the main text yielding Δ = 8.1 × 10 12 cm -2 and = 1.87°. For device D2 it was impossible to reach secondary NPs using electrostatic gating; to find the twist angle in this case we used the fact that NPs can be found by extrapolation of the reciprocal Hall resistivity, ( xy ) −1 as shown in Fig. S1A. The main neutrality point is clearly seen at zero carrier density and two van Hove singularities (vHS) are shown by red arrows. To find the secondary NPs, the linear dependences below and above the vHS were extrapolated as shown in Fig. S1A, with NPs corresponding to the intersections with

Section 2. Device fabrication and measurement details
The heterostructures studied in this work were assembled using the standard dry-transfer technique (24,25), and for the fabrication of the twisted bilayer graphene (TBG) we adapted the tear-and-stack (26,4) method. Details of these methods are outlined below.
First, the top hexagonal boron nitride (hBN) crystal was picked up using a polypropylene carbonate (PPC) polymer spun onto a polydimethylsiloxane (PDMS) film. Then we used a micromanipulator to place the hBN crystal so as to cover only a part of the monolayer graphene located on a SiO 2 /Si substrate. Next, hBN was slowly peeled off the substrate, tearing the graphene flake into two pieces while picking up the part covered with hBN. The remaining part of the graphene flake was rotated by 2˚ and picked up with the first half attached to hBN to produce TBG. The temperature of the substrate was kept at 70˚ C throughout this process in order to reduce thermally induced strain or relaxation of the layers. By carefully controlling the micromanipulator, we ensured that graphene layers had no contact with the PPC polymer, guaranteeing a clean interface between the two graphene monolayers. Finally, the bottom hBN crystal was picked up to encapsulate the TBG, and the whole stack released onto a SiO 2 /Si substrate. To define 1D contacts, we used reactive ion etching to selectively remove the heterostructures areas, followed by deposition of Cr (3 nm) and Au (60 nm). An additional lithography step was used to make a gold top gate, which also served as an etching mask to define the mesa. An example of the final device is shown in Fig. S2A, where we show one of our dual gated TBG samples. For this sample we used p-doped Si as the bottom gate and Au as the top gate.
Resistance measurements were carried out using the standard low-frequency lock-in technique with a small excitation current ~100 nA; this ensured negligible heating effects down to the lowest measurement temperature, T = 2K. The dual-gated geometry allowed us to control the total carrier density and the displacement field independently. The total carrier density is the sum of carrier densities induced by the top and bottom gates: total = 1 ( tg tg + bg bg ), where tg and bg are the top and bottom gate voltages, is the electron charge, tg and bg are the respective capacitances per unit area of the top and bottom gate (here tg and bg were obtained from the Hall measurements). The displacement field is calculated using = 1 2 0 ( tg tg − bg bg ), where 0 is the vacuum permittivity. To achieve a fixed displacement field, tg and bg were varied simultaneously according to the above formula, so that only the total carrier density changed. An example of such measurements is shown in Fig. S2B.

Section 3. Calculations of TBG band structure
To calculate the band structure of TBG, we used the model reported in ref. (31), where the Hamiltonian is given by: This is equivalent to the continuum-model Hamiltonian derived in ref.
(3) up to a gauge transformation. The on-diagonal blocks describe the top and bottom layers of graphene where ν is the Dirac velocity for the monolayer and Δ the energy shift induced by the perpendicular displacement field. The off-diagonal blocks describe the interactions between the two layers, with the interlayer coupling strength given by w=110 meV.
The original band touching point of the top layer is placed at the corner of the mini Brillouin zone (mBZ) κ, and the band touching point of the bottom layer is placed at κ'. The vectorsΔ j accounting for the shift between Brillouin zone corners for the two layers (as illustrated in Fig. S3) are given by: where is the twist angle between the two layers in radians and a=2.46 Å is graphene's lattice constant. The miniband spectrum is calculated by zone folding (32), i.e., by bringing the states in the K valley with momenta connected by the reciprocal lattice vectors of the moiré pattern to the first mBZ. The basis of k-states of the top and bottom layer are formed from the reciprocal lattice vectors ∆ 0 + m 1 + n 2 and 2∆ 0 + m 1 + n 2 , respectively, where 1 = ∆ 0 − ∆ 1 and 2 = ∆ 0 − ∆ 2 . The number of basis states is chosen to ensure convergence of the first three conduction and valence bands. The resulting Hamiltonian that contains the matrix elements between the basis states is diagonalised using a similar method to ref. (33).

Section 4. Numerical simulations of transverse magnetic focusing
To simulate TMF maps shown in Fig. 2 in the main text, we first calculate the band structures to extract the Fermi surfaces and then determine the cyclotron orbits in real space by rotating the Fermi surfaces by 90° and scaling by ħ/eB (the scaling factor is obtained from the quasi-classical equations of motion (32)). Charge carriers propagate either clockwise or anticlockwise, depending on the sign of the effective charge. We assume specular boundary conditions, so that in a magnetic field the carriers travel along the edge of the sample following cyclotron (skipping) orbits and caustics of skipping orbits focus onto equidistant points. The drift direction of the skipping orbits depends on the effective charge of the carriers and the direction of the magnetic field.
To achieve consistency with the experiment we select the states that are moving away from the injection point with energies between ε f and ε f +eV (where ε f is the Fermi energy and V the applied voltage). The group velocity is calculated from the band structure using v ∝ , i.e., it is related to the energy dispersion (for example, the velocity is smaller in flatter parts of the dispersion). Accordingly, as the applied voltage elevates the Fermi level, it results in extra states being occupied such that the available states are populated with a probability proportional to | | −1 and different injection angles should be weighted with a probability proportional to the density of states.
The TMF spectra shown in Fig. 2B, E in the main text are calculated numerically by using a similar method to ref. (12). This is achieved by counting how many electrons enter contact 3 (Fig. 1C in  . Here the subscripts correspond to the device contacts in Fig. 1C and d i is the distance between consecutive skips along the edge of the i th trajectory. To investigate whether the TMF spectra are sensitive to the crystallographic orientation of graphene layers with respect to the skipping direction (the edge of the sample), Fig. S4 compares TMF maps simulated for different edge orientations characterised by an angle Ф. To this end, we fix the orientation of one of the monolayers so that Ф = 0° corresponds to the zig zag edge and Ф = 90° to the armchair edge. The results for parameters of device D1 at |n| = 6.6x10 12 cm -2 give triangular skipping orbits with the distance between the focusing peaks along the sample boundary weakly dependent on Ф -see Fig. S4. Similar results are obtained for all carrier densities where the Fermi surfaces are anisotropic, i.e., for |n| > 3x10 12 cm -2 where the Fermi surface around the  point has a pronounced triangular shape (see Fig. 2C in the main text). Corresponding TMF maps show focusing peaks at slightly shifted positions relative to each other. The 3-fold symmetry of the triangular Fermi surface means that the TMF maps should repeat after 60° as is indeed seen for Ф = 30° and 90° in Fig. S4, where the results are identical. At low carrier densities, near the main neutrality point, the Fermi surfaces are almost isotropic and the TMF maps are independent of Ф. The positions of vHSs are independent of Ф as well, in agreement with ref. (12). Fig. S4. TMF maps and simulated skipping orbits for different edge alignment. The TMF maps are simulated for Ф = °,°,°,° and° for device D1. The orientation of one of the monolayers is fixed such that ° corresponds to the zig zag edge and ° to the armchair edge. The skipping orbits are shown at |n| = 6.6x10 12 cm -2 .

Section 5. Electrostatic screening in a finite displacement field
The TMF map in Fig. 2E in the main text shows the effect of a finite displacement field between the two graphene monolayers. To find the effective electric field for each n in this figure, we need to take into account electrostatic screening. At twist angles ~2° and low carrier densities, the two monolayers are almost decoupled. To take into account electrostatic screening in this case, we include a screening term as proposed in ref. (34,35): where n 1 and n 2 are carrier densities in the two parallel graphene layers separated by a distance d, D the applied displacement field, the Dirac velocity, the band indices s 1 and s 2 are given by s i =n i /|n i |, and the electron charge e < 0. In case of the TBG, we use ≈ 0.34 nm and following Refs. (35,36,37), the dielectric constant for twisted bilayer = 2.7. The total carrier density n is given by Eq. (S5).
To find the effective electric field for each value of n and D used in the experiment, the two equations are solved simultaneously using the Dirac velocity for monolayer graphene, v=10 6 m/s. In our calculations, we take n 1 to be the bottom layer and n 2 to be the top layer. The positive direction of D is from the top to the bottom (pointing downwards).  Fig. 1a in the main text), where is the distance from the current injector to the voltage probe, which is independent of the angle between the crystallographic axes orientation of graphene layers and the sample edge. Near the secondary neutrality point ℎ is sensitive to the relative orientation of the graphene layers and the sample edge, as can be seen in Fig. S4. We have calculated ℎ and extracted the corresponding scattering lengths for several characteristic angles between the device edge and zig zag axis of the top graphene layer: ℎ0°= 1.11 ; ℎ10°= 1.25 ; ℎ20°= 1.58 ; ℎ30°= 1.77 ; ℎ45°= 1.93 . Using these values we extracted the scattering lengths, presuming different crystallographic orientations and compare them in Fig. S6. This showed that in all cases scattering lengths vary between ~100 µm at low T and a few µm at T =30K, indicating the importance of electron-electron scattering at elevated T as discussed in the main text.  S6. Temperature dependence of electron scattering length. Electron scattering lengths corresponding to different relative orientations of the graphene's crystallographic axes and the sample edge were extracted from the temperature dependence of the first focusing peak in Fig. 3B of the main text using ℎ calculated as described in Supplementary section 6.