Dynamic band structure tuning of graphene moiré superlattices with pressure
Matthew Yankowitz Jeil Jung Evan Laksono Nicolas Leconte Bheema L. Chittari K. Watanabe T. Taniguchi Shaffique Adam David Graf Cory R. Dean
cond-mat.mes-hall






 


 
 
 
 

dean@phys.columbia.edu

^1Department of Physics, Columbia University, New York, NY, USA
^2Department of Physics, University of Seoul, Seoul 02504, Korea
^3Centre for Advanced 2D Materials, National University of Singapore, Singapore 117546, Singapore
^4Deparment of Physics, Faculty of Science, National University of Singapore, Singapore 117542, Singapore
^5National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan
^6Yale-NUS College, Singapore 138527, Singapore
^7National High Magnetic Field Laboratory, Tallahassee, FL, USA






Heterostructures of atomically-thin materials have attracted significant interest owing to their ability to host novel electronic properties fundamentally distinct from their constituent layers <cit.>. In the case of graphene on boron nitride, the closely-matched lattices yield a moiré superlattice that modifies the graphene electron dispersion <cit.> and opens gaps both at the primary Dirac point (DP) and the moiré-induced secondary Dirac point (SDP) in the valence band <cit.>. While significant effort has focused on controlling the superlattice period via the rotational stacking order <cit.>, the role played by the magnitude of the interlayer coupling has received comparatively little attention. Here, we modify the interaction between graphene and boron nitride by tuning their separation with hydrostatic pressure.  We observe a dramatic enhancement of the DP gap with increasing pressure, but little change in the SDP gap. Our surprising results identify the critical role played by atomic-scale structural deformations of the graphene lattice and reveal new opportunities for band structure engineering in van der Waals heterostructures.



Dynamic band structure tuning of graphene moiré superlattices with pressure
    Cory R. Dean^1
    December 30, 2023
===========================================================================



Heterostructures fabricated from the mechanical assembly of atomically-thin van der Waals (vdW) crystals represent an exciting new paradigm in materials design. Owing to weak interlayer bonding, 2D crystals with wide ranging characteristics and composition – such as graphene, boron nitride (BN), and the transition metal dichalcogenides – can be readily mixed and matched, without the usual interfacial constraints of conventional crystal growth <cit.>. Moreover, atomic-scale crystalline alignment between the layers often plays a critical role in the resulting device characteristics leading to additional and controllable degrees of freedom. For example, electronic coupling processes that are sensitive to momentum mismatch, such as interlayer tunneling <cit.> or exciton binding <cit.>, can be sensitively tuned by varying the rotational order.  For crystals with closely matched lattice constants, moiré interference at zero-angle alignment can additionally result in long-range superlattice potentials, which in turn can lead to entirely new electronic device characteristics <cit.>. 
  
BN-encapsulated graphene provides a model example of the variety of electronic properties that can be realized on-demand in a single type of vdW heterostructure. At large relative twist angles, BN acts as a featureless dielectric for graphene, minimizing coupling to extrinsic disorder but otherwise remaining effectively inert <cit.>. However, at small twist angle, coupling to the resulting moiré superlattice (MSL) profoundly alters the graphene bandstructure, giving rise to secondary Dirac cones at finite energy <cit.> while also modifying the Fermi velocity near the Dirac point <cit.>. As a consequence several unusual electronic properties have been observed, such as density-dependent topological valley currents at zero magnetic field <cit.> and the fractal Hofstadter butterfly spectrum at high field <cit.> – recently identified to host integer and fractional Chern insulating states <cit.> and charge density waves <cit.>. Additionally, the MSL opens band gaps at both the primary and secondary graphene Dirac points, which are of particular interest for graphene-based digital logic applications.   





Numerous techniques have been developed to control the rotational alignment in graphene/BN and related vdW heterostructures <cit.>. Equally important is the spacing between the layers, which dictates the magnitude of interlayer interactions. However, little  experimental work has been done to characterize or control this parameter, which is often not well-known and considered invariable. Here we demonstrate that by applying hydrostatic pressure to BN-encapsulated graphene, we are able to decrease the interlayer spacing by more than 5%. At small rotation angles the resulting increase in the effective MSL potential substantially modifies the electronic device characteristics. Most dramatically, we observe a divergence in the moiré-induced band gap at the DP with increasing pressure and near doubling over the pressure ranges studied, yielding the largest gap so far demonstrated in pristine monolayer graphene. By contrast, the SDP gap shows little change with pressure. This unexpected result provides new insight into the precise influence of the MSL on the graphene layer and suggests that in addition to electrostatic coupling, lattice-scale deformations play an important role. Our findings reveal that interlayer spacing in vdW heterostructures is both an important and tunable degree of freedom that provides a new route to bandstructure engineering.





Fig. 1a shows a cartoon schematic of our experimental setup. We fabricate BN-encapsulated graphene devices using the vdW assembly technique <cit.> and mount them into a piston-cylinder pressure cell with electrical feedthroughs, capable of reaching temperatures below 1 K and magnetic fields above 18 T (see Methods and Supplementary Information 1 for details). The pressure cell sample space is filled with an oil solution which results in a uniform transfer of pressure to the sample. However, since the Young's modulus in the the out-of-plane stacking direction (c-axis) of the vdW crystal is typically a few orders of magnitude smaller than in-plane <cit.>, the pressure primarily results in a c-axis compression. We first examine the effect of applying pressure on a misaligned heterostructure where no MSL effects are present, and find that the oil environment and application of pressure have virtually no effect on the electronic properties of the graphene. Fig. 1b shows low temperature transport acquired from a misaligned device (> 2^∘ relative alignment) at zero magnetic field in both vacuum (with no oil) and under high pressure (1 GPa). The high-density resistance grows by small amount (less than 25 Ω) with pressure in this device. Other similar devices showed no measurable change or even slightly decreasing resistance under pressure, indicating the field effect mobility of the encapsulated device is largely insensitive to the application of pressure (see Supplementary Information 2). The DP resistance is also not strongly modified, and there is no significant pressure dependence on the contact resistance (Fig. 1b inset). 





Fig. 1c shows similar transport measurements but acquired from an aligned MSL device. Under ambient pressure the device shows excellent transport characteristics with moiré coupling evident by the appearance of two resistance peaks symmetrically located about the DP at roughly ±3.3× 10^12 cm^-2, corresponding to a relative rotation angle of ∼ 0.6^∘. Moreover, the large resistivities exceeding 25 kΩ at the DP and SDP suggest sizable band gaps at those points in the band structure. Notably, as we apply high pressure up to 2.3 GPa, we find that the positions of the secondary peaks move symmetrically towards the DP. 





