Tunable crystal symmetry in graphene-boron nitride heterostructures with coexisting moir\'e superlattices

TL;DR

This study demonstrates how the electronic properties of graphene encapsulated in boron nitride can be precisely tuned by controlling the relative orientations of the layers, leading to complex moiré patterns and band structure modifications.

Contribution

It reveals the ability to manipulate band gaps and electronic structures in graphene-BN heterostructures through tunable moiré superlattices and layer alignment.

Findings

Band gaps depend on BN layer alignment (0° or 60°).

Small misalignments create multiple secondary Dirac points.

Interplay of moiré patterns enables electronic property control.

Abstract

In heterostructures consisting of atomically thin crystals layered on top of one another, lattice mismatch or rotation between the layers results in long-wavelength moir\'e superlattices. These moir\'e patterns can drive significant band structure reconstruction of the composite material, leading to a wide range of emergent phenomena including superconductivity, magnetism, fractional Chern insulating states, and moir\'e excitons. Here, we investigate monolayer graphene encapsulated between two crystals of boron nitride (BN), where the rotational alignment between all three components can be varied. We find that band gaps in the graphene arising from perfect rotational alignment with both BN layers can be modified substantially depending on whether the relative orientation of the two BN layers is 0 or 60 degrees, suggesting a tunable transition between the absence or presence of inversion symmetry in the heterostructure. Small deviations ($<1^{\circ}$) from perfect alignment of all three layers leads to coexisting long-wavelength moir\'e potentials, resulting in a highly reconstructed graphene band structure featuring multiple secondary Dirac points. Our results demonstrate that the interplay between multiple moir\'e patterns can be utilized to controllably modify the electronic properties of the composite heterostructure.

Figures

