Tunable crystal symmetry in graphene-boron nitride heterostructures with coexisting moiré superlattices
Nathan R. Finney Matthew Yankowitz Lithurshanaa Muraleetharan K. Watanabe T. Taniguchi Cory R. Dean James Hone





 

 
 
 

  

^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.)





Tunable crystal symmetry in graphene–boron nitride heterostructures with coexisting moiré superlattices
    James Hone^1†
    December 30, 2023
=======================================================================================================




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 <cit.>, magnetism <cit.>, fractional Chern insulating states <cit.>, and moiré excitons <cit.>. 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^∘) 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 <cit.> and band gaps at the charge neutrality point and valence band SDP <cit.>. In twisted bilayer graphene (tBLG), correlated insulating states and superconductivity emerge at a twist angle of ∼1.1^∘ where the lowest energy moiré bands become exceptionally flat <cit.>. 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 <cit.>. 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, θ_b, and the top BN, θ_t (angles illustrated schematically in Fig. <ref>a). In this work, we fabricate devices in which the graphene is aligned to the bottom BN (θ_b=0^∘±0.15^∘), but θ_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 — θ_t=0^∘ and 60^∘ — 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 θ_t) such that the two interfaces give rise to separate and incommensurate moiré patterns but with similar wavelengths.





As illustrated in Fig. <ref>b, we utilize a dynamically rotatable heterostructure to vary θ_t while maintaining fixed θ_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. <ref>b). The top BN can be mechanically rotated to any angle (0^∘<θ_t<360^∘) using an AFM tip, and electrical transport measurements can be performed simultaneously <cit.>. We are able to measure changes in the twist angle of the top BN (Δθ_t) to better than 0.1^∘ from AFM topographs. 

Fig. <ref>c shows spectra of the graphene 2D Raman peak at different values of θ_t. In Fig. <ref>d we plot the full width at half maximum of the 2D peak (FWHM_2D) versus θ_t over a 60^∘ range. Over most of this range the FWHM_2D exhibits a constant value that is ∼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 <cit.>. When θ_t approaches 0^∘ or 60^∘, the linewidth is broadened by an additional ∼20 cm^-1, with FWHM_2D 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 <cit.>. 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. <ref>a-b, for θ_t=0^∘ the moiré pattern is three-fold rotationally symmetric, whereas for θ_t=60^∘ the moiré pattern is six-fold rotationally symmetric. The lattice structures at the high-symmetry points of the moiré patterns (outer schematics in Figs. <ref>a-b) highlight the important difference between these two cases. For θ_t=0^∘ the top layer B (N) atoms sit atop the bottom layer B (N) atoms and the overall structure breaks inversion symmetry. In contrast, for θ_t=60^∘ 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 <cit.>, suggesting that the electronic properties of our devices may vary substantially with θ_t.

We first investigate room temperature electrical transport of our device, and in particular observe notable differences between θ_t=0^∘ and 60^∘. Figs. <ref>c-f show the device resistivity, ρ, as a function of charge carrier density, n, for various values of θ_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 <cit.>. As θ_t approaches zero (Figs. <ref>c-d), both the PDP and SDP grow, but their relationship inverts with ρ_SDP exceeding ρ_PDP very near θ_t=0^∘. By contrast, near θ_t=60^∘ (Fig. <ref>f), 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. <ref>g-h, where the peak resistance values are plotted over a full 360^∘ rotation. In particular, we note that both the PDP and hole-band SDP are maximal at the angle we label θ_t=0^∘. While translation of the top BN with respect to the bottom would change the overall stacking configuration, we observe nearly equivalent transport every Δθ_t=120^∘. 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 <cit.>. We conjecture that the asymmetry between θ_t=0^∘ and 60^∘ in our device correlates with the transition between the broken inversion symmetry structure at θ_t=θ_b=0^∘ (Fig. <ref>a) — in which the PDP gap is likely to be largest — and the inversion symmetric structure at θ_t=60^∘ (Fig. <ref>b) — in which no PDP gap would be expected within a single-particle model. 





Figs. <ref>a-d compare the temperature dependence of the PDP and SDP at θ_t=0^∘ and 60^∘. Strongly insulating behavior is observed in both at θ_t=0^∘, whereas only weakly insulating behavior is observed for θ_t=60^∘. A linear fit to the thermally activated regime (black dashed lines in Figs. <ref>e-f) gives a measure of the activation gap, Δ. In Fig. <ref>g we plot the gaps for 5 different values of θ_t (square markers). Both the PDP and SDP gaps are notably enhanced at θ_t=0^∘, whereas at θ_t=60^∘ both are notably reduced. The gaps have little dependence on θ_t at all other angles. Our observation of Δ_PDP>60 meV is so far the largest gap observed in a pristine graphene device, and may potentially be significantly further enhanced under pressure <cit.>. Although the gaps extracted from an Arrhenius fit are not identically zero at θ_t=60^∘, 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. <ref>h, where we anticipate the largest gaps for θ_t=θ_b=0^∘ owing to the doubled moiré potential, while at θ_t=60^∘ inversion symmetry protects the Dirac crossings. Following this simple expectation, the measured gaps corroborate our previous assignments of θ_t=0^∘ and 60^∘. 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 FWHM_2D Raman peak. The gaps are similarly either enhanced or suppressed (diamond and triangle markers in Fig. <ref>c), suggesting this effect is generic for samples in which both BN layers are aligned to graphene.