To understand this effect, we track the back gate capacitance per unit area C_g = en/V_g = ϵϵ_0/t, where e is the charge of the electron, ϵ is dielectric constant of BN, ϵ_0 is the vacuum permittivity, t is the thickness of the BN, V_g is the applied gate bias, and n is the carrier density. We measure the density in two ways; both from the Hall effect and magnetoresistance oscillations (Figs. 2a and b, see Methods). Notably, all devices exhibit a universal increase of C_g with pressure of roughly 9% per GPa independent of their relative alignment (Fig. 2c), which must arise due to a decrease in t and/or an increase in ϵ. To deconvolve the two, we have performed LDA ab initio simulations of bulk BN multilayers under pressure and find approximately 2.5% compression per GPa (green curve in the top inset of Fig. 2c), in remarkable quantitative agreement with previous X-ray diffraction measurements from Ref. <cit.> (reproduced in the red curve). The remaining increase in gate capacitance is therefore attributed to an increase in the dielectric constant of the BN of roughly 6% per GPa (blue curve in the bottom inset of Fig. 2c). Our ab initio simulations also confirm that the BN dielectric constant should grow with pressure, though the effect is predicted to be a few times smaller (green curve). Taken together, this suggests that pressure is able to sensitively tune both the interlayer spacing between layered 2D materials and their dielectric properties <cit.>. Returning to the transport measurements of the MSL device, we find that the three resistance peaks align exactly at all pressures when replotted against charge carrier density n (Fig. 1c inset), suggesting that while the graphene and BN layers move closer together, the relative rotation angle and moiré period remains fixed under pressure.

A second notable feature of the MSL transport is that the resistance at the DP grows strongly with pressure, increasing by roughly 100 kΩ between 0 GPa and 2.3 GPa. We investigate the DP response in more detail by measuring its temperature dependence. Fig. 3a shows an Arrhenius plot of the DP conductivity, σ_DP, versus inverse temperature for various pressures, where for each pressure a linear fit to the simply activated regime (red dashed lines) gives a measure of the activation gap (see Methods). The resulting DP gap, Δ_p, is shown versus pressure in Fig. 3b (square markers). The gap is found to diverge with increasing pressure, and is enhanced by nearly a factor of 2 at the highest pressure studied for the ∼ 0.6^∘ device presented here. Similar gap enhancement was observed in other nearly-aligned samples as well, and we find that the order and direction of pressure cycling do not influence the measurement (see Supplementary Figure 3a). We similarly measure the pressure dependence of the valence band SDP gap, Δ_s, plotted with circle markers in Fig. 3b. In significant contrast to the DP gap, the SDP gap is nearly unresponsive to pressure. The insets of Fig. 3a schematically illustrate the inferred band structure modifications with interlayer spacing, showing a growing Δ_p but a fixed Δ_s. 

As a simple approximation, we may consider the increasing DP gap with pressure to result from increasing MSL coupling owing to the decreasing interlayer spacing. In this case, we expect that the SDP behavior might respond in a similar way, and from this perspective its insensitivity to pressure is surprising. We note, however, that despite considerable effort <cit.> consensus is still lacking as to the exact origin of these band gaps. Lattice scale deformations (in-plane strains and out-of-plane corrugations) of the graphene layer are expected to play an important role <cit.>, but the exact equilibrium structure of graphene in contact with BN remains poorly understood, including whether these deformations even exist in fully BN-encapsulated devices <cit.>. In an effort to understand our observed behavior, we therefore first consider a rigid graphene lattice and examine the effects of pressure on the heterostructure using a combination of ab initio and analytical models (see Supplementary Information 4 - 6). The interlayer electronic coupling between the graphene and BN, V_0 = Ve^-β (z_0-z_r), is highly sensitive to pressure, with ab initio predictions indicating that the average interlayer spacing z_0 should decrease by ∼0.07 Å/GPa. Here V≈ 18 meV is the interlayer electronic coupling, z_r ≈ 3.35 Å  is the equilibrium average interlayer spacing between the graphene and BN, and β≈ 3.2 Å^-1 quantifies the rate of increase of interlayer coupling when the spacing is reduced. Generically, both Δ_p and Δ_s should scale proportionally to V_0 under applied pressure, however this is in stark contrast to our experimental observations where Δ_s in particular exhibits little pressure dependence. Even at 0 GPa, a rigid model does not properly capture the experimentally observed hierarchy of the gaps, predicting a large Δ_s and a Δ_p ≈ 0 (Fig. 4a).

Next we consider the potential role of atomic-scale graphene deformations. By considering realistic values of the graphene elastic deformation potential, our atomic structure models suggest that it is favorable for the graphene lattice to both corrugate out-of-plane and strain in-plane on the moiré scale to minimize the overall stacking potential with the BN substrate (Fig. 4b). These deformations break the sublattice symmetry of the graphene, resulting in a finite average mass term in the Hamiltonian which opens a sizable Δ_p. There are also moiré-induced electrostatic potentials and pseudomagnetic fields which contribute to Δ_s in addition to the mass term. Including corrugations alone opens a small Δ_p at 0 GPa, but still does not recover the observed gap hierarchy of Δ_p > Δ_s (Fig. 4a). By additionally including strains, we are able to recover the observed 0 GPa gap hierarchy (Fig. 4a) as well as the diverging Δ_p and flat Δ_s response with pressure. One such example, whose corrugations and strains are obtained for BN-encapsulated graphene by a model partly informed by ab initio calculation inputs (see Supplementary Information 4 - 6), is shown in Fig. 4c, and predicts gaps in remarkably good quantitative agreement with our experimental data (dashed curves in Fig. 3b).

While more experimental and theoretical work is necessary to understand the exact equilibrium structure of aligned graphene on BN, the varying evolution of the gaps with pressure directly rule out the possibility of a rigid graphene lattice in these structures, and moreover demonstrates that these gaps are of fundamentally different origin as the SDP cone is not a true replica of the DP cone. This suggests the possibility to independently control the magnitude of the two gaps, as well as other features of the moiré band structure and magnetoresponse, by selectively engineering specific lattice deformations. Additionally, sufficiently strong enhancement of the interlayer coupling could drive a phase-transition to the fully commensurate (lattice-matched) stacking configuration <cit.>, marked by the absence of a MSL but the emergence of strong sublattice-symmetry breaking in the graphene and a gap many times room temperature at the Dirac point <cit.>. More generally, our results indicate that a wide variety of vdW heterostructure properties may be tunable by controlling the interlayer coupling strength with pressure. For example, the band structures of Bernal-stacked and twisted bilayer graphene depend critically on interlayer hopping terms <cit.>, as do the strength of proximity-induced spin-orbit interactions in graphene resting on heavy transition metal dichalcogenides <cit.>.



