In situ twistronics of van der Waals heterostructures

Yaping Yang, Jidong Li, Jun Yin, Shuigang Xu, Ciaran Mullan, Takashi Taniguchi, Kenji Watanabe, Andre K. Geim, Konstantin S. Novoselov, Artem Mishchenko School of Physics and Astronomy, University of Manchester, Oxford Road, Manchester, M13 9PL, UK National Graphene Institute, University of Manchester, Oxford Road, Manchester, M13 9PL, UK State Key Laboratory of Mechanics and Control of Mechanical Structures and MOE Key Laboratory for Intelligent Nano Materials and Devices, College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, 210016, China National Institute for Materials Science, 1-1 Namiki, Tsukuba, 305-0044, Japan Centre for Advanced 2D Materials, National University of Singapore, 117546, Singapore *e-mails: ypyang0916@gmail.com, artem.mishchenko@gmail.com

Graphene placed on hexagonal boron nitride (hBN) represents a prototypical system of a twisted heterobilayer. In this twisted system, when at a small twist angle, mono-or few-layer graphene in a quantizing magnetic field exhibits a fractal energy spectrum known as Hofstadter's butterfly [5][6][7] . This twisted system has been reported to bear topological bands and strong correlations, exhibiting a number of fascinating many-body phenomena, such as correlated insulating states and superconductivity reported in trilayer graphene on hBN 8,9 , and fractional Chern insulators observed in bilayer graphene on hBN 10 . In twisted homobilayers, strongly-correlated phenomena including correlated insulating states, unconventional superconductivity and ferromagnetism, have been observed in twisted bilayer graphene [11][12][13][14] and twisted bilayer-bilayer graphene systems [15][16][17][18] . A variety of other systems, such as twisted transition metal dichalcogenide layers [19][20][21][22] , graphene/WS2 bilayer 23 , and twisted bilayer graphene aligned to hBN 24 , have also manifested many intriguing phenomena.
In twistronics, accurate positioning, rotation, and manipulation of 2D materials are needed to fabricate a system with desired twist angles. To this end, several techniques have been developed, including optical alignment of crystal edges 7,25 , tear and stack technique for twisted homobilayers 26 , and in situ rotation mediated by atomic force microscope (AFM) tips 27 . Here we present a new experimental strategy to dynamically manipulate layered heterostructures in situ with precise control, allowing investigation of optical, mechanical, and electronic properties of a system with arbitrary twist angles between individual layers.
For our in situ twistronics technique, we use a glass slide with a droplet of polydimethylsiloxane (PDMS) as a manipulator, which is cured and naturally shaped into a hemisphere geometry. For a carefully fabricated PDMS hemisphere, the contact area between the manipulator and a 2D crystal can be as small as a few tens of micrometers 26,28 , which depends on the hemisphere radius and is highly sensitive to the contact force, making it difficult to precisely control the target flake. To solve this problem, we developed a critical step, where we intentionally deposit an epitaxial polymethyl methacrylate (PMMA) patch on top of a target flake through a standard electron-beam lithography (EBL) process. The concept is illustrated in Fig. 1. The epitaxial PMMA patch could be designed into an arbitrary shape that fits the target flake and is normally a few hundred nanometers thick, thereby the contact area between the PDMS hemisphere and the flake is precisely limited to the area of the epitaxial PMMA patch without PDMS touching beyond the target flake at a much larger contact force threshold.   S1). Figure 1e shows the top view when PDMS hemisphere touches PMMA patch, where the interference rings imply that PDMS hemisphere is in the proximity of the flake without touching it. The friction between two incommensurately stacked 2D materials is dramatically reduced thanks to the superlubricity 29,30 , and, therefore, the top hBN, as well as the underneath graphene, can not only slide but also rotate freely on the surface of the bottom hBN under the control of the PDMS hemisphere (Fig. 1c, d, f and g, and Supplementary Video 1 and 2). Figures 1f and g show that the PMMA patch is in good adhesion to the top hBN and does not delaminate during the manipulation process.
The manipulation technique presented here enables continuous modification of the twist angle between the layers even after the heterostructure is already assembled. In order to perform in situ optical measurements such as Raman spectroscopy after each manipulation step, the PMMA patch was intentionally designed not to cover the graphene so that it can give strong enough signal ( Fig. 1f and g). Compared to previous in situ rotation technique mediated by an AFM tip 27 , our manipulation technique is convenient and reproducible since the PMMA patch can be easily washed away by acetone and re-patterned by EBL. In addition, our technique can manipulate flakes regardless of their thickness, whereas an AFM tip might destroy thin flakes.
To demonstrate the potential of the manipulation technique in twistronics, we fabricated another hBN/graphene/hBN heterostructure (sample 1) where the graphene layer was aligned to both top and bottom hBN layers. We firstly assembled the heterostructure intentionally with 0 < t, b < tb < 60⁰ (Fig. 2d), where t, b and tb are twist angles between graphene and top hBN, graphene and bottom hBN, as well as top hBN and bottom hBN, as shown in Fig. 2a and b. Then we used our technique to rotate the top hBN layer, by making sure that the rotating PMMA patch is patterned within the top layer (Fig. 2d). The twist angles t and b varied simultaneously as the rotation progressed until either of them reached 0⁰ where the graphene layer was locked to one of the hBN layers, that is graphene/hBN interface went through a transition from incommensurate to commensurate state 31,32 . At a certain moment the smooth rotation stacks and the PMMA delaminated from the top hBN layer (Fig. 2e). We associate such conditions with t and b were both ≈0⁰, indicating that graphene was aligned to both hBN layers. At this stage, the two hBN layers were aligned to each other as well, with tb also reaching 0⁰, as shown in Fig. 2e. The delamination of PMMA after all the 2D layers are locked to each other indicates that the interaction between the 2D materials at commensurate state is stronger than that between the PMMA and the top 2D crystal (see section 2 in SI).
In principle, there are two types of alignment in a doubly-aligned hBN/graphene/hBN heterostructure, that is t=b=tb=0⁰ and t=0⁰, b=tb=60⁰. The state of the resulting stack is determined by the initial settings of t and b and the rotation direction. Scenarios for various t and b are presented in Supplementary Fig. S2, where we also discuss in details the motion of incommensurately stacked 2D materials (section 2 in SI). The initial settings of the twist angles rely on the known crystal orientation of graphene and hBN layers, which is normally distinguished by straight long edges of the flake. Since hBN crystals are three-fold rotationally symmetric, using hBN flakes from the same exfoliated crystal can help to secure the crystal orientations.
However, for hBN, the symmetry between the top and bottom atomic layers depends on whether the number of layers in hBN is odd or even (Fig. 2c). This ambiguity can, in principle, be circumvented by using the same surface, either top or bottom, of the original hBN crystal for alignment. To demonstrate how this can be done we fabricated sample 2 with t=0⁰, b=tb=60⁰ in the final stack (see section 2 and Figs. S3, S4 in SI). To confirm the alignment of graphene to both hBN layers, we carried out Raman characterization, as shown in Fig. 2g and h. The graphene layer in sample 1 originally contained monolayer and bilayer regions ( Fig. 2f). For the monolayer region, the full width at half maximum of 2D peak (FWHM2D) increases from 17 cm -1 to 55 cm -1 after rotation; G peak width also broadens with the emergence of lower frequency components (Fig. 2g). These results are consistent with previous reports 33,34 . The broadening of 2D peak by near 40 cm -1 and the appearance of lower frequency components in G peak signify near-perfect alignment between the hBN and graphene crystals and arise from graphene coupling to the moiré potentials from both top and bottom hBN crystals 27,34 . For the bilayer region, the splitting of 2D peak is strongly enhanced when graphene is misaligned to both hBN layers, whereas after rotation, these components broaden and shift towards each other, resulting in a prominent change in the line shape of 2D peak. In contrast, G peak splits into two distinct components after rotation. These results for doubly-aligned bilayer graphene have not been reported before and we ascribe the change in Raman spectra to the periodic strain field induced by moiré patterns (see section 3 in SI). Figure 2h shows the 2D peak width Raman map of the graphene layer after rotation. It clearly demonstrates the homogeneous distribution of FWHM2D among different regions of the graphene layer, indicating a spatially uniform twist angle in the stack. The AFM topography (Fig. 2f) here shows that the rotation process did not damage or crease the graphene layer.
To further quantify the existence of two moiré superlattices at both sides of graphene layer in this heterostructure, we made it into a device and investigated its transport properties. We will now focus on our bilayer device, Fig. 3. Figure 3a shows the longitudinal resistivity ρxx and transverse resistivity ρxy as a function of charge carrier density n at non-quantizing magnetic field of B = 0.03 T. We observed both broadened resistivity peak at the primary Dirac point (PDP) and satellite resistivity peaks at finite densities situated symmetrically with respect to PDP. The Hall resistivity ρxy near each satellite peak changes sign, suggesting that they are moiré-induced secondary Dirac points (SDPs).
In quantizing magnetic fields, the presence of moiré potential splits the Landau levels into a fractal structure: the Hofstadter minibands separated by a hierarchy of self-similar minigaps, to which the In the fractal superlattice spectra, the self-similar minigaps occur at  = 0p/q, where  = BA is the magnetic flux per superlattice unit cell, A is the superlattice unit cell area, 0 is the magnetic flux quantum, p and q are co-prime integers. The strongest minigaps arise at  = 0/q, resulting in Brown-Zak (BZ) magneto-oscillations 36 with the periodicity of 1/B = A/0. These minigaps also correspond to the intersections between the central and satellite fans in the Landau fan diagram 5,6 . We observe Brown-Zak oscillations (horizontal streaks visible in Fig. 3b) belonging to the two sets of satellite fans, as shown by the horizontal blue and magenta dashed lines in Fig. 3d, which provides independent confirmation of the two distinct moiré periods.
Brown-Zak oscillations are more clearly seen in ∂/∂B(B,n) maps in Fig. S8b, e and f, in which we observed other p/q fractions. High temperature suppresses cyclotron oscillations, making the BZ oscillations more visible, Fig. S9. At electron doping where BZ oscillations are more pronounced 36 , two sets of maxima, corresponding to the two periodicities, 15.2 nm and 14.7 nm, are clearly seen. This independently verifies the existence of two moiré superlattices with different wavelengths at both sides of the graphene layer. 0.0510 12 cm −2 (right inset of Fig. 3a). Notably, in the plot of ∂xx/∂B(n,B) within |n| < ±210 12 cm −2 we observe a prominent horizontal streak at B = 0.45 T (Fig. 3c), which perfectly matches the magnetic field of the first-order magnetic Bloch state originating from the composite moiré pattern when BA = 0 (where q = 1). The results for the bilayer graphene region discussed above are similar to what was found in the monolayer graphene region (Figs. S10-S12). The manipulation technique presented here proves successful in making the hBN/graphene/hBN heterostructures with perfect alignment between graphene and the two hBN layers with a high twist angle homogeneity. In contrast, optical alignment of crystal edges or heating of the final stack cannot consistently guarantee perfect alignment with the twist angle < 1⁰ in a single stack 5-7,10,37,39 .
In conclusion, we introduce an in situ manipulation of van der Waals heterostructures mediated by patterning a polymer resist patch onto target flakes, which can precisely and dynamically control the rotation and positioning of 2D materials by a simple gel stamp manipulator. Using this technique, we realized a perfect alignment of graphene with respect to hBN layers with a much higher success rate compared to conventional optical alignment of crystal edges during micromechanical transfer. Our technique, which can be easily generalized to other 2D material systems (as shown in Fig. S5), opens up a new strategy in device engineering for twistronics, finding its applications in the design of 2D quasicrystals 42,43 , magic-angles flat bands 11-14 , devices with AB/BA domain walls 44 and other topologically nontrivial systems.