Finally, we investigate the small angle limit (0^∘<θ_t<1^∘) where the top and bottom BN yield moiré patterns with only slightly different period. Fig. <ref>b shows the hole-doped transport for three values of θ_t≠0^∘ 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 θ_t and corresponds to the SDP arising from the moiré potential of the perfectly aligned graphene and bottom BN (θ_b≈0^∘), 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 θ_t is shown the inset of Fig. <ref>b, and is in good quantitative agreement with the anticipated dispersion (see Supplementary Section 4) <cit.>. 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 ∼20 meV and ∼5 meV, respectively (see Supplementary Fig. 5). This implies the reconstructed band structure shown schematically in Fig. <ref>c in which an isolated narrow band emerges between the two superlattice gaps Δ_SDP,1 and Δ_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 <cit.>. 

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. <ref>d-e, for θ_t=59.4^∘). 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ρ/dT versus T and n in Fig. <ref>f. 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. <ref>d, 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. <ref>a), inducing an additional resistive state at low density <cit.>. However, the position of the low density resistive feature does not agree precisely with quantitative theoretical estimates for all measured θ_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 <cit.>, and rest atop a Si/SiO_2 substrate. The device fabrication of the rotatable devices largely follows the techniques developed in Refs. <cit.>. Device R1 consists of a graphene Hall bar on a ∼44 nm thick BN resting on a ∼11 nm thick graphite gate. A ∼56 nm thick BN flake capped by ∼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 FWHM_2D 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 SiO_2 thickness of ∼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=ν eB/h, where ν is the filling factor, h is Planck's constant, and e is the elementary charge. The moiré wavelength, λ, is calculated using the geometric relation λ^2 = 8/(n_SDP√(3)) where n_SDP 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.

naturemag












 




§ SUPPLEMENTARY INFORMATION












 §.§ 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/SiO_2 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^∘ 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. <ref>a-b). The use of a graphite gates has previously been demonstrated to significantly improve the charge homogeneity of graphene devices <cit.>.

Devices with good rotational alignment of the graphene and BN are next identified using Raman spectroscopy <cit.>. 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. <ref>c). In order to realize more robust electrical contacts, the BN surrounding the graphene contacts is further etched using a CHF_3/O_2 plasma through a PMMA mask (Fig. <ref>d). 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/SiO_2 substrate and etch an array of BN “rotators” using a protective HSQ mask (Figs. <ref>f-g). The rotators are picked up using a separate PPC/PDMS stamp (Fig. <ref>h) and transferred onto the graphene Hall bar (Fig. <ref>i). 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. <ref>j). The graphene channel is further cleaned by sequentially pushing a BN rotator across the surface to collect and remove interfacial contamination.



 §.§ 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. <ref>a). 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. <ref>b-c), leaving the remaining piece on the Si/SiO_2 wafer. A flake of graphene is then added to the heterostructure, and is chosen and aligned to fit completely within the BN area (Fig. <ref>d). 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. <ref>e). The heterostructure is transferred to a Si/SiO_2 substrate and annealed at T=350^∘ 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 <cit.>. 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. <ref>f), and the heterostructure is transferred to a separate Si/SiO_2 substrate. The device is then shaped into a Hall bar geometry and electrically contacted following the procedures described for devices R1 and R2 (Figs. <ref>g-h).







 §.§ Analysis of Raman spectroscopy


Fig. <ref>a 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 FWHM_2D. The reported values of FWHM_2D are averaged across a few (3-6) measured background spots. Fig. <ref>b 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. <ref>c plots the graphene FWHM_2D at several alignment configurations for two different rotatable devices. The open triangle markers are taken from Ref. <cit.>, in which the graphene is misaligned from the bottom BN. The closed circle markers reproduce the data from device R1 shown in Fig. <ref>d of the main text. The dashed lines show FWHM_2D 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). 













 §.§ Determination of moiré wavelength and zero top BN angle assignment


We extract the moiré wavelength in our devices via the geometric relation λ^2 = 8/(n_SDP√(3)). We accurately determine n_SDP by projecting the quantum oscillations from the SDP Landau fan to B=0. Table <ref> lists n_SDP and λ for several values of θ_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 θ_b, and extract a corresponding λ = 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 ±0.1 nm.

We observe identical n_SDP for all θ_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 θ_t=0^∘ and the calculated λ = 15.2 nm, we use the geometric relation derived in Ref. <cit.> to extract the lattice mismatch between the graphene and BN to be δ≈ 1.65%. We then estimate the misalignment in θ_b to be less than ±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. <ref>b in the main text.