§ METHODS




 §.§ Application of pressure


In order to control the interlayer spacing in vdW heterostructures, we fabricate graphene devices encapsulated in BN with a graphite bottom gate on Si/SiO_2 wafers (Fig. 1a) and make one-dimensional electrical contact using standard reactive ion etching and electron beam patterning and deposition techniques <cit.>. Alignment to one of the encapsulating BN layers is achieved through either optical matching of crystalline edges <cit.> or through thermal self-rotation <cit.>. We then mount the device on a stage with a pre-wired feed-through and affix wires onto the electrodes by hand using silver paint. The stage is then enclosed by a Teflon cup filled with a hydrostatic pressure medium (Daphne 7373 or 7474 oil) <cit.>, and the cup is fit into the inner bore of a piston-cylinder pressure cell and loaded to the desired pressure using a hydraulic press (Fig. 1a). Finally, the pressure cell is affixed to a probe for electrical characterization at low temperature and high magnetic field (see Supplementary Information 1 for further details of the experimental setup). The pressure was determined by measuring the photoluminescence of a ruby crystal at both room- and low-temperature. 



 §.§ Extraction of capacitance and band gaps


We extract n(V_g) for each pressure through either the dispersion of the quantum Hall states in high magnetic field as n = ν e B/h, where ν is the filling factor (Fig. 2a), or from the low field Hall resistance R_xy before the onset of strong Shubnikov-de Haas oscillations as n = 1/eR_xy (Fig. 2b). Quantum Hall states move symmetrically closer to the DP with increasing pressure, and the slope of the Hall resistance decreases with pressure, both implying a growing n(V_g) with increasing pressure. In Fig. 2c, C_g is normalized to its measured value at 0 GPa to account for the different thicknesses of the bottom BNs across the different devices. Additionally, the ab initio value of ϵ at 0 GPa is normalized to match the average experimental value of ≈ 3. 

We extract the band gaps Δ of the DP and SDP at each pressure according to the thermally activated response where σ_DP(T) ∝ e^-Δ/2 k T, where k is the Boltzmann constant.



§ ACKNOWLEDGEMENTS


We thank Pablo San-Jose, Justin Song, Andrey Shytov, Leonid Levitov, John Wallbank, and Pilkyung Moon for valuable theoretical discussions. This work was supported by the National Science Foundation (DMR-1462383). CRD acknowledges partial support from the David and Lucille Packard foundation. We acknowledge Stan Tozer for use of his 16 T PPMS which is partially supported as part of the Center for Actinide Science and Technology (CAST), an Energy Frontier Research Center (EFRC) funded by the Department of Energy, Office of Science, Basic Energy Sciences under Award Number DE-SC0016568. A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-0654118, the State of Florida and the U.S. Department of Energy and additionally provided support for pressure cell development through User Collaboration Grant Program (UCGP) funding. JJ and NL have been supported by the Korean NRF through the grant NRF-2016R1A2B4010105 and the Korean Research Fellowship grant NRF-2016H1D3A1023826, and BLC has been supported by the grant NRF-2017R1D1A1B03035932. EL and SA are supported by the National Research Foundation of Singapore under its Fellowship program (NRF-NRFF2012-01). K.W. and T.T. acknowledge support from the Elemental Strategy Initiative
conducted by the MEXT, Japan and JSPS KAKENHI Grant Numbers JP15K21722.



§ AUTHOR CONTRIBUTIONS


M.Y. and C.R.D. conceived the experiment. M.Y. fabricated the samples, analyzed the data and wrote the paper. M.Y. and D.G. performed the experiments. J.J., E.L. and S.A. developed the theory. N.L. and B.L.C. calculated the ab initio potentials. K.W. and T.T. grew the hBN crystals. C.R.D. advised on the experiments.








 




§ SUPPLEMENTARY INFORMATION




 §.§ Details of pressure experiments


All devices in this study consist of monolayer graphene encapsulated between two layers of boron nitride. The BNs were generally between 20 - 60 nm thick, though our results did not depend on this in any noticeable way. The encapsulated stack sits on a flake of graphite which acts as a local back gate. Fig. <ref>a shows an image of a completed stack on a transfer slide. For aligned samples, the graphene (outlined with a dotted white line) was either intentionally aligned to one of the BNs using straight edges, or rotated to an aligned position with thermal heating during the transfer process. The final stack was partially etched into a Hall bar geometry, leaving some of the bottom BN unetched to prevent the metal contacts from shorting to the graphite gate (Fig. <ref>b). The Hall bar was intentionally kept small (≈ 6 μm by 2 μm) to keep the pressure as uniform as possible across the entire device. The entire device sits on a Si/SiO_2 wafer, which must be diced to approximately 2 mm by 2 mm to fit into the inner bore of the pressure cell (Fig. <ref>c).





To prepare the pressure cell, a clean metal stage (Fig. <ref>d) is threaded with insulated copper 75 μm wires with tinned ends, which are epoxied in place using Stycast 2850 FT using 24LV catalyst (Fig. <ref>e). A ruby crystal is glued to the tip of a thin optical fiber that is fixed in place by the wires for in situ pressure calibration. The sample is then glued above the fiber (Fig. <ref>f). Flexible 15 μm Pt wires are soldered to the copper wires and affixed by hand to the gold sample contacts using Dupong 4929N silver paste (Fig. <ref>g). Next, a Teflon cup is filled with the pressure medium (Daphne 7373 or 7474 oil) and carefully fitted over the sample and onto the stage (Fig. <ref>h). The sample is now completely encapsulated in oil. For the special case of vacuum measurements (0 GPa), no oil is loaded into the Teflon cup. The stage/teflon cup is then fitted into the inner bore of a piston cylinder cell (Fig. <ref>i), and a hydraulic press is used to compress the top of the Teflon cup (Fig. <ref>j). The oil medium hydrostatically increases in pressure as the cup is compressed. A top locking nut is tightened as more force is applied on the hydraulic press until the desired load is reached, at which point the hydraulic press is backed off, and the threads of the locking nut hold the pressure in the cell. The cell is then affixed to the end of a probe (Fig. <ref>k) for low-temperature, high-field measurements. Changing the pressure requires warming the cell to room temperature, adding or removing load with the hydraulic press, and re-cooling the cell.