Van der Waals assembly
All heterostructures were assembled using standard dry-transfer technique 45 using a polymethyl methacrylate (PMMA) coated on a polydimethyl siloxane (PDMS) stamp, and a SiO2 (290 nm)/Si as a substrate.
The devices made from the monolayer and bilayer graphene regions are in a Hall bar geometry, with electrical contacts made by Cr/Au (3 nm/80 nm).

Raman characterization
The Raman spectra were acquired by Renishaw Raman system with 1800 lines/mm grating, using linearly polarized laser radiation at the wavelength of 532 nm. The laser power was kept below 5 mW. The Raman spatial maps were taken with a step size of 0.5 µm. To extract the peak width of 2D band, we fit the spectrum at each pixel in the spatial mapping to a single Lorentzian function.

Transport measurements
Transport measurements were carried out in a four-terminal geometry with a low-frequency ac current excitation of 100 nA using standard lock-in technique at the base temperature of 0.3 K (Brown-Zak oscillations were measured at 70 K). The devices were gated to control a charge carrier density (n) and a displacement

Calculation of moiré wavelength, twist angle and super-moiré wavelength
The moiré wavelength  is calculated using the geometric relation ns = 4/ √3 2  2 , where ns is the charge density at full filling of the moiré miniband (at the density of four electrons per superlattice unit cell). The relation between twist angle  and moiré superlattice wavelength  in hBN/graphene/hBN heterostructure is given by 46 : where a is the graphene lattice constant, δ is the lattice mismatch between hBN and graphene which is 1.65% (we used the value calculated by Ref 33,38 ), ϕ is the relative orientation of the two moiré patterns with respect to the graphene lattice (φ1 and φ2, respectively).
Using the method described previously 37 , for super-moiré wavelength Λ, in analogy to the equation above, the relation between Λ and the twist angle  between the constituent moiré patterns (the wavelengths being 1 and 2, 1≥ 2) is given by: where  is the mismatch between the constituent moiré patterns, given by: and the twist angle  is given by: The modular division used here to find  reproduces the output of the piecewise conditional reported in Ref