e), we observe both a sharp resistance peak at the primary Dirac point (PDP) and broad resistance peaks at finite density corresponding to the moiré-induced SDPs. We note that the resistivity of the hole-band SDP is larger than the electron-band SDP but smaller than the PDP, consistent with previous room temperature studies of graphene aligned to a single BN layer [Wang2015 ](#bib.bib16). As $\theta_{t}$ approaches zero (Figs. [2](#S0.F2)c-d), both the PDP and SDP grow, but their relationship inverts with $\rho_{SDP}$ exceeding $\rho_{PDP}$ very near $\theta_{t}=0^{\circ}$. By contrast, near $

f of the main text at higher densities which remain unexplained. While these may be related to replica features of the low density state around SDPs, their density spacing from SDP features is not identical to $N$.

Tables

Table 1. Table S1: Moiré wavelength for several top BN configurations in device R1.

${\theta }_t$	nSDP (1012 cm-2) (doubly-aligned)	nSDP,1 (1012 cm-2)	nSDP,2 (1012 cm-2)	$\lambda$ (nm) (doubly-aligned)	${\lambda }_1$ (nm)	${\lambda }_2$ (nm)
unencapsulated	n/a	2.05	n/a	n/a	15.0	n/a
-0.7∘	n/a	2.05	3.22	n/a	15.0	12.0
0∘	2.00	n/a	n/a	15.2	n/a	n/a
9.9∘	n/a	2.05	n/a	n/a	15.0	n/a
59.1∘	n/a	2.05	3.00	n/a	15.0	12.4
59.4∘	n/a	2.05	3.72	n/a	15.0	11.1
60∘	2.00	n/a	n/a	15.2	n/a	n/a

Equations

$$\Lambda =\frac{(1+D)a^2}{\sqrt{2(1+D)(1-cos\Theta )+D^2}}$$

$$D=\frac{{\lambda }_b-{\lambda }_t}{{\lambda }_t},{\lambda }_t\le {\lambda }_b$$

$$\tan {\varphi }_{t,b}=\frac{-\sin {\theta }_{t,b}}{(1+D)-\cos {\theta }_{t,b}}$$

$$\Theta =\left(|{\varphi }_t-{\varphi }_b|-30^{\circ }\right)\text{mod}(60^{\circ })-30^{\circ }$$

Full text

Tunable crystal symmetry in graphene–boron nitride heterostructures with coexisting moiré superlattices

Nathan R. Finney1∗

  

Matthew Yankowitz2∗

  

Lithurshanaa Muraleetharan1

  

K. Watanabe3

  

T. Taniguchi3

  

Cory R. Dean2†

  

James Hone1†

1Department of Mechanical Engineering, Columbia University, New York, NY, USA

2Department of Physics, Columbia University, New York, NY, USA

3National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan

*∗*These authors contributed equally to this work.

† cd2478@columbia.edu (C.R.D.); jh2228@columbia.edu (J.H.)

In heterostructures consisting of atomically thin crystals layered on top of one another, lattice mismatch or rotation between the layers results in long-wavelength moiré superlattices. These moiré patterns can drive significant band structure reconstruction of the composite material, leading to a wide range of emergent phenomena including superconductivity Cao2018b ; Yankowitz2019 ; Chen2019 , magnetism Sharpe2019 , fractional Chern insulating states Spanton2018 , and moiré excitons Seyler2019 ; Tran2019 ; Jin2019 ; Alexeev2019 . Here, we investigate monolayer graphene encapsulated between two crystals of boron nitride (BN), where the rotational alignment between all three components can be varied. We find that band gaps in the graphene arising from perfect rotational alignment with both BN layers can be modified substantially depending on whether the relative orientation of the two BN layers is 0 or 60 degrees, suggesting a tunable transition between the absence or presence of inversion symmetry in the heterostructure. Small deviations ($<1^{\circ}$) from perfect alignment of all three layers leads to coexisting long-wavelength moiré potentials, resulting in a highly reconstructed graphene band structure featuring multiple secondary Dirac points. Our results demonstrate that the interplay between multiple moiré patterns can be utilized to controllably modify the electronic properties of the composite heterostructure.

The ability to combine diverse vdW materials into a heterostructure enables engineering of new properties not observed in the constituent materials alone. A unique degree of freedom within these vdW heterostructures is the twist angle between layers, and changing this angle can strongly modify the material properties owing to the formation of moiré patterns. In graphene–BN heterostructures, the moiré pattern introduces a spatially-periodic effective potential that modifies the graphene band structure, giving rise to emergent secondary Dirac points (SDPs) at finite energy Yankowitz2012 ; Ponomarenko2013 ; Dean2013 ; Hunt2013 and band gaps at the charge neutrality point and valence band SDP Hunt2013 ; Woods2014 ; Chen2014 ; Wang2015 ; Wang2016 ; Song2013 ; Bokdam2014 ; Moon2014 ; Wallbank2015 ; Jung2015 ; SanJose2014 ; Slotman2015 ; Jung2017 ; Yankowitz2018 . In twisted bilayer graphene (tBLG), correlated insulating states and superconductivity emerge at a twist angle of $\sim$1.1*∘* where the lowest energy moiré bands become exceptionally flat Cao2018a ; Cao2018b ; Yankowitz2019 ; Sharpe2019 ; Lu2019 . However, typical vdW heterostructures comprising many flakes possess numerous crystal interfaces, and in principle multiple long-wavelength moiré patterns may coexist within a single heterostructure, likely with profound consequences on moiré-driven physics. For example, topological bands have been shown to potentially arise in tBLG aligned to BN Sharpe2019 . So far, little has been done to controllably tune the alignment of multiple pairs of crystals within a single device, and it is not well understood how multiple moiré patterns interact to influence the properties of the vdW heterostructure.

In a heterostructure where graphene is encapsulated on both sides by BN, there are a number of qualitatively distinct stacking orders that can be realized by independently controlling the twist angle of the graphene relative to the bottom BN, $\theta_{b}$, and the top BN, $\theta_{t}$ (angles illustrated schematically in Fig. 1a). In this work, we fabricate devices in which the graphene is aligned to the bottom BN ($\theta_{b}=0^{\circ}\pm 0.15^{\circ})$, but $\theta_{t}$ can be tuned to arbitrary angle. In this structure there are two unique positions in which the top BN layer may be “aligned” to the graphene — $\theta_{t}=0^{\circ}$ and $60^{\circ}$ — which have distinct symmetry due to the inequivalence of the boron and nitrogen triangular sublattices of the BN unit cell. Additionally, we study the case of small misalignment of the top layer (i.e. small but non-zero values of $\theta_{t}$) such that the two interfaces give rise to separate and incommensurate moiré patterns but with similar wavelengths.

As illustrated in Fig. 1b, we utilize a dynamically rotatable heterostructure to vary $\theta_{t}$ while maintaining fixed $\theta_{b}$. Graphene is first aligned to the bottom BN layer and shaped into a Hall bar geometry. The rotational alignment of the graphene and BN is initially determined by Raman spectroscopy prior to making electrical contact, and subsequently confirmed precisely by measuring the charge carrier density of the SDPs in electrical transport measurements (see Methods). The top BN is separately patterned into a circular shape with rectangular “handles,” and is subsequently transferred onto the graphene Hall bar such that it covers the entire active area of the channel (Fig. 1b). The top BN can be mechanically rotated to any angle ($0^{\circ}<\theta_{t}<360^{\circ}$) using an AFM tip, and electrical transport measurements can be performed simultaneously Chari2016 ; Ribeiro2018 . We are able to measure changes in the twist angle of the top BN ($\Delta\theta_{t}$) to better than 0.1*∘* from AFM topographs.

Fig. 1c shows spectra of the graphene 2D Raman peak at different values of $\theta_{t}$. In Fig. 1d we plot the full width at half maximum of the 2D peak (FWHM2D) versus $\theta_{t}$ over a 60*∘* range. Over most of this range the FWHM2D exhibits a constant value that is $\sim$20 cm*-1* larger than isolated graphene. This behavior is consistent with the presence of a single long-wavelength moiré pattern resulting from the fixed zero-angle alignment between the graphene and bottom BN Eckmann2013 ; Ribeiro2018 . When $\theta_{t}$ approaches $0^{\circ}$ or $60^{\circ}$, the linewidth is broadened by an additional $\sim$20 cm*-1*, with FWHM2D near 55 cm*-1*. We interpret this additional broadening to result from rotational alignment of the top BN layer to the graphene/bottom BN, and therefore provides an absolute measure of the top BN layer orientation. In previous studies, the broadening of the 2D mode was understood to arise from moiré-scale relaxations of the graphene lattice Eckmann2013 . Our observation of an approximate doubling of this broadening indicates that the graphene couples similarly to the moiré patterns from both the top and bottom encapsulating BN layers.

While the 2D Raman peak is 60*∘* periodic, the BN–graphene–BN trilayer lattice structure is not precisely equivalent upon 60*∘* rotations of the top BN crystal. As illustrated schematically in Figs. 2a-b, for $\theta_{t}=0^{\circ}$ the moiré pattern is three-fold rotationally symmetric, whereas for $\theta_{t}=60^{\circ}$ the moiré pattern is six-fold rotationally symmetric. The lattice structures at the high-symmetry points of the moiré patterns (outer schematics in Figs. 2a-b) highlight the important difference between these two cases. For $\theta_{t}=0^{\circ}$ the top layer B (N) atoms sit atop the bottom layer B (N) atoms and the overall structure breaks inversion symmetry. In contrast, for $\theta_{t}=60^{\circ}$ the top layer B (N) atoms sit atop N (B) atoms and the structure hosts inversion symmetry. The nature of inversion symmetry in graphene–BN heterostructures has previously been tied to band structure reconstruction of the graphene Hunt2013 ; Wang2015 , suggesting that the electronic properties of our devices may vary substantially with $\theta_{t}$.

We first investigate room temperature electrical transport of our device, and in particular observe notable differences between $\theta_{t}=0^{\circ}$ and 60*∘*. Figs. 2c-f show the device resistivity, $\rho$, as a function of charge carrier density, $n$, for various values of $\theta_{t}$. When the top BN is far from alignment (Fig. 2e), we observe both a sharp resistance peak at the primary Dirac point (PDP) and broad resistance peaks at finite density corresponding to the moiré-induced SDPs. We note that the resistivity of the hole-band SDP is larger than the electron-band SDP but smaller than the PDP, consistent with previous room temperature studies of graphene aligned to a single BN layer Wang2015 . As $\theta_{t}$ approaches zero (Figs. 2c-d), both the PDP and SDP grow, but their relationship inverts with $\rho_{SDP}$ exceeding $\rho_{PDP}$ very near $\theta_{t}=0^{\circ}$. By contrast, near $\theta_{t}=60^{\circ}$ (Fig. 2f), the SDP again grows, but less so, and the PDP appears slightly suppressed. The angle dependence of the PDP and hole-band SDP peaks are shown in more detail in Figs. 2g-h, where the peak resistance values are plotted over a full $360^{\circ}$ rotation. In particular, we note that both the PDP and hole-band SDP are maximal at the angle we label $\theta_{t}=0^{\circ}$. While translation of the top BN with respect to the bottom would change the overall stacking configuration, we observe nearly equivalent transport every $\Delta\theta_{t}=120^{\circ}$. This suggests translation does not play a significant role, and that the stacking configuration with no translational offset is the structural ground state at the aligned positions.

In graphene aligned to a single BN layer, the staggered sublattice potential of the BN breaks inversion symmetry in the graphene layer for both 0*∘* and 60*∘* “aligned” orientations, resulting in a band gap at the PDP whose value is expected to scale with the magnitude of effective superlattice potential Hunt2013 ; Woods2014 ; Chen2014 ; Wang2015 ; Wang2016 ; Song2013 ; Bokdam2014 ; Moon2014 ; Wallbank2015 ; Jung2015 ; SanJose2014 ; Slotman2015 ; Jung2017 ; Yankowitz2018 . We conjecture that the asymmetry between $\theta_{t}=0^{\circ}$ and 60∘ in our device correlates with the transition between the broken inversion symmetry structure at $\theta_{t}=\theta_{b}=0^{\circ}$ (Fig. 2a) — in which the PDP gap is likely to be largest — and the inversion symmetric structure at $\theta_{t}=60^{\circ}$ (Fig. 2b) — in which no PDP gap would be expected within a single-particle model.

Figs. 3a-d compare the temperature dependence of the PDP and SDP at $\theta_{t}=0^{\circ}$ and 60*∘*. Strongly insulating behavior is observed in both at $\theta_{t}=0^{\circ}$, whereas only weakly insulating behavior is observed for $\theta_{t}=60^{\circ}$. A linear fit to the thermally activated regime (black dashed lines in Figs. 3e-f) gives a measure of the activation gap, $\Delta$. In Fig. 3g we plot the gaps for 5 different values of $\theta_{t}$ (square markers). Both the PDP and SDP gaps are notably enhanced at $\theta_{t}=0^{\circ}$, whereas at $\theta_{t}=60^{\circ}$ both are notably reduced. The gaps have little dependence on $\theta_{t}$ at all other angles. Our observation of $\Delta_{PDP}>60$ meV is so far the largest gap observed in a pristine graphene device, and may potentially be significantly further enhanced under pressure Yankowitz2018 . Although the gaps extracted from an Arrhenius fit are not identically zero at $\theta_{t}=60^{\circ}$, the device only exhibits activated behavior over well less than a decade change in conductance, hence we expect our reported gaps to be an upper bound at this angle.

The band structure modification anticipated from symmetry considerations is illustrated schematically in Fig. 3h, where we anticipate the largest gaps for $\theta_{t}=\theta_{b}=0^{\circ}$ owing to the doubled moiré potential, while at $\theta_{t}=60^{\circ}$ inversion symmetry protects the Dirac crossings. Following this simple expectation, the measured gaps corroborate our previous assignments of $\theta_{t}=0^{\circ}$ and $60^{\circ}$. We additionally measure the band gaps in two “stationary” devices (i.e. without the ability to dynamically rotate the top BN, see Supplementary Section 2), which also exhibit large broadening of the FWHM2D Raman peak. The gaps are similarly either enhanced or suppressed (diamond and triangle markers in Fig. 3c), suggesting this effect is generic for samples in which both BN layers are aligned to graphene.

Finally, we investigate the small angle limit ($0^{\circ}<\theta_{t}<1^{\circ}$) where the top and bottom BN yield moiré patterns with only slightly different period. Fig. 4b shows the hole-doped transport for three values of $\theta_{t}\neq 0^{\circ}$ in this regime. In addition to a peak in the resistivity at the PDP, we observe two sizable peaks at finite hole-doped densities denoted $n_{SDP,1}$ and $n_{SDP,2}$. $n_{SDP,1}$ is independent of $\theta_{t}$ and corresponds to the SDP arising from the moiré potential of the perfectly aligned graphene and bottom BN ($\theta_{b}\approx 0^{\circ}$), while the resistance peaks at $n_{SDP,2}$ correspond to the moiré potential from the top BN. The position of $n_{SDP,2}$ as a function of $\theta_{t}$ is shown the inset of Fig. 4b, and is in good quantitative agreement with the anticipated dispersion (see Supplementary Section 4) Yankowitz2012 . We note that in the small angle limit, thermally activated behavior is observed at both $n_{SDP,1}$ and $n_{SDP,2}$, suggesting band gaps associated with each SDP with typical values of $\sim$20 meV and $\sim$5 meV, respectively (see Supplementary Fig. 5). This implies the reconstructed band structure shown schematically in Fig. 4c in which an isolated narrow band emerges between the two superlattice gaps $\Delta_{SDP,1}$ and $\Delta_{SDP,2}$. Further theoretical and experimental effort is necessary to fully explore the consequence of this band reconstruction, however the appearance of a flat band whose width varies with rotation angle provides the intriguing possibility of hosting tunable correlated states at low temperature Cao2018a ; Chen2019 .

In a magnetic field, we observe sequences of quantum oscillations emerging from the Dirac points at $n=0$, $n_{SDP,1}$, and $n_{SDP,2}$ (Fig. 4d-e, for $\theta_{t}=59.4^{\circ}$). We further observe weak signatures of resistive states adjacent to each of the three Dirac points, marked by orange arrows. To identify these resistive states more clearly, we plot $d\rho/dT$ versus $T$ and $n$ in Fig. 4f. In addition to the insulating states previously discussed at the PDP and the two SDPs, we observe a negative temperature dependence at the same densities marked by orange arrows in Fig. 4d, suggesting the presence of new insulating-like states. The two coexisting moiré patterns may in principle interfere to produce a second-order moiré pattern with a very long period (Fig. 4a), inducing an additional resistive state at low density Wang2019 . However, the position of the low density resistive feature does not agree precisely with quantitative theoretical estimates for all measured $\theta_{t}$ (see Supplementary Section 6). The disagreement may arise due to unexpected structural reconstructions of the second order moiré pattern, however at this point we are unable to understand these new resistive states in detail.

In conclusion, we demonstrate the ability to induce and control multiple moiré patterns within a BN–graphene–BN heterostructure. In particular, we are able to dynamically tune the crystal symmetry of the composite material by realizing distinct stacking configurations of the three layers, and further induce coexisting moiré patterns which combine to strongly modify the graphene band structure. Our techniques for engineering multiple moiré patterns within a single vdW heterostructure are quite general and can easily be extended to a wide variety of other device structures, motivating a new class of experiments in which the twist angle of multiple crystal interfaces can be tuned to realize novel material properties.

Methods

All heterostructures are assembled using standard dry-transfer techniques with a poly-propylene carbonate (PPC) film on a polydimethyl siloxane (PDMS) stamp Wang2013 , and rest atop a Si/SiO2 substrate. The device fabrication of the rotatable devices largely follows the techniques developed in Refs. Chari2016 ; Ribeiro2018 . Device R1 consists of a graphene Hall bar on a $\sim$44 nm thick BN resting on a $\sim$11 nm thick graphite gate. A $\sim$56 nm thick BN flake capped by $\sim$40 nm of hydrogen silsesquioxane (HSQ) is subsequently transferred onto the graphene Hall bar. The HSQ cap acts as the etch mask to shape the rotating BN, and also provides additional durability during AFM pushes. “Stationary” devices consist of a graphene Hall bar fully encapsulated by BN, all atop a graphite gate. Electrical contact to all devices is made by Cr/Au (2 nm/100 nm). Supplementary Sections 1-2 and Supplementary Figs. 1-2 provide full details of the device fabrication.

Raman spectroscopy measurements are performed at room temperature in air. All measurements are acquired using a 532 nm wavelength laser with a power $<$ 1 mW. Measurements to extract the graphene/BN twist angle in stationary devices are performed before the addition of the graphite back gate. In the rotatable device, the encapsulated graphene region sits atop the graphite gate. To isolate the Raman response from the graphene, we separately acquire a Raman spectrum from a nearby region of the graphite gate without the graphene and subtract this background response. We extract FWHM2D from a Lorentzian fit (see Supplementary Section 3 and Supplementary Fig. 3 for full details).

Transport measurements are conducted in a four-terminal geometry with ac current excitation of 10-100 nA using standard lock-in technique at 17.7 Hz. The graphene contact regions (which extend beyond the graphite bottom gate) are gated to a high carrier density by applying a bias to the silicon substrate (typically 5-50 V for SiO2 thickness of $\sim$285 nm) to reduce the contact resistance. We extract $n(V_{G})$ by fitting the dispersion of the quantum Hall states in high magnetic field as $n=\nu eB/h$, where $\nu$ is the filling factor, $h$ is Planck’s constant, and $e$ is the elementary charge. The moiré wavelength, $\lambda$, is calculated using the geometric relation $\lambda^{2}=8/(n_{SDP}\sqrt{3})$ where nSDP is the density at full filling of the moiré unit cell.

acknowledgments

The authors thank Rebeca Ribeiro-Palau, Changjian Zhang, and Shaowen Chen for technical support, as well as Jeil Jung, Mikito Koshino and Chris Marianetti for helpful discussions. This work was primarily supported by the NSF MRSEC program through Columbia in the Center for Precision Assembly of Superstratic and Superatomic Solids (DMR-1420634). Sample device design and fabrication was partially supported by DoE Pro-QM EFRC (DE-SC0019443). N.F. acknowledges support from the Stewardship Science Graduate Fellowship program provided under cooperative agreement number DE-NA0002135. CRD acknowledges the support of the David and Lucile Packard Foundation. K.W. and T.T. acknowledge support from the Elemental Strategy Initiative conducted by the MEXT, Japan and the CREST (JPMJCR15F3), JST.

Author contributions

N.R.F. and L.M. fabricated the devices. N.R.F. and M.Y. performed the measurements and analyzed the data. K.W. and T.T. grew the hBN crystals. C.R.D. and J.H. advised on the experiments. The manuscript was written with input from all authors.

Competing interests

The authors declare no competing interests.

Supplementary Information

S1 Fabrication of rotatable devices

To fabricate rotatable devices R1 and R2, we begin by sequentially picking up flakes of BN, graphite, BN, and graphene using a PPC/PDMS stamp. The PPC film is then mechanically removed from the PDMS stamp and placed onto a Si/SiO2 substrate, such that the graphene layer is exposed at the top of the heterostructure. The underlying PPC is removed by vacuum annealing at $T=350^{\circ}$ C. In all processed devices the silicon is used as a global gate for the graphene contacts, and the graphite is used as a local gate for the graphene channel (Fig. S1a-b). The use of a graphite gates has previously been demonstrated to significantly improve the charge homogeneity of graphene devices Zibrov2017 .

Devices with good rotational alignment of the graphene and BN are next identified using Raman spectroscopy Eckmann2013 ; Ribeiro2018 . Devices in which the graphene and BN are well aligned are then shaped into a Hall bar geometry using an oxygen plasma etch through a PMMA mask defined by standard electron-beam lithography (Fig. S1c). In order to realize more robust electrical contacts, the BN surrounding the graphene contacts is further etched using a CHF3/O2 plasma through a PMMA mask (Fig. S1d). All electrical contacts in device R1 show no significant change over 8 thermal cycles ($T$ = 300 K - 1.5 K), while many of the contacts in device R2 (in which no BN etch is performed) became open after thermal cycling.

We exfoliate BN onto a separate Si/SiO2 substrate and etch an array of BN “rotators” using a protective HSQ mask (Figs. S1f-g). The rotators are picked up using a separate PPC/PDMS stamp (Fig. S1h) and transferred onto the graphene Hall bar (Fig. S1i). A single rotator is aligned as desired onto the graphene Hall bar by positioning with an AFM tip in contact mode. The device is annealed a second time before deposition of Cr/Au (2 nm/100 nm) metal electrodes (Fig. S1j). The graphene channel is further cleaned by sequentially pushing a BN rotator across the surface to collect and remove interfacial contamination.

S2 Fabrication of stationary devices

To fabricate stationary devices S1-S4, we begin by identifying flakes of BN which have fractured into two pieces during the mechanical exfoliation procedure (Fig. S2a). These fractured BN pieces share common edge profiles which can be easily optically aligned. One piece of the fractured BN is first picked up by a PPC/PDMS stamp (Fig. S2b-c), leaving the remaining piece on the Si/SiO2 wafer. A flake of graphene is then added to the heterostructure, and is chosen and aligned to fit completely within the BN area (Fig. S2d). Next, the two BN flakes are rotationally aligned optically (with precision better than 0.5*∘*) and the remaining BN is added to the heterostructure (Fig. S2e). The heterostructure is transferred to a Si/SiO2 substrate and annealed at $T=350^{\circ}$ C. No effort is made to align crystal edges of graphene to the BN, however, because the graphene rests entirely within both BN flakes it tends to rotate to alignment with the BN during the transfer or annealing Wang2015 . Devices with good rotational alignment of the graphene and BN are identified with Raman spectroscopy. A new PPC/PDMS stamp is then used to pick up the heterostructure, followed by a graphite flake (Fig. S2f), and the heterostructure is transferred to a separate Si/SiO2 substrate. The device is then shaped into a Hall bar geometry and electrically contacted following the procedures described for devices R1 and R2 (Figs. S2g-h).

S3 Analysis of Raman spectroscopy

Fig. S3a plots Raman spectra near the 2D graphene peak in device R1. The measured signal (purple curve) features contributions from both the graphene channel and the graphite back gate. To isolate the graphene response (black curve), we separately acquire a Raman spectrum on the graphite gate alone (green curve) and subtract that contribution. The insets show the position of the laser (green dot) during acquisitions of the purple and green curves. We fit the graphene 2D peak with a Lorentzian (red dashed curve) to extract FWHM2D. The reported values of FWHM2D are averaged across a few (3-6) measured background spots. Fig. S3b shows the 2D Raman peak for stationary BN encapsulated graphene devices in three distinct stacking configurations. These spectra are acquired prior to the addition of graphite back gates.

Fig. S3c plots the graphene FWHM2D at several alignment configurations for two different rotatable devices. The open triangle markers are taken from Ref. Ribeiro2018 , in which the graphene is misaligned from the bottom BN. The closed circle markers reproduce the data from device R1 shown in Fig. 1d of the main text. The dashed lines show FWHM2D of stationary devices S1 (graphene aligned to both BN layers), S3 (graphene aligned to a single BN layer), and S4 (graphene not aligned to either BN layer).

S4 Determination of moiré wavelength and zero top BN angle assignment

We extract the moiré wavelength in our devices via the geometric relation $\lambda^{2}=8/(n_{SDP}\sqrt{3})$. We accurately determine $n_{SDP}$ by projecting the quantum oscillations from the SDP Landau fan to $B=0$. Table S1 lists $n_{SDP}$ and $\lambda$ for several values of $\theta_{t}$ in device R1. We additionally measure $n_{SDP}$ in a region of the graphene Hall bar which is not encapsulated by the top BN rotator, such that we are sensitive only to $\theta_{b}$, and extract a corresponding $\lambda$ = 15.0 nm. This implies nearly perfect alignment of the graphene and bottom BN. The maximum uncertainty in moiré wavelength calculated at each configuration is estimated to be $\pm 0.1$ nm.

We observe identical $n_{SDP}$ for all $\theta_{t}$ except at the angles labeled 0*∘* and 60*∘, in which $n_{SDP}$ becomes marginally smaller. We take the smallest measured value of $n_{SDP}$ to correspond to perfect alignment of the graphene and top BN. Using $\theta_{t}=0^{\circ}$ and the calculated $\lambda$ = 15.2 nm, we use the geometric relation derived in Ref. Yankowitz2012 to extract the lattice mismatch between the graphene and BN to be $\delta\approx$ 1.65%. We then estimate the misalignment in $\theta_{b}$ to be less than $\pm$0.15∘*, although we are not sensitive to the sign of the misalignment. We use these values to generate the black dashed curve in the inset of Fig. 4b in the main text.

We note that if we instead assume $\theta_{b}\equiv 0^{\circ}$, then the corresponding $n_{SDP,1}=2.05\times 10^{12}$ cm*-2* sets a lattice mismatch of $\delta\approx$ 1.67%. The larger 15.2 nm wavelength computed for the $\theta_{t}=0^{\circ}$ and 60*∘* configurations (where $n_{SDP}=2.00\times 10^{12}$ cm*-2*) could instead be understood to arise owing to a decrease in the lattice mismatch ($\sim$1.3%), indicating the onset of a commensurate transition, which has not yet been observed in exfoliated graphene–BN heterostructures. However, we are not able to reliably distinguish between these two possibilities within our experimental resolution.

The values of $\theta_{b}$ and $\theta_{t}$ are more challenging to precisely determine in the stationary devices owing to the lack of dynamic rotational control. Our method of using fractured BN flakes for both the top and bottom encapsulating layers ensures that the relative twist angle of the two BN layers can be determined optically to better than 1*∘. We use the broadened Raman 2D peak to first identify that the graphene is well aligned to both BN layers. However, because stacking sequence of bulk BN is AA’ (in which boron and nitrogen atoms alternate between stacked layers) Constantinescu2013 , there is ambiguity in whether $\theta_{t}=0^{\circ}$ or 60∘* in our stationary devices depending on whether the number of layers in the BN flake is even or odd. We therefore rely on the measured band gap — in particular whether it is enhanced or suppressed from the typical value of a graphene device aligned with a single BN — to distinguish between $\theta_{t}=0^{\circ}$ and 60*∘*.

Figs. S4a-b plot Landau fan diagrams for stationary devices S1 and S2. Graphene is aligned to both BN layers in these devices, with $\theta_{t}=0^{\circ}$ in device S1 and 60*∘* in device S2. Despite the difference in stacking configuration, we observe qualitatively similar features in both fan diagrams, with sequences of quantum oscillations emerging from the PDP and from the SDPs. Although there may be quantitative differences in the Landau level gaps at the SDPs compared to devices in which graphene is aligned to only a single BN layer due to the doubled moiré potential strength, we have not investigated these effects in detail in this study.

S5 Band gaps for small top BN twist angles

Fig. S5 shows the band gaps measured by thermal activation at the PDP and both valence band SDPs in the case of small, non-zero $\theta_{t}$. In the case of perfect alignment, only a single SDP gap is observed, as discussed in the main text. We find that the $\Delta_{SDP,2}$ are much smaller than the corresponding $\Delta_{SDP,1}$, however we have not measured these gaps at enough values of $\theta_{t}$ to conclusively determine a trend with twist angle.

S6 Second order moiré patterns

Coexisting moiré patterns with different wavelengths may in principle interfere to produce a second-order moiré pattern. We determine the wavelength, $\Lambda$, using the dispersion in Ref. Yankowitz2012 with lattice constants taken as the moiré wavelength of the graphene and bottom BN, $\lambda_{b}$, and the graphene and top BN, $\lambda_{t}$. The twist angle between the constituent moiré patterns, $\Theta$, is a function of the twist angle of the top and bottom moiré patterns relative to the graphene ($\phi_{t}$ and $\phi_{b}$, respectively), determined by $\theta_{t}$, $\theta_{b}$, and $\delta$ Yankowitz2012 . Values of $|\phi_{t}-\phi_{b}|$, must correctly map to values of $\Theta$ on a 60*∘* interval centered at 0*∘* for the dispersion to predict $\Lambda$ Wang2019 . Accounting for this, we have

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

The modular division used here to find $\Theta$ reproduces the output of the piecewise conditional reported in Wang2019 , with a negligible error produced at the endpoints of the interval.

Fig. S6a shows second order moiré patterns for two combinations of top and bottom BN twist angles, with values of $\Lambda$ calculated using the scheme described above. We note that there is an inequivalence in second order moiré wavelength for $\theta_{b}=+$0.15*∘* and $\theta_{b}=-$0.15*∘, owing to the difference in moiré twist angle $\phi_{b}$ for each case. Therefore for the top BN twist angle configuration shown ($\theta_{t}=$ 0.6∘*), there are two possible second order moiré wavelengths possible if the sign of the bottom BN twist angle is unknown.

Fig. S6b plots the density at full filling of the second order moiré unit cell, $N$, determined by $N$ = 8/($\Lambda^{2}\sqrt{3}$) as a function of $\theta_{t}$ and $\theta_{b}$, with an assumed graphene/BN lattice mismatch of $\delta\approx$ 1.65%. Our transport measurements exhibit dominant resistive peaks corresponding to densities $n_{SDP,1}$ and $n_{SDP,2}$ (Figs. 4b-f and  S6c-h) with at least weakly activated behavior and corresponding sequences of quantum oscillations. Furthermore, the measured top BN twist angles corresponding to $n_{SDP,2}$ are consistent with top BN twist angles predicted by the dispersion in Ref. Yankowitz2012 , and we observe that $n_{SDP,1}$ is fixed by $\theta_{b}$ for all values of $\theta_{t}$. We therefore conduct our analysis of possible second order moiré effects on the basis that the densities labeled $n_{SDP,1}$ and $n_{SDP,2}$ correspond to the density at full filling of the bottom and top moiré unit cells, respectively.

To compute the expected values of $N$ for both $\theta_{b}>0^{\circ}$ and $\theta_{b}<0^{\circ}$, we first take the uncertainty in $\delta$, $n_{SDP,1}$, and $n_{SDP,2}$ into account. We estimate the uncertainty in $n_{SDP,1}$ and $n_{SDP,2}$ to be $\sim\pm 0.02\times 10^{12}$cm*-2*, corresponding to the uncertainty in fitting the trajectories of the quantum oscillations in the Landau fan diagram. Using a standard differential error propagation we estimate the range of values of $\theta_{b}$ and $\theta_{t}$ over a range of possible $\delta$ (ranging from $\sim$1.637% to $\sim$1.654%). This sets a range of possible values of $N$ at each value of $\theta_{t}$.

Figs. S6c-e show $d\rho/dT$ as a function of $T$ and $n$ for $\theta_{t}=$ 59.1*∘, -0.7∘* and 59.4*∘*. We observe regions of negative $d\rho/dT$ (orange arrows, indicating insulating-like behavior) in the vicinity of the expected values of $N$. Figs. S6f-h show $\rho$ as a function of $n$ at $T$ = 1.7 K for the same values of $\theta_{t}$. The bounds of the negative $d\rho/dT$ regions observed in Figs. S6c-e are marked by orange dashed lines. The computed range of possible values of $N$ at each value of $\theta_{t}$ are shown for $\theta_{b}<0^{\circ}$ (red shaded region) and $\theta_{b}>0^{\circ}$ (blue shaded region).

We find that there exists some overlap between the computed values of $N$ and the observed regions of negative $d\rho/dT$, suggesting a possible correspondence with second-order moiré features. However, for fixed $\delta$ and $n_{SDP,1}$, there appears to be no combination of choices for $n_{SDP,2}$ that generate values of $N$ matching negative $d\rho/dT$ features for all three configurations ($\theta_{t}=59.1^{\circ}$, -0.7*∘, and 59.4∘*) simultaneously. Furthermore, there are additional regions of negative $d\rho/dT$ marked by orange arrows in Fig. 4f of the main text at higher densities which remain unexplained. While these may be related to replica features of the low density state around SDPs, their density spacing from SDP features is not identical to $N$.

We additionally plot the transport over a full range of accessible gate voltage in Fig. S7. We observe three features for electron-type doping roughly corresponding to the two SDPs and a lower density feature (black arrows). However, unexpectedly the exact densities of these features do not exactly correspond with their hole-doped counterparts. Therefore, a more detailed study with small $\theta_{t}$ is necessary to understand all of these extra resistive states in more detail, and in particular how they relate to second order moiré patterns.

S7 Room temperature transport in device R2

We measure the room temperature transport of a second rotatable device, R2. Fig. S8a-b shows $\rho_{PDP}$ and $\rho_{SDP}$ as a function of $\theta_{t}$. Similar to device R1, we observe extrema in both $\rho_{PDP}$ and $\rho_{SDP}$ with 60*∘* periodicity. Following the same procedure, we define $\theta_{t}=0^{\circ}$ as the angle corresponding to the maximum values observed in both $\rho_{PDP}$ and $\rho_{SDP}$ over a given 120*∘* interval. We observe minima in $\rho_{PDP}$ as well as suppressed enhancement of $\rho_{SDP}$ near $\theta_{t}=60^{\circ}$ in both devices.

Fig. S8c plots gate sweeps at various $\theta_{t}$ in device R2. Fig. S8d shows a zoomed-in view of the PDP peaks from Fig. S8c. We observe a slight suppression in $\rho_{PDP}$ very near $\theta_{t}=180^{\circ}$ compared with its value at large misalignment (151.4*∘*). Comparable behavior is observed in device R1, however at present we do not have a complete model to understand the room-temperature transport response, which is complicated by scattering from acoustic phonons in the graphene and from polar optical phonons in the BN substrate.