Although the oil does not influence the electronic properties of the device, special care must be taken to account for its presence during measurements. The oil freezes at around 200 K as it is cooled at ambient pressure, and the freezing point moves to higher temperatures at higher starting loads. The pressure in the cell also drops as the oil cools inside the cell, with a larger relative drop in pressure at smaller initial loads. For example, a starting load below ∼ 0.3 GPa at room temperature will result in nearly ambient conditions at low temperature, while at pressures above 2 GPa there is virtually no change between the room temperature and low temperature pressures. We have found that the primary consequence of this effect is that starting at room temperature loads of roughly 0.5 GPa or less is dangerous for the device, as the stack is typically torn when the oil freezes below these pressures. However, above these pressures the devices always survive the cooling, and we have not noticed any effect of the oil freezing in transport measurements. Finally, special care is taken to account for the large thermal load of the pressure cell when performing temperature sweeps to measure band gaps. The temperature is swept slowly to keep the sample as close to equilibrium as possible, ideally at 0.5 K/min, and no faster than 1.5 K/min.



 §.§ Pressure dependent transport in other devices


The behavior of high density resistance of devices under pressure varies slightly across devices. For Device P2 shown in Fig. 1b of the main text, an increase of 15-25 Ohms is observed at high density under pressure. However, this is not always the case, as is shown for Device P3 in Fig. <ref>. This device exhibits virtually no pressure dependent hole-side resistance, and shows a decrease in the electron-side resistance at high pressure. 





We have measured a total of four gapped devices as a function of pressure. Fig. <ref>a plots the gap of the primary DP, where the data from Device P1 is copied from the main text. Error bars are left off for clarity, but are similar in magnitude to those in Fig. 3b of the main text. Devices P3 and P4 are slightly misaligned, such that the SDPs are outside the accessible density range in the device - therefore the rotation angle is unknown, but must be larger than ∼ 2^∘. Previous work has demonstrated that devices with misalignment angles as large as ∼ 5^∘ can exhibit band gaps in transport <cit.>. In that study, the misalignment angle was determined by scanning tunneling topography measurements, however this is not possible in our devices due to the encapsulating top BN layer. Nevertheless, these devices show clear activation gaps, with insulating behavior at low temperatures and a decrease of the device conductivity by nearly an order of magnitude in the thermally activated regime at high temperatures (Fig. <ref>a - c). The magnitude of the gap grows with pressure in all three samples.





The fourth gapped device (P5) was very well aligned, with an estimated misalignment angle of ∼0.1^∘. However, the device was significantly more disordered than all the other devices in this study, and as a result the onset of the variable range hopping regime occurred at significantly higher temperature. Consequentially, the simply-activated regime spanned only a small range of temperature, and therefore the gap extraction is significantly more uncertain than the other gapped devices examined. Nevertheless, the primary DP showed an unambiguous enhancement with pressure, while the SDP appeared relatively insensitive to pressure (Fig. <ref>b). This behavior is completely consistent with the clean aligned device (P1) presented in the main text, providing further evidence of the fundamental difference between the primary and secondary DPs. Fig. <ref> additionally shows the gap extraction for the SDP of the aligned device P1 in the main text Fig. 3b, demonstrating that this gap does not change with pressure.









The enhancement of the band gap does not depend on the history of the pressure the device has been exposed to, rather it responds directly to the pressure the device is under at the time of measurement. To illustrate this, Fig. <ref>a annotates the order in which the gaps were acquired. In no case does the measured gap fall out of the anticipated sequence. This is an especially robust effect, as it does not depend on whether the device becomes more or less insulating with pressure. Fig. <ref> tracks the CNP conductivity across a number of devices down to the lowest temperatures measured (∼ 2 K). Typically the devices become more insulating at low temperatures for higher pressures, however this is not always the case. For Devices P1 and P3, the low temperature conductivity becomes out of sequence upon unloading pressure, or in reloading pressure after previously pressure cycling and returning to 0 GPa (see the annotations in Fig. <ref>a for the order of curve acquisition). However, even in Device P4, where the pressure was loaded up uniformly, the low temperature conductivity still depends non-monotonically on pressure. This suggests that small changes in the amount of disorder in the device and potentially even the details of its microscopic organization can ultimately influence how insulating the device becomes. However, no matter how the DP conductivity behaves at low temperature, the band gap is always enhanced by pressure, suggesting that the gap enhancement is a robust property of the heterostructure and depends most critically on the interlayer interaction strength between the graphene and the BN rather than small changes in device disorder.



 §.§ Effects of pressure in the quantum Hall regime


The disorder in these devices can also be characterized by considering the amount of density n necessary to switch the Hall resistance R_xy from positive to negative (Fig. <ref>a), as this gives a measure of the effective magnitude of the electron-hole puddles in the bulk. Fig. <ref>b shows that the device disorder is not strongly dependent on pressure, growing by less than the extraction uncertainty over the applied pressure range. Curiously, pressure serves to qualitatively improve the quantum Hall response at high magnetic field. Fig. <ref>c shows one such example, where both integer and fractional quantum Hall states develop more clearly under pressure. This effect is observed in every device examined. It is even more surprising that this improvement persists even after the pressure is released (Fig. <ref>d), as this suggests the need for an explanation describing an irreversible effect. At present, we may only speculate as to the nature of this effect. 





First, it is important to note that despite the high electronic quality of the devices, only the main sequence of IQH gaps at ν = -2, -6, -10, etc., are typically clearly developed in initial measurements in vacuum, with the symmetry broken gaps still developing (Fig. <ref>d). This is not an intrinsic property of the graphene – rather it is symptomatic of the metal contacts sitting above the graphite gate in the device geometry employed here (Fig. <ref>b). While the reason for this is currently not well understood, the situation may be analogous to previous observations in GaAs, where a partial reduction of the carrier concentration in the graphene just in front of the contacts can lead to a depletion region along the boundary, which may impede the ability to observe well-developed quantum Hall plateaus <cit.>. This problem has been addressed previously by leaving Si-gated regions of the graphene Hall bar leads to act as contacts, where the density at the boundary can vary more smoothly <cit.>. However, the devices in this study were intentionally kept small to ensure the pressure is as uniform as possible across the device, and so this geometry was not utilized.

Applying pressure may, for instance, provide a way to reduce this depletion region at the contacts, permitting better coupling to the QH edge modes. This could arise due to self-cleaning of contaminants along the contact boundaries, or if the metal making edge contact to the graphene forms a better bond to the graphene edge and lowers the work function mismatch at the boundary. Both of these effects could in principle persist even after the pressure is released. We find that the contact resistance is not significantly modified by pressure at B = 0 T (Fig. 1b of the main text), pointing to an improvement arising from the detailed electrostatics at the contact barrier in high magnetic field, or to a reduction of non-local contamination which may be most relevant in the QH regime. Further study is necessary to understand the exact nature of this effect.