Identification of the initial contact between polymer gel and epitaxial polymer
In our manipulation technique, the critical step to precisely control the movement of 2D flakes is the identification of the initial contact between the polymer gel (PDMS) and the epitaxial polymer patch (PMMA).
The signature of the contact is a colour change in the PMMA patch and an appearance of a sharpen edge of the contact area between the PDMS gel and PMMA patch, which is easily distinguished under optical microscope, Fig. S1.

The motion of incommensurately stacked 2D materials
In our manipulation technique, the motion of 2D materials in the van der Waals heterostructure depends on the balance between external driving force and the adhesion between PMMA patch and the top 2D layer, as well as the kinetic friction originating from the atomic shear force between the adjacent 2D layers. The kinetic friction between 2D layers is dramatically reduced in the incommensurate state, the so called superlubricity where the 2D flakes can move smoothly. Inversely, in the commensurate state, the two adjacent 2D layers are locked together due to pinning of the boundaries that separate local regions of the commensurate phase 1-3 . Such boundaries, either topological defects or misfit dislocations, are caused by the lattice mismatch and twist angle between adjacent layers.
For graphene/hBN and hBN/hBN interfaces, the commensurate state occurs at a small twist angle between the layers, since these two interfaces already satisfy the condition of a small lattice mismatch (1.65% 4,5 for graphene/hBN interface and 0% for hBN/hBN interface). The locking of the graphene/hBN and hBN/hBN interfaces limits further motion of the flakes at commensurate state, thereby when the external driving force is large enough, the PMMA patch/hBN/graphene/hBN stack will break at the weaker coupled PMMA patch/hBN interface, in agreement with the delamination of PMMA patch when graphene is aligned to both hBN layers after rotation (Fig. 2e in the main text and Fig. S3c). Rotation of the top two layers of hBN/graphene/hBN heterostructure with different initial settings of t and b should result in the final stack with commensurate or incommensurate states at the two interfaces, as shown in Fig. S2 and Table S1. Here we consider both the 0⁰ and 60⁰ alignment of graphene and hBN layer, since these two types of alignment have distinct symmetry and affect electronic structure of graphene differently 4,6 .
We fabricated hBN/graphene/hBN heterostructure (sample 2) with the manipulation process matching the cases in Fig. S2b (initial alignment) and Fig. S2j (final alignment). To make sure that the crystal orientation of the flakes is known, we first identified hBN fakes which have fractured into two pieces during the mechanical exfoliation procedure (Fig. S3a). with t=0⁰, b=tb=60⁰ (Fig. S3c), which is the case of Fig. S2j.
We carried out Raman spectroscopy to prove the alignment of graphene and hBN in sample 2, as shown in Fig. S3d and e. The behaviour of Raman spectra in monolayer and bilayer graphene regions are similar to those observed in sample 1. Whereas the full width at half maximum of 2D peak (FWHM2D) increased from 17 cm -1 to 48.5 cm -1 after rotation, smaller than that in the perfect alignment case, which means that the twist angles t and b are larger than those in sample 1. Figure S3e shows the distribution of FWHM2D among the graphene flake after rotation, indicating a spatially uniform twist angle in both monolayer and bilayer regions. The AFM topography (Fig. S3f) here shows that the rotation process did not damage or crease the graphene layer.
The assembly process of sample 2 is illustrated in Fig. S4. The details of the process are as follows.
Step 1, we used PMMA coated PDMS block mounted on a glass slide to pick-up one of the fractured hBN pieces as the top hBN (hBN piece 1).
Step 2, we used polypropylene carbonate (PPC) coated PDMS block mounted on another glass slide to pick-up the other fractured hBN piece (hBN piece 2), which serves as the bottom hBN.
Step     For other interfaces where the adjacent layers have larger lattice mismatch, such as MoS2/hBN interface, a small or even zero twist angle will not result in the commensurate state. The kinetic friction between the 2D layers always remains at an extremely low value, thereby the 2D flake can move and rotate smoothly even when passing through the point where twist angle is zero without the delamination of PMMA patch, as shown in Fig. S5.