We note that if we instead assume θ_b≡0^∘, then the corresponding n_SDP,1=2.05×10^12 cm^-2 sets a lattice mismatch of δ≈ 1.67%. The larger 15.2 nm wavelength computed for the θ_t = 0^∘ and 60^∘ configurations (where n_SDP=2.00×10^12 cm^-2) could instead be understood to arise owing to a decrease in the lattice mismatch (∼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 θ_b and θ_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) <cit.>, there is ambiguity in whether θ_t=0^∘ 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 θ_t=0^∘ and 60^∘.

Figs. <ref>a-b plot Landau fan diagrams for stationary devices S1 and S2. Graphene is aligned to both BN layers in these devices, with θ_t=0^∘ 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.







 §.§ Band gaps for small top BN twist angles


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











 §.§ 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, Λ, using the dispersion in Ref. <cit.> with lattice constants taken as the moiré wavelength of the graphene and bottom BN, λ_b, and the graphene and top BN, λ_t. The twist angle between the constituent moiré patterns, Θ, is a function of the twist angle of the top and bottom moiré patterns relative to the graphene (ϕ_t and ϕ_b, respectively), determined by θ_t, θ_b, and δ <cit.>. Values of |ϕ_t -ϕ_b|, must correctly map to values of Θ on a 60^∘ interval centered at 0^∘ for the dispersion to predict Λ <cit.>. Accounting for this, we have

    Λ = (1+D)a^2/√(2(1+D)(1-cosΘ) + D^2)

where

    D = λ_b-λ_t/λ_t,  λ_t ≤λ_b


    tanϕ_t,b = -sinθ_t,b/(1+D)-cosθ_t,b


    Θ = (|ϕ_t-ϕ_b|-30^∘) mod (60^∘) - 30^∘


The modular division used here to find Θ reproduces the output of the piecewise conditional reported in <cit.>, with a negligible error produced at the endpoints of the interval.

Fig. <ref>a shows second order moiré patterns for two combinations of top and bottom BN twist angles, with values of Λ calculated using the scheme described above. We note that there is an inequivalence in second order moiré wavelength for θ_b=+0.15^∘ and θ_b=-0.15^∘, owing to the difference in moiré twist angle ϕ_b for each case. Therefore for the top BN twist angle configuration shown (θ_t = 0.6^∘), there are two possible second order moiré wavelengths possible if the sign of the bottom BN twist angle is unknown.

Fig. <ref>b plots the density at full filling of the second order moiré unit cell, N, determined by N = 8/(Λ^2√(3)) as a function of θ_t and θ_b, with an assumed graphene/BN lattice mismatch of δ≈ 1.65%. Our transport measurements exhibit dominant resistive peaks corresponding to densities n_SDP,1 and n_SDP,2 (Figs. <ref>b-f and  <ref>c-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. <cit.>, and we observe that n_SDP,1 is fixed by θ_b for all values of θ_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 θ_b>0^∘ and θ_b<0^∘, we first take the uncertainty in δ, n_SDP,1, and n_SDP,2 into account. We estimate the uncertainty in n_SDP,1 and n_SDP,2 to be ∼±0.02 × 10^12cm^-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 θ_b and θ_t over a range of possible δ (ranging from ∼1.637% to ∼1.654%). This sets a range of possible values of N at each value of θ_t.

Figs. <ref>c-e show dρ/dT as a function of T and n for θ_t = 59.1^∘, -0.7^∘ and 59.4^∘. We observe regions of negative dρ/dT (orange arrows, indicating insulating-like behavior) in the vicinity of the expected values of N. Figs. <ref>f-h show ρ as a function of n at T = 1.7 K for the same values of θ_t. The bounds of the negative dρ/dT regions observed in Figs. <ref>c-e are marked by orange dashed lines. The computed range of possible values of N at each value of θ_t are shown for θ_b<0^∘ (red shaded region) and θ_b>0^∘ (blue shaded region). 

We find that there exists some overlap between the computed values of N and the observed regions of negative dρ/dT, suggesting a possible correspondence with second-order moiré features. However, for fixed δ and n_SDP,1, there appears to be no combination of choices for n_SDP,2 that generate values of N matching negative dρ/dT features for all three configurations (θ_t=59.1^∘, -0.7^∘, and 59.4^∘) simultaneously. Furthermore, there are additional regions of negative dρ/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. <ref>. 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 θ_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. 







 §.§ Room temperature transport in device R2


We measure the room temperature transport of a second rotatable device, R2. Fig. <ref>a-b shows ρ_PDP and ρ_SDP as a function of θ_t. Similar to device R1, we observe extrema in both ρ_PDP and ρ_SDP with 60^∘ periodicity. Following the same procedure, we define θ_t=0^∘ as the  angle corresponding to the maximum values observed in both ρ_PDP and ρ_SDP over a given 120^∘ interval. We observe minima in ρ_PDP as well as suppressed enhancement of ρ_SDP near θ_t=60^∘ in both devices.

Fig. <ref>c plots gate sweeps at various θ_t in device R2. Fig. <ref>d shows a zoomed-in view of the PDP peaks from Fig. <ref>c. We observe a slight suppression in ρ_PDP  very near θ_t=180^∘ 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.