Finally, we examine the effects of pressure on the Hofstadter butterfly spectrum of the aligned Device P1. Fig. <ref> shows the high field response at ambient conditions (panel a) and at 2.3 GPa (panel b). As pressure enhances the effective strength of the moiré potential, it should also influence the relative size of the LL gaps originating from the secondary DPs. The qualitative behavior of the two butterfly maps seems to be in qualitative agreement with this expectation, with stronger features originating from the secondary DPs under pressure. However, this must be deconvolved with the overall change in the QHE behavior with pressure, and a full investigation of this effect is outside the scope of this work. 



 §.§ First principle calculations and parametrizations


The atomic and electronic structures of graphene and hexagonal boron nitride (G/BN) interfaces are captured using input from ab-initio density functional theory (DFT) calculations.
We rely on exact exchange and random phase approximation (EXX+RPA) input for atomic structure, while our treatment of the electronic structure is based on local-density approximation (LDA). 



  §.§.§ Calculation of the dielectric constant from DFT

One of the striking features of the experiment is an increase of the dielectric constant as a function of pressure. 
Our DFT calculations for bulk confirm this behavior showing increases on the order of 3% for the applied pressures up to 2.5 GPa as shown in Fig. <ref>.

When the pressure is modified from P=0 to P=2.5 GPa,
the EXX+RPA equilibrium interlayer distance is squeezed by ∼ 0.2 as we will discuss in the following subsection.







First principles calculations based on density functional theory (DFT) were performed using Quantum Espresso <cit.> under LDA  
within the plane wave basis set. Norm-conserving Vanderbilt pseudopotentials were used. 
The structures are fully optimized without any 
symmetry constraint by using conjugate gradient method.
The convergence criterion for the force on each ion is
taken to be less than 0.005 eV/ while the energy is converged with a tolerance of 10^-5 eV.


In an insulating BN, the optical properties are dominated with the direct inter band contributions to the absorptive or imaginary part of dielectric function <cit.>. The matrix elements for a given interband transition β→α for a set of plane wave Bloch function: 



    | ψ_k,n⟩  = e^iG · r u_k,n  = 1/√(V)∑_G a_n,k,G e^i(k + G) · r


is given by:


    M̂_α ,β  = ( ∑_G a_n,k,G^* a_n',k,G G_α)( ∑_G a_n,k,G^* a_n',k,G G_β)



These matrix elements account only for the interband transitions, i.e the electric-dipole approximation where the momentum transfer is zero. The imaginary part of the dielectric tensor ε_2α ,β is a response function that comes from a perturbation theory within adiabatic approximation. All of the possible transitions from the occupied to the unoccupied states without local field effects are given by the imaginary dielectric function:

 

    ε _2α ,β  = 4π e^2 /Ω N_k m^2 ∑_n,n'∑_k M̂_α ,β/(E_k,n'  - E_k,n )^2 {f(E_k,n )/E_k,n'  - E_k,n  + ħω  + iħΓ + f(E_k,n )/E_k,n'  - E_k,n  - ħω  - iħΓ},



where Γ is the adiabatic parameter which must be zero due to conservation of energy. So, the equation can be rewritten in terms of Dirac delta functions: 



    ε_2α ,β  = 4π e^2 /Ω N_k m^2 ∑_n,n'∑_k M̂_α ,β f(E_k,n )/(E_k,n'  - E_k,n )^2 [ δ (E_k,n'  - E_k,n  + ħω ) + δ (E_k,n'  - E_k,n  - ħω )],



However, interaction with the electromagnetic field which also take place even in the absence of photons, i.e spontaneous emission, brings intrinsic broadening to all excited states that leads to a finite life time (Γ > 0). In the limit of non-vanishing Γ, the dielectric tensor takes the Drude-Lorentz form <cit.>:
 


    ε _2α ,β  = 4π e^2 /Ω N_k m^2 ∑_n,kdf(E_k,n )/dE_k,nηωM̂_α ,β/ω^4  - η ^2 ω ^2  + 8π e^2 /Ω N_k m^2 ∑_n  n'∑_k M̂_α ,β/E_k,n'  - E_k,nΓω f(E_k,n )/[(ω _k,n'  - ω_k,n )^2  - ω^2 ]^2  + Γ^2 ω ^2 ,



and the real part of the dielectric function is derived from the Kramers-Kronig transformation:



    ε_1α ,β  = 1 + 2/π∫_0^∞ω 'ε _2α ,β (ω ')/ω'^2  - ω ^2 dω'.



Calculations of the dielectric function were carried out using the EPSILON package 
distributed with Quantum Espresso <cit.>.


















  §.§.§ DFT deformation with pressure




The mechanical behavior of nearly aligned graphene and boron nitride heterostructure (G/BN) is influenced by the encapsulating top BN layer.
The interlayer interaction of nearly aligned G/BN is described using the assumption that they can be viewed as a collection of different stacking registries to be treated independently, from which the atomic structure of the whole sample can be obtained. 
We obtain the equilibrium spacing for each stacking point z_𝐬 (𝐬∈{AA, AB, BA}, see Fig. <ref>) by minimizing the total potential energy as a function of the
separation distance between the top and bottom BN layer d. We also assume that the elastic resistance of the graphene sheet is negligible in the out-of-plane direction. 
The interlayer potential energy of the graphene sheet considering the effects of bottom and top BN layers is given by:



    V(z_AA,z_AB,z_BA,d) = 1/3∑_𝐬[V_𝐬(z_𝐬) + V̅(d - z_𝐬)]



where V̅(d-z) = 1/3∑_𝐬 V_𝐬(d-z) accounts for the averaging of the interlayer potential between the top BN and the graphene sheet due to their misalignment.

We use for the interlayer interaction potentials the exact exchange and random phase approximation (EXX+RPA) for G/BN and BN/BN interactions as implemented in VASP. <cit.>


In equilibrium, the local pressure exerted by the encapsulating layer is equally balanced by that coming from the underlying substrate, i.e. P_𝐬(z_𝐬) = P̅(d - z_𝐬). Defining the interlayer distance between graphene and the top layer z'_𝐬 = d - z_𝐬, we can write the equilibrium condition as 1/3∑_𝐬 P_𝐬(z_𝐬) = 1/3∑_𝐬P̅(z'_𝐬) = P, where P is the external pressure applied on the system. Fig. <ref> illustrates how the corrugation profiles of a graphene sheet change with pressure under the influence of two BN layers.









  §.§.§ Electronic coupling parameters and elastic properties of G/BN