Raman spectra of graphene in the monolayer and bilayer region before and after rotation
Raman spectra of graphene are strongly modified when graphene is subjected to double moiré superlattices (Fig. 2g). In our sample, for monolayer graphene region, G peak changes its shape from a narrow peak centred at 1582 cm -1 with FWHMG ≈14.5 cm -1 to a broadened asymmetric peak after rotation. The presence of a low energy shoulder at ≈1558 cm -1 of G peak after rotation is attributed to a TO phonon, similar to what was observed in graphene aligned to only one hBN crystal 7 . For 2D peak, the broadening effect is strongly enhanced in double-moiré superlattice compared to a simple aligned graphene/hBN heterostructure: ≈20 cm -1 increase in FWHM for graphene in a single moiré superlattice 7,8 and ≈40 cm -1 increase in FWHM for our doubly-aligned sample, indicating that the double moiré superlattices induce a much stronger periodic inhomogeneity originated from charge accumulation, strain, etc.
For bilayer graphene region, G peak splits into two components in the presence of double moiré superlattices. For 2D peak, the four components broadened and their position differences reduced, which is similar to what was reported in bilayer graphene aligned to only one hBN crystal 9 .
We found an overall downshift of around 2cm -1 for both G peak and 2D peak among the whole flake after rotation (Fig. S6 a-d), which is contradictory to the results in previous study 4,8 . Similar behaviour is found in Sample 2 (Fig. S7). Near the corner of the folded region (as indicated by the arrows), before rotation, the peak positions of G and 2D peaks are slightly lower than those of the nearby region. We attribute this phonon softening to the strain caused by folding 10 . Whereas after rotation, peak position near the corner of the folded region shows an upshift behaviour, and is higher than that of the nearby region. Note that during the rotation, the folded graphene rotated as well, therefore we expect an enhanced strain effect near this region. However the abnormal behaviour of the peak position implies a mechanism beyond strain effect.
To further look into the line shape of G peak after rotation, we fit G peak with two Lorentzian peaks for both bilayer and monolayer graphene regions and plot the maps of the peak position difference and intensity ratio of the two components, as shown in Fig. S6e and f. The maps show that the G peak splitting is homogeneous for both bilayer and monolayer graphene regions, which is around 13 cm -1 and 22 cm -1 , respectively. The intensity ratio of the two components is around 1 for the bilayer region, which means that they have comparable weight, and is around 0.6 for the monolayer region.
The overall downshift and broadening of the peak positions for G and 2D peaks and splitting of G peak remind one of the effects of strain. Under uniaxial strain, the doubly degenerate E2g mode splits into two components, one along the strain direction and the other perpendicular, which leads to the splitting of G peak 11 . Tensile strain usually gives phonon softening, while compressive strain usually leads to phonon stiffening. In the presence of moiré superlattice, the strain spatial distribution is periodically modulated, and