In our treatment of the moiré bands in G/BN, we consider modifications in the low-energy regime of graphene electronic structure which are introduced by the interlayer tunneling. The modification is captured by two interlayer coupling terms V and Ṽ whose values can be obtained from first-principle calculations. V and Ṽ can be calculated from the following expressions <cit.>:  



    V = 1/2[t_NC^2/|ϵ_N|-t_BC^2/|ϵ_B|]≈ 0.04  meV,Ṽ = √(3)/2[t_NC^2/|ϵ_N|+t_BC^2/|ϵ_B|]≈ 10  meV,
 


where t_NC (t_BC) stands for the interlayer hopping term between C–N (C–B) atoms which highly depends on the interlayer separation, and ε_N (ε_B) are the on-site energies of N (B) atoms. Both parameters are calculated from the parametrizations of the interlayer terms in Ref. <cit.>, at a separation of z_r = 3.35. 
To fully account for the modification in the electronic couplings, ab-initio method in Ref. <cit.> introduces three exponential factors {β_AA,β_AB,β_BB}, each quantifies the rate of increase of each sublattice term (β_i > 0). The values of β_i's are determined by the amount overlap between carbon orbitals with the orbitals of the underlying BN atoms, for which it ranges between 3.0–3.3 ^-1 <cit.>.

As sample deformation is expected to play an important role in determining the band features, we model the structural relaxation within Born-von Karman plate theory in which the elastic properties of graphene are fully characterized by two Lamé constants: λ_g = 3.25  eV ^-2 and μ_g = 9.57  eV ^-2 <cit.>. One mechanism through which deformation leaves an imprint on the bands is the intralayer effects which are induced by the change of the carbon on-site energy and nearest-neighbor hoppings as each carbon atom is displaced with respect to its neighbors. They give rise to additional terms in the moiré couplings which are equivalent to the electrostatic potentials and the pseudomagnetic fields, quantified by γ_V≈ 4.0 eV and γ_B≈ 4.5 eV respectively <cit.>.



 §.§ Theoretical Model


Graphene and BN are two-dimensional materials with hexagonal crystal structures, but with different lattice constants. When graphene is deposited on a BN substrate with near perfect alignment (small twist angles), the slight lattice mismatch  (ε = a_G-a_BN/a_BN = -1.7%), leads to the formation of a moiré structure which is characterized by a wavelength that is one or two orders of magnitude larger than graphene's unit cell. As a consequence, the dynamics of the low energy Dirac electrons of graphene is governed by the corresponding long-wavelength component of the electronic couplings. Similarly, the interlayer potentials and deformations in G/BN are also largely characterized by this long-wavelength theory.  It is by now established that the dominant contributions are captured by the first harmonics of the moiré structure
characterized by a set of vectors G⃗_m, that is related with the original graphene's lattice vectors g⃗_m:  



    g⃗_m = 
    R̂_2π (m-1)/6 (0,g) ,G⃗_m = [(1+ε)-R̂_θ]g⃗_m≈εg⃗_m-θẑ×g⃗_m,



where m ∈{1,2,...,6}, R̂_θ denotes a rotation by θ, g = 4π/3a is the length of graphene lattice vectors with a≈ 1.42 stands for the carbon-carbon distance, and the approximate sign indicates the approximation within small twist angle limit (θ≪ 1). 
We also show in Fig. <ref>, that G⃗_m defines the moiré Brillouin Zone (mBZ), whose lateral dimension is scaled by ε̃=√(ε^2+θ^2)  with respect to graphene BZ that amounts to values ≲ 5% for twist angles θ≲ 2^∘. In G/BN system which possesses triangular symmetry, we can define two periodic functions within the first harmonics: (1) f_1(r⃗) = ∑_mexp(iG⃗_m·r⃗) which satisfies inversion and hexagonal symmetries, and (2) f_2(r⃗) = -i∑_m(-1)^mexp(iG⃗_m·r⃗) which is asymmetric under inversion.





  §.§.§ Definition of symmetric and antisymmetric deformations




Deformation of the sample results from the minimization of the energy functional, leading to the following equation of motion for the in-plane components:



    ε+θ(ẑ×1)/ε̃^2A_g ∇⃗[U_1f_1(r⃗) + U_2 f_2(r⃗)]
    ≈ 2μ_g[∇⃗^2u⃗-1/2(ẑ×∇⃗)
    (∂_xu_y-∂_yu_x)]
    +λ_g∇⃗(∇⃗·u⃗),



Similar to the electronic couplings, the interlayer potential consists of symmmetric and asymmetric components which are quantified by U_1 and U_2. We treat the deformation within the first harmonics, which dominates the relaxation profile for small atomic displacements (|u⃗|≪ a) which accounts for the approximate sign in <ref>. Since the interlayer potential is completely defined by the first harmonics (Eq. <ref>, LHS), the solution for the displacement vector u⃗(r⃗) = u_x x̂ + u_y ŷ + h(r⃗) ẑ is also constrained by moiré periodicity, containing only the first harmonics components. The displacement and the corresponding in-plane expansion of carbon site at r⃗ is thus given by:



    u⃗(r⃗) =
    k_s|ε^3|/G^2ε̃^2 [∇⃗+χ_Rθ/ε(ẑ×∇⃗)][U_1f_1(r⃗)+U_2f_2(r⃗)] + [z_1f_1(r⃗) + z_2f_2(r⃗)] ẑ, 
    
    a/a_0≈ 1+1/2Tr[u_ij(r⃗)] = 1 + 1/2[∂ u_x/∂ x+∂ u_y/∂ y] = 
    1 - 
    k_s|ε^3|/2ε̃^2 [U_1f_1(r⃗)+U_2f_2(r⃗)],



where U_1 (U_2) parametrize the inversion-symmetric (asymmetric) components of in-plane deformation, and z_1 (z_2) quantifies similar components for the deformation in the out-of-plane direction, i.e. corrugation. Meanwhile, the coefficient k_s which relates u⃗(r⃗) with the potentials, and the factor χ_R which sets the magnitude of the curl-like terms are as follows 



    k_s = 2/3√(3)ε^2a^2(λ_g+2μ_g)
    = 0.029
    
      meV^-1,
    χ_R = [2+λ_g/μ_g]
    ≈ 2.34,



It can be seen in the the subsequent analysis that k_s serves as the coefficient which encodes the contribution of in-plane deformation which is parametrized by U_1 and U_2, on the resulting gaps. We also assume the following definition for the symmetric and asymmetric components of the deformation:



    [ U_0; U_1; U_2 ]
    = M[ U_AA; U_AB; U_BA ]
    = 1/18[ 6(U_AA+U_AB+U_BA);   2U_AA-U_AB-U_BA;   √(3)(U_AB-U_BA) ],



in which we take AA stacking point as the origin of our coordinate system, and the matrix elements of M are obtained in a way consistent with our definitions of f_1(r⃗) and f_2(r⃗), and in accordance with the convention in the literature <cit.>. The symmetric and antisymmetric profiles are plotted in Fig. <ref>. We also use similar definition of symmetric and antisymmetric components for the corrugation and moiré electronic couplings. Within this choice of coordinates, we can expect from symmetry consideration alone that a sizable primary gap only opens up when there is a relatively large symmetric component in the deformation (the asymmetric terms vanish after doing the mBZ average).  In general, since in our choice of coordinate system the inversion-asymmetric electronic couplings dominate (Ṽ≫ V), this implies that the dominant contribution to the  global mass term results from symmetric deformations.  Since the secondary Dirac cone gap is not a simple average over the mBZ, we expect both the symmetric and asymmetric contributions to be equally important.



  §.§.§ Moiré couplings and gaps in G/BN


The dynamics of low-energy electrons in G/BN can be understood as a relativistic Dirac particle subject to three types of perturbations: (1) H_0(r⃗), describing the periodic electrostatic potential, (2) H_z(r⃗), acting as a local mass term which arises from the sublattice asymmetry, and (3) H_xy(r⃗), quantifying the gauge field due to asymmetry in hoppings induced by the interaction with the BN substrate. For graphene electrons on valley K, the equation of motion can be expressed as follows



    H_K(r⃗) 
       =υp⃗·σ⃗ +
    
    
    
    H_0(r⃗)1 +  H_z(r⃗)σ_3+ [ε(ẑ×1)-θ/ε̃G]∇⃗H_xy(r⃗)·σ⃗,



It is important to note that only the H_xy contribution is scaled by factors which depend on the twist angle. The term in Eq. <ref> which is proportional to θ can be gauged away so that the influence of H_xy on the electronic structures is scaled by χ_θ = ε/ε̃. The three pseudospin components can be written in a compact manner using the following matrix notation:



    [  H_0(r⃗) H_xy(r⃗)  H_z(r⃗) ]
    =
    [       1 f_1(r⃗) f_2(r⃗) ]W,





where W stands for the electronic coupling matrix. Typically each pseudospin term for zero pressure has oscillating amplitudes on the order of ∼50 meV. Around AA, DFT parametrizations of the interlayer couplings lead to a dominating inversion-asymmetric electronic coupling (Ṽ≫ V). In addition, the exponential coefficients for all psuedospin terms also agree within 10%. Therefore, in this work we develop an analytical model containing only the inversion-asymmetric term Ṽ, which is scaled according to a single exponential coefficient β when pressure is applied on the system. Treating the problem at θ = 0^∘, W is given by the following expression:



    W = 
    
    
    
    
    
    
    [ 1+
    k_s[    0 6U_1 6U_2;    0  U_1 -U_2;    0 -U_2 -U_1 ]]
    CW_0
    +
    [        0        0        0;  -Γ_VU_1 -Γ_B U_2        0;  -Γ_VU_2  Γ_B U_1        0 ],



which is related to the corrugation matrix C and the coupling matrix in a rigid sample W_0:



    C = 
    M^-1[     e^-3β z_1             0             0;             0 e^-2√(3)β z_2             0;             0             0  e^2√(3)β z_2;               ]M,
    W_0
    =
    e^-β (z_0-z_r)[     0     0     0;     0     2 -√(3);     1     0     0 ]Ṽ/2,



The coupling terms appearing in W_0 are based on the special ratio between psuedospin components <cit.>, while the prefactors for z_1 and z_2 are estimated from the full numerical treatment of the Fourier components within the first harmonics. In Eq. <ref>, the first change in the electronic couplings as captured by the first term, is induced by shifts in the positions of the carbon atoms with respect to the underlying BN substrates. In-plane deformation of graphene brings additional contributions to the electrostatic potentials (Γ_V) and pseudomagnetic fields (Γ_B), both independent on the interlayer spacings.

The primary gap arises primarily from the mass term Δ_p∼ 2|ω_z|, while the the secondary gap that appears on the valence band on valley K' at the edges of the mBZ (see Fig. <ref>): K⃗' = 1/3(G⃗_2 + G⃗_3) = 1/3(G⃗_4 + G⃗_5) = 1/3(G⃗_1 + G⃗_6) results from an interplay between the pseudospin terms H_0, H_z and H_xy <cit.>.  These gaps have the following analytical form in a fully relaxed sample:



    Δ_p   
    = 2√(3)Ṽe^-β (z_0-z_r)[cosh(2√(3)β z_2)-e^-3β z_1 + k_s(2U_1e^-3β z_1+U_1cosh(2√(3)β z_2) + √(3)U_2sinh(2√(3)β z_2))], 
       Δ_s   
    = √(3)Ṽ/6e^-β (z_0-z_r)[7e^-3β z_1 + 2cosh(2√(3)β z_2) + 12sinh(2√(3)β z_2) ] - Γ_G(√(3)U_1-U_2) +δ_s,



where  δ_s is a small additional contribution to the secondary gap which results from the interplay between corrugation and in-plane deformation (assumed here to be zero).  We define Γ_G = (√(3)/2)(Γ_V + 2Γ_B), for Γ_V,B= |ε^3|k_sγ_V,B/ε̃^2.  The first principle calculations at θ = 0^∘ discussed above, gives a sizable Γ_G≈ 5.6  meV^-1.  This is important, because for generic band parameters, in the case of Γ_G→ 0 the bands that form the sDC are possibly mixed with the surrounding moiré minibands such that no gap could be observed experimentally.  Since graphene sites have an affinity towards the BA configuration, where one carbon atom sits on top of the boron atom, both U_1 and U_2 are expected to acquire positive values ∼3 meV, such that the BA domain is enlarged compared to the rigid case. 