Transport properties of graphene in double moiré superlattices.
In commensurate magnetic field and periodic moiré potential, the fractal energy spectrum of electron in graphene exhibits a self-similar recursive, where the superlattice minibands show a gapped Dirac-like spectrum at  = 0p/q and exhibit Landau levels (LLs) characteristic of Dirac fermions 15  Based on the Raman spectra of the hBN/graphene/hBN heterostructure, we estimate that the twist angles t and b should be both close to 0⁰, thereby in the transport properties, we expect to see two sets of SDPs close to each other apart from the PDP. For both bilayer and monolayer graphene devices, in longitudinal resistivity (ρxx) vs carrier density (n) dependence, we observed broadened (and with a complex structure) satellite peaks located around n = 210 12 cm −2 (Fig. 3a in the main text and Fig. S10a). At small B = 0.03 T where the Landau quantisation is not yet developed, the transversal resistivity ρxy changes sign in the carrier density regions of these satellite peaks, indicating that they are moiré superlattice induced SDPs. These SDPs are more prominent in the hole side compared to the electron side, which is consistent with previous studies 16,17 .
In the Landau fan diagram at high magnetic fields, by tracing the linear trajectories we observed two sets of SDPs induced by the two independent moiré patterns for both bilayer and monolayer device. For bilayer device, the SDPs are at ns1 = 2.1510 12 cm −2 , and ns2 = 2.3410 12 cm −2 (Fig. 3a In the fan diagram ∂/∂B(n,B), we observed Brown-Zak (BZ) oscillations originating from the first order magnetic Bloch states (/0=1/q) for the two sets of SDPs in both bilayer and monolayer devices (Fig. S8b, e, f and Fig. S11b, e, f). The first order magnetic Bloch states formed by the two moiré superlattices are also prominent at high temperature (T = 70 K, where the Landau quantisation is suppressed), as indicated by the arrows in Figs. S9 and S12. We also observed high order magnetic Bloch states (/0 = 3/q) belonging to the two moiré superlattices, as shown in Fig. S8e, f and Fig. S11e, f.
In principle, the super-moiré pattern generated by the two original moiré patterns should have six possible reciprocal lattice vectors, as described in the previous study 5  3.1410 12 cm −2 , 5.810 12 cm −2 , 7.510 12 cm −2 and 8.910 12 cm −2 . We observed satellite peaks in ρxx near most of these carrier densities nsm (Fig. 3a in the main text and Fig. S10a). Note that near nsm corresponding to these resistivity peaks, most of the low-B ρxy regions have sign reversal. In addition, similar to bilayer device, in the plot of ∂xx/∂B(n,B) of monolayer device, we observe prominent horizontal streaks at B = 1.14 T and 0.57 T, which perfectly matches the magnetic field of the first and second order magnetic Bloch state originating from the super-moiré pattern when BA = 0 (where q = 1) and BA = 0/2 (where q = 2), respectively.
These features allow us to attribute the origin of the observed carrier densities nsm to the super-moiré pattern with different wavelengths. (1/q) 0 with q = 5 to 7. At low magnetic fields, the resolution in B is not sufficient to distinguish the two sets of oscillations.