The exponential prefactor e^-β (z_0-z_r) in the primary gap gives the increase in gap arising from the increased coupling between the graphene and BN when the two layers are closer together under pressure.  This multiplies two terms, the first depends only on the corrugations and to leading order is given by the symmetric component of the out-of-plane deformations ∼  3 β z_1.  The second term proportional to k_s depends on both in-plane and out-of-plane deformations, but for realistic parameters depends only on the symmetric deformation ∼ 3 k_s U_1.  This shows that to leading order, the primary gap is determined mostly by the symmetric in-plane and out-of-plane deformations and vanishes in the case of rigid graphene.

The secondary gap is finite even for the rigid case without any deformations Δ_s = (√(27)/2) Ṽ e^-β (z_0-z_r) that grows with increasing pressure.  The experimental observation that the primary gap is larger than the secondary gap even at zero pressure implies that the second (negative) term in Eq. <ref> proportional to Γ_G must be of comparable magnitude.  In contrast to the primary gap, the secondary gap is controlled by both the symmetric and antisymmetric components of the deformation i.e. ∼ (U_2 - √(3)U_1) for in-plane deformations, and ∼ ((8 √(3)/7) z_2 - z_1) for corrugations.  In both cases, for the secondary gap to remain constant under pressure, the difference between the asymmetric deformations and symmetric deformations must increase with pressure.       

We note that these experimental observations consistent with DFT expectations.  The EXX-RPA calculations discussed above show that with increasing pressure, (i) The difference between the asymmetric and symmetric components of in-plane deformations increases, while (ii) The asymmetric part of out-of-plane corrugation decreases.  This suggests that the dominant contribution to the flatness of the secondary gap comes from corrugations.  Physically, the AA points (carbon-carbon on BN) start further away form the underlying BN substrate. With increased pressure, they compress more than either the AB (carbon on nitrogen) or BA (carbon on boron).  While all potentials increase with pressure, the BA increases the least making this region expand.  This naturally explains the growing primary gap.  Similarly, under pressure the G/BN spacing decrease at all stacking points (so that the graphene becomes flatter). But the AA points which were further away to begin with decrease faster than either the BA or AB points, which is consistent with the Δ_s being constant.




 §.§ Reverse Engineering




  §.§.§ Fixed corrugation, free in-plane deformation
 
Since the observed secondary gaps are smaller than the primary gaps at all pressures, this implies that Γ_G term originating from the in-plane strain-induced pseudomagnetic fields and potentials is large. 

We rely on DFT (EXX+RPA) calculations to determine the corrugation, assuming that the forces exerted by a top BN layer completely fixes the shapes of graphene corrugations under pressure. Since the corrugation is known at all pressures, the in-plane deformation parameters U_1 and U_2 are then uniquely determined by the observed primary and secondary gaps.  Here, the graphene sample is relatively flat δ z ≲ 0.1 and the insensitivity of the secondary gap is accounted by an increase in U_1 with respect to U_2 as shown in the main text and in the top row of Fig. <ref>.





  §.§.§ Fixed in-plane deformation, free corrugation
 

To explore the role of corrugation in the experimental observation, we consider the opposite scenario where the in-plane deformation is assumed to be constant under pressure.  We fix U_1 = 3 meV and U_2 = 2 meV, comparable to ab-initio estimates <cit.>. The insensitivity of the secondary gap under pressure constrains both the symmetric (z_1) and asymmetric components (z_2) of the corrugation, resulting in a significant transformation of the corrugation profile as shown in the middle row Fig. <ref>.  One striking feature is that the experimental observations can be explained by a decrease in the asymmetric corrugations i.e. the separation between AB and BA becomes smaller under pressure.  Since this contradicts DFT expectations we think this scenario is unlikely.



  §.§.§ Self-consistent corrugation and in-plane deformation


Another possible approach to address the gap behaviors is to construct interlayer potentials on G/BN special stacking points, which simultaneously define in-plane and out-of-plane deformation in G/BN and capture the increasing and constant trend in the primary and secondary gaps respectively. We assume a uniform pressure across the sample, and it is expected that there exists a certain potential profile which can account for the observed gaps. Inspired by the DFT results which suggest Lennard-Jones (LJ)-like potential profiles around the minima, we limit our search for the self-consistent potentials U_𝐬 which assume the LJ form at each stacking point (𝐬∈{AA,AB,BA}), which also leads to an interlayer distance of z_𝐬(P) at pressure P:



    U_𝐬(z) =
    σ_𝐬[
    1/z^2n_𝐬-1-
    (2n_𝐬-1/n_𝐬-1)
    1/z_0,𝐬^n_𝐬 z^n_𝐬-1]
    + U_0,𝐬,
    
    
    z_𝐬(P) = [1/2z_0,𝐬^n_𝐬 + √(1/2z_0,𝐬^n_𝐬+4PA_g/(2n_𝐬-1)σ_𝐬) ]^1/n_𝐬



The coefficient σ_𝐬 accounts for the strength of the potentials, and z_0,𝐬 denotes the equilibrium distance at zero pressure. The power coefficient n_𝐬 determines the behaviors of the profile on the small and large z tails, while A_g denotes the area of graphene's unit cell. We therefore try to obtain the best estimates for the potentials which account for the observed experimental trend. The LJ-like potentials which lead to remarkable agreement with the experiment can be found, i.e. <10% error for each data point, with the fitting parameters are shown in Table <ref>. 

As expected, this solution suggests the need for asymmetric components of the corrugation to decrease rapidly under pressure.  However, this leads to a striking conclusion: Although BA is initially at a minimum spacing at zero pressure, AB overtakes BA as the closest point from the underlying BN substrate at P ≳ 1.5 GPa. This is accompnied by a slight increase in the asymmetric component of the in-plane deformation. Again, in disagreement with the DFT, the self-consistent solution is dominated by the rapid change in the asymmetric component of the corrugations (z_AB-z_BA).  The EXX-RPA ab initio potentials (which can also be fit with Lennard-Jones-like functions) are shown in the inset for comparison with the potentials determined self-consistently.



  §.§.§ Generic features of the deformation




While the analysis was done in the simplified semi-analytical model, we have checked that the conclusions remain robust in a full numerical simulation of the first-Harmonic approximation.  The relatively large primary gap suggests there is a strong presence of either in-plane relaxation or corrugation.   The fact that secondary gap is always smaller than the primary gap requires large in-plane deformation at all pressures. Since the increase in the moiré couplings as the layers get closer would give larger secondary gaps, the experimental observation of a flat secondary gap behavior implies that one or both of the following must be happening: (1) There is a significant decrease in the corrugation asymmetry z_2, which places AB closer to BA in the out-of-plane deformation when pressure is applied on G/BN, (2) There is a rapid transformation of the in-plane deformation profile under pressure such that the difference between U_1 and U_2 balances the effects of stronger electronic couplings. 

