Elastic screening of pseudo gauge fields in graphene
Christophe De Beule Robin Smeyers Wilson Nieto Luna E. J. Mele Lucian Covaci





Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA

Department of Physics and NANOlab Center of Excellence, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium

Department of Physics and NANOlab Center of Excellence, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium

Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, Pennsylvania 19104, USA

Department of Physics and NANOlab Center of Excellence, University of Antwerp, Groenenborgerlaan 171, 2020 Antwerp, Belgium




§ ABSTRACT 
Lattice deformations in graphene couple to the low-energy electronic degrees of freedom as pseudo scalar and gauge fields. Using molecular dynamics simulations, we show that the optical component of the displacement field, i.e., the relative motion of different sublattices, contributes at equal order as the acoustic component and effectively screens the pseudo gauge fields. In particular, we consider twisted bilayer graphene and corrugated monolayer graphene. In both cases, optical lattice displacement significantly reduce the overall magnitude of the pseudo magnetic fields. For corrugated graphene, we further show that the corresponding electronic bands are significantly modified by the screened pseudo magnetic field. Previous studies based on continuum elasticity, which ignores this effect, have therefore systematically overestimated the strength of the strain-induced pseudo magnetic field. Our results have important consequences for the interpretation of experiments and design of straintronic applications.


Elastic screening of pseudo gauge fields in graphene
    Lucian Covaci
    September 9, 2024
====================================================



It is well known that lattice deformations in graphene couple to the low-energy electronic degrees of freedom as effective scalar and gauge fields <cit.>. Intuitively, this can be understood from local symmetry breaking which shifts the Dirac cones near charge neutrality both in energy and momentum, respectively. Indeed, the microscopic 𝒞_3z symmetry of pristine graphene pins the two Dirac points at the zone corners (valleys) of the Brillouin zone <cit.>. Atomic displacements that break this symmetry, i.e., shear strain, therefore result in a spatially-varying shift of the Dirac point. This is the action of a pseudo vector potential with opposite sign in the two valleys which preserves time-reversal symmetry. For specific strain configurations, the corresponding pseudo magnetic fields give rise to pseudo Landau levels <cit.> with field strengths that can exceed several hundreds of Tesla <cit.>. Moreover, similar pseudo gauge fields also arise in strained 2D semiconductors <cit.>, as well as in 3D topological semimetals <cit.>. In fact, pseudo gauge fields were first considered in semiconductors in the 1980s <cit.>. In the context of graphene, they were first discussed in carbon nanotubes where the curvature of the tube gives rise to a pseudo gauge field that results in a band gap for nanotubes that should otherwise be metallic <cit.>.

In most electronic continuum theories, the pseudo gauge field is derived from deformations obtained from continuum elasticity which accounts for acoustic displacements. Since graphene has two sublattices [see Fig. <ref>(a)], there are both center-of-mass (acoustic) and relative (optical) displacements, and both are important in the long-wavelength limit. In this paper, we show that the optical displacement field contributes at the same order of magnitude to the pseudo vector potential as the acoustic shear strains. This finding resolves a long-standing conundrum from earlier works that combine molecular dynamics (MD) simulations and electronic tight-binding models, for which the pseudo magnetic field is found to be much smaller as predicted by elastic theory <cit.>.

We demonstrate our theory for two graphene systems: twisted bilayer graphene and corrugated graphene, see Fig. <ref>(c). The displacement fields are computed from MD simulations with the lammps code <cit.>. In particular, for twisted bilayer graphene near the magic angle we show that optical displacements induced by lattice relaxation reduce the pseudo magnetic field by almost one order of magnitude. 


Pseudo gauge fields in graphene revisited.—Consider a sheet of graphene subject to atomic displacement fields that vary slowly relative to the lattice. These may be induced by external stresses or lattice relaxation, e.g., in a graphene moiré. In the long-wavelength limit, we define fields u_σ( r) and h_σ( r) for sublattice σ = A,B [see Fig.<ref>(a)] projected on and normal to the nominal graphene xy plane, respectively. The long wavelength acoustic and optical displacement fields are defined as

    u     = ( u_A + u_B ) / 2,            h = ( h_A + h_B ) / 2, 

    v     = u_A - u_B,            w = h_A - h_B. 


Atomic displacements modulate the electronic hopping amplitude, which couples to the low-energy Dirac electrons as effective scalar and gauge fields <cit.>. We revisit this theory, starting from a tight-binding description for the p_z electrons of graphene in the nearest-neighbor approximation. We show that there is an important contribution to the pseudo gauge field from optical displacements [v and w in Eq. (<ref>)] that was not considered in previous works and which cannot be obtained from continuum elasticity <cit.>. Here we do not consider the pseudo scalar field or deformation potential <cit.> which to leading order only depends on acoustic displacements as it enters through on-site potentials. Moreover, the deformation potential gives rise to electron-hole puddles and is strongly screened <cit.>.

In the presence of atomic displacements, the Hamiltonian for p_z electrons in graphene can be written as

    H = -∑_ r∑_n=1^3 t_n( r) c_A^†( r) c_B[  r + δ_n( r) ] + h.c.,

where the first sum runs over cells r and the second one over nearest neighbors. Here c_σ^†( r) [c_σ( r)] are electron creation (annihilation) operators. The position of A atoms is given by r +  u_A( r) + h_A( r) ẑ and δ_n( r) = δ_n^0 +  u_B( r + δ_n^0) -  u_A( r) + [ h_B( r + δ_n^0) - h_A( r) ] ẑ are the nearest-neighbor bond vectors with δ_n^0 those of pristine graphene, see Fig. <ref>(a). Taking the continuum limit, one obtains an effective low-energy Hamiltonian <cit.>

    H_eff = ħ v_F ∑_τ∫ d^2  r  ψ_τ^†( r) [ -i ∇ + eτ/ħ A( r) ] ·σψ_τ( r),

with ħ v_F = √(3)t_0a/2, σ = (τσ_x, σ_y), and field operators ψ_τ( r) = [ ψ_τ A( r), ψ_τ B( r) ]^t where τ = ± 1 is the valley index. See Supplementary Material (SM) for details [See Supplementary Material at [insert link] for more details on the derivation of the continuum model and the pseudo gauge field, the lammps calculations for twisted bilayer graphene and corrugated graphene, as well as continuum elasticity for the corrugated graphene]. Here we take the zigzag direction along the x axis. In this case, the pseudo vector potential is defined as

    A_x( r) - i A_y( r) = -1/ev_F∑_n=1^3 δ t_n( r) e^i K ·δ_n^0,

where δ t_n = t_n - t_0 is the change in hopping with t_0 ≈ 2.8  eV and K_+ = 4π/(3a) x̂ is the zone corner [see Fig. <ref>(b)] with a ≈ 2.46  Å the graphene lattice constant <cit.>. In lowest order of displacements and their gradients,

    δ_n( r)- δ_n^0 ≃( δ_n^0 ·∇) (  u+ h ẑ) - (  v + w ẑ),

where the first and second term give acoustic and optical contributions, respectively.
From Eqs. (<ref>) and (<ref>),

    A = √(3)ħβ/2ea[ [ u_yy - u_xx; u_xy + u_yx ] + 2 √(3)/aẑ×(  v + w ∇ h ) ],

with β = -a/√(3)t_0. ∂ t/∂ d|_0 ≈ 3 
<cit.> 
and u_ij = [∂_i u_j + ∂_j u_i + ( ∂_i h ) ( ∂_ j h ) ]/2
the strain tensor for acoustic displacements. Note this expression only holds for the orientation of the unstrained graphene shown in Fig. <ref>(a). See SM for the general case <cit.>. The first term in Eq. (<ref>) is the usual acoustic contribution (A_ac) while the second term (A_op) is new and gives a contribution from the relative motion between sublattices. This is one of the main results of our work. It is reminiscent of the renormalization of the electron-phonon coupling by optical phonons in pristine graphene <cit.>. In general, the optical component gives a pseudo magnetic field (PMF)

    B_op = ẑ·( ∇× A_op) = 3 ħβ/ea^2∇·(  v + w ∇ h ).


To demonstrate our theory, we performed lammps molecular dynamics simulations for two graphene systems: moiré graphene and corrugated graphene. In all cases, we find that A_op acts to reduce the overall magnitude of the PMF.

Graphene moirés.—We first consider twisted bilayer graphene (TBG) <cit.>, a twist moiré formed by stacking two layers of graphene with a relative twist, see Fig. <ref>(c). In moirés, the atomic stacking between layers varies spatially. Certain stackings are favorable, and the systems relaxes to minimize the total elastic and adhesion energy. This gives rise to atomic displacements and concomitant pseudo gauge fields <cit.>.

As before, one defines displacement fields projected on the xy plane, u_l and v_l with l = 1,2 the layer index, and similar for out-of-plane displacements. These displacements are calculated for a relaxed structure with lammps using the REBO intralayer <cit.> and DRIP interlayer potential <cit.>. As expected for a twist moiré, the acoustic in-plane displacement field is almost entirely solenoidal, consistent with previous theory <cit.> and experiment <cit.>. On the other hand, the optical field is mostly irrotational yielding a finite optical PMF. See SM for details on the displacement fields <cit.>. The acoustic, optical, and total PMF are shown in Fig. <ref>(a). We only show the PMF for one layer since the D_6 symmetry of the moiré yields A_2(x,y) = diag(-1,1)  A_1(-x,y) from 𝒞_2y rotation symmetry. Moreover, the PMF is odd under r ↦ - r due to 𝒞_2z symmetry. Consequently, we have B_2(x,y) = B_1(x,-y). We find that the acoustic and optical PMFs have similar shapes but opposite signs such that B_tot = B_ac + B_op is about five times smaller than B_ac. We quantify this by plotting the root mean square (RMS) as a function of twist angle in Fig. <ref>(b). The magnitudes of both B_ac and B_tot are nearly constant above the magic angle <cit.> and the ratio √(<B_ac^2>/<B_tot^2>)≈ 4.8 for twists from 1^∘ to 4^∘. We also show the PMF from Eq. (<ref>) in Fig. <ref>(c), using the atomic positions from lammps and intralayer hopping

    t(d) = V_ppπexp[ - β( √(3) d / a - 1 ) ],

that only depends on the interbond distance d with V_ppπ = 2.8  eV and β = 3.37 <cit.>. This is a good approximation because the nearest neighbors lie approximately in the same plane such that the overlap integral is dominated by V_ppπ.

Our results show that optical displacements, i.e., the relative motion of different sublattices, effectively screen the PMF in twisted graphene, yielding much lower values as predicted by continuum elasticity <cit.>. Physically, we find that the relative displacements act to restore the microscopic 𝒞_3z symmetry by reducing changes in the bond length.


Corrugated graphene.—As a second example, we consider monolayer graphene subjected to a long-wavelength corrugation, see Fig. <ref>(c). This setup may be realized by engineering a suitable substrate <cit.> or through a buckling transition as was observed for graphene on NbSe_2 <cit.>. In particular, we consider a periodic corrugation with C_3v symmetry <cit.> commensurate with graphene and given by h_sub( r) = h_0 ∑_n=1^3 cos( _n · r + θ) where g_1 = 4π / ( √(3) L) ( 0, 1 ) and g_2,3 = 4π / ( √(3) L ) ( ∓√(3)/2, -1/2 ). These are the three shortest nonzero reciprocal vectors of the corrugation related by 𝒞_3z rotations. Here h_0 and θ control the amplitude and shape of the corrugation and L = Na is the superlattice constant with integer N≫1. 

The corrugated substrate is modeled in lammps with a dense honeycomb lattice to ensure commensuration and mimic a smooth substrate. Interactions in the graphene are modeled with REBO <cit.> and the interaction with the substrate by a 12-6 Lennard-Jones potential <cit.>. The resulting in-plane displacements are analyzed using a Helmholtz decomposition, u = ∑_ g( α_ gẑ× g + β_ g g ) e^i g · r / (ig^2) and similar for v. Here the transverse and longitudinal coefficients α_ g and β_ g are c-numbers corresponding to the curl and divergence, respectively. These coefficients are constrained by symmetry <cit.> and together with h_ g and w_ g there are four complex coefficients for each reciprocal star. For example, for θ = 0 (modulo π/3) the corrugation has C_6v symmetry. In this case, 𝒞_2z implies u( r) = - u(- r) and h( r) = h(- r) but v( r) =  v(- r) and w( r) = -w(- r) since a 𝒞_2z rotation exchanges the sublattices. 

As a first approximation, one can use continuum elasticity for the acoustic in-plane displacement field. In the limit where graphene is pinned to the substrate, we have h=h_sub <cit.> and the only nonzero in-plane coefficients are β_1 = (1 - 3ν) π^2 h_0^2 e^-2iθ / (3L^2), β_2 = (3-ν) π^2 h_0^2/ (3L^2), and β_3 = 2π^2 h_0^2  e^2iθ/ (3L^2) where the subscript indexes the star and ν is the Poisson ratio of graphene. From lammps, we also find that volumetric components are dominant for both u and v. While some rotational components are symmetry-allowed, they are at least one order of magnitude smaller <cit.>. Hence, unlike for twist moirés, u is mostly irrotational for corrugations. 

As was the case for TBG, optical displacements are significant even in the regime where continuum elasticity accurately describes the acoustic displacements (h_0/L < 0.02). While generally v is at least one order of magnitude smaller than u, the former contributes to the pseudo gauge field at zeroth order and the latter via the strain tensor, see Eq. (<ref>). Hence the acoustic part is suppressed by a factor a/L such that ultimately both fields contribute at the same order. In Fig. <ref>(a) we show the PMFs for a corrugation with h_0 = 5  Å and θ = 15^∘. We see that the total PMF is more concentrated and its magnitude is halved compared to B_ac. This is further illustrated in Fig. <ref>(b) where the RMS of the PMFs is plotted as a function of h_0/L. Note that for h_0/L < 0.02, the PMF from Eq. (<ref>) matches Eq. (<ref>) only if we include optical displacements. For large amplitudes h_0/L > 0.02, there are higher-order corrections which further reduce the magnitude of the PMF.


The suppression of the pseudo gauge field due to optical displacements strongly modifies the electronic minibands near charge neutrality. In Fig. <ref>, we show the bands calculated with the continuum model together with the bands from the tight-binding model in Eq. (<ref>) using atomic positions from lammps with the hopping amplitude from Eq. (<ref>). As expected, optical contributions reduce the minigaps between these bands and increase the bandwidth, while the topology of the bands remains unchanged. The latter is given by the valley Chern number 𝒞 = ( 𝒞_+ - 𝒞_- ) / 2 since the magnetic point group in a single valley breaks time-reversal symmetry. For a C_3v corrugation, it is given by 
C_3v(C_3) = < 𝒞_3z, ℳ_x 𝒯> where ℳ_x (x ↦ -x) is an in-plane mirror and 𝒯 is spinless time reversal. Importantly, the band manifold near charge neutrality is well reproduced by the continuum model only if we include optical contributions. However, the continuum model cannot reproduce the higher bands from tight binding. This discrepancy is likely due to higher-order corrections in the continuum theory such as a position-dependent Fermi velocity <cit.>. Moreover, using Eq. (<ref>) directly by taking its Fourier transform changes only slightly the continuum bands <cit.>.

The reduced band flattening and gaps from the elastically screened gauge field, has important consequences for the feasibility of symmetry-broken phases <cit.> and fractional Chern insulators <cit.> in periodically corrugated graphene. However, in the presence of an electric field normal to the xy plane, which couples to the height modulation <cit.>, one can still obtain isolated and flattened minibands. This results from the sublattice polarization induced by the PMF in real space <cit.> such that a scalar potential V( r) = V_0 h( r) / h_0 effectively acts as a staggered sublattice potential on the superlattice scale. This is shown in Fig. <ref>(c) where we plot the bandwidth W and gaps Δ of the highest valence band versus V_0 for h_0 = 6  Å and θ = 15^∘. As the gap opens at the Dirac point, the 2nd valence band closes and reopens to a topological phase with a minimum bandwidth of 20 meV for a field strength 40  mV/nm. Moreover, in this case the trace condition <cit.> is violated on the order of 10%, see Fig. <ref>(d). Hence, the reduced sublattice polarization from the screened PMF makes this system less favorable for hosting a fractional Chern insulator <cit.>.


Conclusions.—We developed a theory of pseudo gauge fields in graphene that takes into account contributions from both acoustic and optical displacement fields. The latter correspond to the relative motion of different sublattices. Using molecular dynamics simulations, we have shown that the optical displacements significantly modify the resulting pseudo magnetic field. Specifically, we applied our theory to moiré graphene and corrugated graphene. In all cases studied, optical contributions screen the acoustic contribution resulting in an overall reduction of the pseudo magnetic field. A simple explanation is given by that fact that the internal relaxation tends to restore the bond lengths to their pristine value.

Our theory elucidates the origin of discrepancies between continuum and tight-binding calculations that use microscopic theories to model lattice relaxation. It also introduces a novel way to engineer pseudo magnetic fields through the optical displacement though this may be difficult to achieve in practice and requires microscopic theories that go beyond continuum elasticity. Furthermore, our theory will help to understand experiments that probe pseudo gauge fields. For example, in twisted bilayer graphene, elastic screening of pseudo magnetic fields may explain why initial predictions based on continuum elasticity <cit.> were not observed.

We conclude that continuum elasticity, which only yields the acoustic displacement field, cannot fully describe pseudo gauge fields in graphene and most likely other low-dimensional materials with nonprimitive lattices, e.g., transition metal dichalcogenides. For displacements that are mostly in the nominal graphene plane, a reduction factor <cit.> may be used for qualitative results which reduces the strength of the pseudo magnetic field by a constant factor. However, the reduction factor is not universal and may vary depending on the amount of strain and microscopic details. Moreover, when out-of-plane displacements are significant, a reduction factor is inadequate even for qualitative results because both the magnitude and shape of the pseudo magnetic field are modified.


 

CDB and EJM are supported by the U.S. Department of Energy under Grant No. DE-FG02-84ER45118. RS, WNL and LC acknowledge support from Research Foundation-Flanders (FWO) research project No. G0A5921N. We thank Shaffique Adam and Mohammed M. Al Ezzi for interesting and helpful discussions. We further acknowledge Gayani N. Pallewela for sharing lammps data of twisted bilayer graphene for different molecular dynamics potentials.









Supplemental Material


 





















§ ELECTRONIC THEORY
 



 §.§ Continuum limit


The tight-binding Hamiltonian of graphene in the nearest-neighbor approximation in the presence of strain can be written as

    H = -∑_ r∑_n=1^3 t_n( r) c_A^†( r) c_B[  r + δ_n( r) ] + h.c.,

where the first sum runs over cells r and the second one over nearest neighbors. Here t_n( r) > 0 is the hopping amplitude between nearest neighbors that is modulated by strain, and c_σ^†( r) [c_σ( r)] are creation (annihilation) operators for sublattice σ = A,B. The position of A atoms is given by r +  u_A( r) + h_A( r) ẑ where u_A is the in-plane displacement, i.e, projected on the original graphene xy plane, and h_A is the out-of-plane displacement in the z direction, and similar for B atoms. 

To study the low-energy physics near the K_± point, we take the continuum limit. This amounts to the replacement

    c_σ( r) →√(A_c)∑_τψ_τσ( r) e^i  K_τ· r,

where A_c = A/N is the unit cell area. Here, ψ_τσ^†( r) [ψ_τσ( r)] are field operators that create (annihilate) a fermion of sublattice σ at position r composed of small momentum components | q|a ≪ 1 near valley K_τ with τ = ± 1 the valley index, and which obeys the usual fermionic anticommutation relations. The effective Hamiltonian becomes


    H_eff    = -∑_τ∑_n=1^3 ∫ d^2  r   t_n( r) e^i  K_τ·δ_n( r)ψ_τ A^†( r) ψ_τ B[  r + δ_n( r) ]  + h.c.
       ≈ -∑_τ∫ d^2  r  ψ_τ A^†( r) ∑_n=1^3 e^i  K_τ·δ_n^0[ t_0 δ_n^0 ·∇_ r + δ t_n( r) ] ψ_τ B( r) + h.c.,


where we let ∑_ r→ A_c^-1∫ d^2  r. In the second line, we expanded everything up to lowest order in gradients and displacements with t_n( r) = t_0 + δ t_n( r), where we take t_0 = 2.8 eV <cit.>. We further defined the nearest-neighbor bond vectors of pristine graphene δ_n^0. For example, for the orientation shown in Fig. <ref>, we have δ_1^0 = ( 0, a_0 ), δ_2^0 = a_0 ( -√(3)/2, -1/2 ), and δ_3^0 = a_0 ( √(3)/2, -1/2 ) where a_0 = 1.42  Å is the nearest-neighbor distance.


Moreover, we have assumed that intervalley coupling is negligible. This is justified in the limit L ≫ a where the strain field varies slowly with respect to the graphene lattice. In the following, we consider a general orientation of the graphene where φ is the angle between the zigzag direction and the x axis. For example, we have φ = 0 for zigzag orientation, shown in Fig. <ref>, and armchair orientation would correspond to φ = -π/2. We now take K_± = ± R(φ) (4π/3a, 0) with R(φ) the standard 2 × 2 rotation matrix for a counterclockwise rotation by an angle φ. One finds <cit.>

    -t_0 ∑_n=1^3 e^i  K_τ·δ_n^0δ_n^0 ·∇_ r = -i ħ v_F  e^iτφ( τ∂_x - i ∂_y ),

with ħ v_F = √(3) t_0 a/2 and we define the pseudo vector potential A = (A_x, A_y) as

    -∑_n=1^3 δ t_n( r) e^i  K_τ·δ_n^0≡ e v_F e^iτφ[ A_x( r) - i τ A_y( r) ],

with -e the electron charge. Explicitly,

    A( r) = R(φ)/2ev_F[ δ t_2 + δ t_3 - 2 δ t_1;   √(3)( δ t_3 - δ t_2 ) ].

The effective low-energy Hamiltonian thus becomes


    H_eff = ħ v_F ∑_τ∫ d^2  r  ψ_τ^†( r) {[ -i ∇_ r + eτ/ħ A( r) ] · (τσ_x, σ_y) }ψ_τ( r),

with ψ_τ( r) = [ e^-iτφ/2ψ_τ A( r), e^iτφ/2ψ_τ B( r) ]^t. 

To compute the change in hopping amplitude δ t_n due to strain, we first consider the change in the nearest-neighbor bond vectors. In lowest order of displacements and their gradients, we have

    δ_n( r)- δ_n^0     =  u_B( r + δ_n^0) -  u_A( r) + [ h_B( r + δ_n^0) - h_A( r) ] ẑ
       ≈( δ_n^0 ·∇) [  u( r) + h( r) ẑ] - [  v( r) + w( r) ẑ],


where

    u     = ( u_A + u_B ) / 2,            h = ( h_A + h_B ) / 2, 

    v     = u_A - u_B,            w = h_A - h_B,

are the center-of-mass (acoustic) and relative (optical) displacements, respectively. For example, in a classical microscopic theory of the in-plane phonon modes of a pristine graphene sheet, one can show that in the long-wavelength limit <cit.>

    v( r) = ( κ - 1 ) a/2√(3)  R(3φ) [ ∂_x u_y + ∂_y u_x; ∂_x u_x - ∂_y u_y ],

with κ∼ 1/3 the so-called reduction factor, whose precise value depends on microscopic details. We do not attempt to find a relation between the optical and acoustic displacements for the cases we consider, e.g., moiré graphene and corrugated graphene. Instead, our microscopic theory is given by a molecular dynamics simulation which yields the fields u_A,B and h_A,B. 

Next we calculate the change in the hopping amplitude. To this end, we assume that t( d) = t(d) only depends on the bond distance d (two-center approximation) and expand them in lowest order of the displacements:


    δ t_n     = t( δ_n ) - t_0 
       ≈. ∂ t/∂ d_i|_0 ( δ_n - δ_n^0 )_i + 1/2. ∂^2 t/∂ d_z^2|_0 ( δ_n - δ_n^0 )_z ( δ_n - δ_n^0 )_z 
        = -3β t_0/a^2[ δ_n^0 ·( δ_n - δ_n^0 ) + 1/2( δ_n - δ_n^0 )_z ( δ_n - δ_n^0 )_z ], 
       ≡ -3β t_0/a^2[ δ_ni^0 δ_nj^0 u_ij - δ_ni^0( v_i + w ∂_i h ) + w^2/2],

where β = -a/√(3)t_0. ∂ t/∂ d|_nn,0≈ 3 <cit.>, u_ij is the strain tensor, and we used that d_z = 0 in the absence of strain. We further used (i = x, y)

    ∂ t/∂ d_i    = d_i/d∂ t/∂ d, 
    ∂ t/∂ d_z    = d_z/d∂ t/∂ d, 
    ∂^2 t/∂ d_z^2    = 1/d∂ t/∂ d - d_z^2/d^3∂ t/∂ d + d_z^2/d^2∂ t^2/∂ d^2, 
    ∂^2 t/∂ d_z ∂ d_i    = d_i d_z/d^2∂^2 t/∂ d^2 - d_i d_z/d^3∂ t/∂ d.


Plugging the result for δ t_n into Eq. (<ref>), we obtain the pseudo vector potential

    A =  A_ac +  A_op = √(3)ħβ/2ea[ R(3φ) [ u_yy - u_xx; u_xy + u_yx ] + 2 √(3)/aẑ×(  v + w ∇ h ) ],

which is preserved under a global 𝒞_3z rotation, as both the graphene lattice and the valley are left unchanged. Note that we omitted a leading-order contribution proportional to

    ∑_n=1^3 e^i  K_τ·δ_n ( r)   ≈ i  K_τ·∑_n=1^3 [ δ_n ( r) - δ_n^0 ] e^i  K_τ·δ_n^0
       ≈ i  K_τ·∑_n=1^3 e^i  K_τ·δ_n^0[ ( δ_n^0 ·∇)  u( r) -  v( r) ] 
        = -2π e^iτφ/√(3)( ∂_x u_x - i τ∂_y u_x ),

since it corresponds to a gauge transformation which does not affect the long-wavelength physics <cit.>.


For completeness, we also consider the second nearest-neighbor hopping t_m'( r) between atoms of the same sublattice. This yields an additional term V( r) ψ_τ^†( r) ψ_τ( r) in the effective Hamiltonian of Eq. (<ref>) where V( r) is the deformation potential. In lowest order of the displacements and up to a constant energy shift, we find 

    V( r)     = -∑_m=1^6 δ t_m'( r) e^i  K_τ·η_m^0
        = δ t_1' + δ t_2' + δ t_3' 
        = 3a/2. ∂ t/∂ d|_nnn,0( u_xx + u_yy),

where η_1^0 =  a_1, η_2^0 =  a_2, η_3^0 = - a_1 -  a_2, η_4^0 = - a_1, η_5^0 = - a_2, and η_6^0 =  a_1 +  a_2 are the six second nearest-neighbor bond vectors of pristine graphene with a_1/2 = δ_1^0 - δ_2/3^0 the primitive lattice vectors. Here we also used the two-center approximation such that δ t' is independent on reversal of the bond vector. We can estimate the prefactor in front of tr(u) = u_xx + u_yy in Eq. (<ref>) with a simple model for the hopping amplitude,

    t(d) = t_0 exp[ - β( d / a_0 - 1 ) ],

which yields a prefactor -3 √(3)β t(a) / 2. The deformation potential has full rotational symmetry and is finite even in the presence of microscopic 𝒞_3z such as for biaxial strain <cit.>. Moreover, it does not depend on the relative displacements in lowest order because it couples equal sublattices.



 §.§ Tight binding
 

In this section, we describe the tight-binding method used to calculate the electronic properties. The advantage of this approach is that it allows us to construct a model directly from the atomic configurations obtained from the molecular dynamics simulations, fully taking into account the microscopic details of the lattice relaxation. The tight-binding Hamiltonian is given by

    H = -∑_< i, j > t_ij c_i^† c_j + ∑_i ξ_i c_i^† c_i,

where the atomic sites are labeled by i and j and the first sum runs over nearest neighbors. The nearest-neighbor hopping amplitude t_ij is approximated as

    t_ij = t_0 exp[- β( | r_i -  r_j|/a_0 - 1 ) ],

where we use the values t_0 = 2.8 eV, a_0 = a / √(3) = 0.142 nm, and β = 3.37 <cit.>. This is expected to be a good approximation since the nearest-neighbors lie approximately in the same plane, which generally differs from the xy plane in the presence of corrugation. Hence, the relevant overlap integral is still given by V_ppπ = t_0.

The second term in Eq. (<ref>) is an on-site electrostatic potential, which can originate from an external electric field or the pseudo scalar field <cit.>. 



 §.§ Pseudo magnetic field on a discrete grid
 

We can further calculate the pseudo magnetic field (PMF) directly from the atomic positions, allowing us to determine the PMF that effectively enters the tight-binding calculations. For the zigzag orientation shown in Fig. <ref>, the vector potential is given by

    A_x( r) - i A_y( r) = -1/e v_F∑_n=1^3 δ t_n( r) e^i K ·δ_n( r),

where δ t_n = t_ij - t_0 is the change in hopping energy, K_+ = 4 π / (3a) x̂ the zone corner of the graphene BZ and 

    δ_n( r) = δ_n^0     +  u_B( r + δ_n^0) + h_B( r + δ_n^0) ẑ
        -  u_A( r) - h_A( r) ẑ,

the modified bond vector. Please note that the definition in Eq. (<ref>) differs from our previous definition since we also take into account changes in the phase factor. These only modify the PMF at next-to-leading order as the lowest order contribution from the phase factor corresponds to a gauge transformation. We prefer this definition for this section since all contributions to (<ref>) are automatically taken into account in tight-binding calculations, whereas in the continuum theory we only consider leading-order terms from an expansion in displacements and momentum <cit.>.

The PMF is given by B  = ∇× A and calculated on a discrete (strained) atomic grid using Stokes' theorem. In two dimensions, this relates the surface integral of the curl of a vector field A to the contour integral around the boundary of the same field:

    ∬_Σ (∇× A) · dΣ = ∮_∂Σ A · d  r.

If we consider B to be slowly varying within the graphene unit cell, the surface integral in (<ref>) becomes trivial:

    ∇× A( r)    ≈1/S∮_∂Σ A · d  r 
        = 1/S( ∫_ r + δ_1^ r + δ_2 + ∫_ r + δ_2^ r + δ_3 + ∫_ r + δ_3^ r + δ_1)  A · d  r 
       ≈1/S[  A( r + δ_1) +  A( r + δ_2) ] ·δ_2 - δ_1/2
        + 1/S[  A( r + δ_2) +  A( r + δ_3) ] ·δ_3 - δ_2/2
        + 1/S[  A( r + δ_3) +  A( r + δ_1) ] ·δ_1 - δ_3/2
        =  A( r + δ_1) ·δ_2 - δ_3/2S +  A( r + δ_2) ·δ_3 - δ_1/2S
        +  A( r + δ_3) ·δ_1 - δ_2/2S,

where we approximated the line integral along the triangle contour shown in Fig. <ref>. The enclosed area S is calculated by taking the cross product of any two sides,

    S = |(δ_1 - δ_2) × (δ_3 - δ_2)|/2.

Finally, we note that this formula is only gauge invariant up to the same order of approximation. In particular, if we send A ↦ A + ∇χ and approximate ∇χ( r+δ_n) ≃∇χ( r) + ( δ_n ·∇) ∇χ( r) one can show that the final expression in Eq. (<ref>) remains unchanged.



§ MOIRÉ GRAPHENE
 

In this section, we discuss the displacement fields from lattice relaxation in twisted bilayer graphene (TBG) obtained from lammps simulations. In particular, we consider commensurate structures that have the periodicity of the moiré lattice with twist angles defined by <cit.>

    cosϑ_m = 3m^2 + 3m + 1/2/3m^2 + 3m + 1.

Moreover, we place the twist center at the center of a graphene hexagon such that ϑ = 0 corresponds to AA stacking. These structures have point group D_6 = < 𝒞_6z, 𝒞_2x> where 𝒞_6z is a rotation by π/3 about the z axis and 𝒞_2x is a π rotation about the x axis <cit.>, as illustrated in Fig. <ref>(a).

We then define the displacement fields as

    r_lσ =  r_lσ^0 +  u_lσ( r_lσ^0) + h_lσ( r_lσ^0) ẑ,

where l= 1,2 is the layer, σ = A, B the sublattice, and 
r_lσ^0= R[(-1)^l+1ϑ/2] ρ_σ are the rigid coordinates in the absence of relaxation. Here ρ_σ = n_1  a_1 + n_2  a_2 + δ_σ are the atomic positions of monolayer graphene with n_1,n_2 ∈ Z and δ_σ the sublattice position in the graphene cell.



 §.§ Displacement fields


Assuming the moiré periodicity is preserved after lattice relaxation, we define the smooth fields

    u_lσ( r) = ∑_ g u_lσ, g e^i  g · r,

and similar for out-of-plane displacements. Here g are moiré reciprocal vectors and u_ g =  u_- g^* are complex Fourier components. In practice, the Fourier components are obtained by taking a discrete Fourier transform of the lammps data. We can now define the acoustic and optical displacement fields for each layer,

    u_l     = ( u_lA + u_lB ) / 2,            h_l = ( h_lA + h_lB ) / 2, 

    v_l     = u_lA - u_lB,            w_l = h_lA - h_lB.

One can now make similar superpositions between layers in terms of homo and hetero displacements. For example, one distinguishes between out-of-plane buckling (homo) and breathing (hetero) displacements. The hetero displacements are given by

    u     = u_1 - u_2,            h = h_1 - h_2, 

    v     = v_1 - v_2,            w = w_1 - w_2,

which are shown in Fig. <ref> for ϑ≈ 1.018^∘. We see that u is mostly solenoidal, i.e., ∇· u ≈ 0. It gives rise to local co-twisting near the AA stacking center (origin) and counter-twist near AB and BA stacking centers <cit.>. This reduces the size of AA regions and increases the size of AB and BA regions. Similarly, the acoustic out-of-plane hetero displacements, i.e., the interlayer distance, conforms to the in-plane stacking. The in-plane optical displacement field v has significant volumetric contributions and is over one order of magnitude smaller than u, while the out-of-plane field w is about six orders of magnitude smaller than h and can be safely neglected.


The in-plane components can be written using a Helmholtz decomposition. For example,

    u_ g = α_ gẑ× g + β_ g g/ig^2,

for g = | g| ≠ 0 and where α_ g = α_- g^* and β_ g = β_- g^* are complex numbers. Note that u_ 0 is a constant relative shift between layers which does not affect the long-wavelength physics for small twists. These coefficients are related to the divergence and curl:

    ∇× u     = ∑_ g i g × u_ g e^i  g · r = ẑ∑_ gα_ g e^i g · r, 
    ∇· u     = ∑_ g i g · u_ g e^i  g · r = ∑_ gβ_ g e^i g · r,

which are the rotational and in-plane volumetric components of the displacement gradient. From the lammps simulations, we find that u is dominated by real rotational coefficients while v is mostly given in terms of imaginary volumetric coefficients. In Fig. <ref>, we show these Fourier coefficients as function of twist angle for the first six reciprocal stars. The scaling of the first Fourier component of the acoustic displacement field yields an estimate of the V_1 / μ where V_1 is the first Fourier coefficient of the adhesion energy <cit.> and μ is the shear Lamé constant. From this, we estimate that relaxation in our molecular dynamics simulations is about 1.5 stronger as compared to density-functional theory calculations using the local-stacking approximation <cit.>. We further find that the Re(  g × u_ g) [Im(  g · u_ g)] can be fitted to a polynomial odd (even) in 1/θ.


All results obtained from lammps molecular dynamics simulations are consistent with the emergent D_6 symmetry of twisted bilayer graphene. To illustrate how symmetries constrain the displacement fields, consider an in-plane symmetry 𝒮. The displacement fields then satisfy 𝒮 u( r) =  u(𝒮 r) and h( r) = h(𝒮 r). In reciprocal space, one then finds that α_ g and β_ g transform as a pseudoscalar and scalar, respectively. Explicitly,

    h_𝒮 g    = h_ g, 
    β_𝒮 g    = β_ g, 
    α_𝒮 g    = (𝒮) α_ g.

The optical displacements transform similarly except for the fact that any transformation that interchanges the sublattices gives an extra minus sign. For example, in the presence of 𝒞_2z rotation symmetry we have u(- r) = - u( r) and h(- r) = h( r), while v(- r) =  v( r) and w(- r) = -w( r). The symmetry-allowed coefficients for the first five reciprocal stars are listed in Table <ref>. Here each reciprocal star consists of six reciprocal vectors closed under 𝒞_6z rotations.




 §.§ Pseudo magnetic fields


The valley-preserving symmetries of the emergent moiré lattice in small-angle twisted bilayer graphene form the dichromatic group 6'2'2 = < 𝒞_6z𝒯, 𝒞_2y> also denoted as D_6(D_3) = D_3 + (D_6 ∖ D_3 ) 𝒯 <cit.> where 𝒯 is spinless time-reversal symmetry with 𝒯^2 = 1. These symmetries yield the following constraints on the pseudo gauge fields,

    A_1,2( r)     =  A_1,2(- r) + ∇χ, 
     A_1,2( r)     = 𝒞_3z^-1 A_1,2(𝒞_3z r) + ∇χ, 
     A_1(x,y)     = [ -1  0;  0  1 ] A_2(-x,y),

where the subscript is the layer index and χ is a scalar function. This implies that the pseudomagnetic field (PMF) satisfies

    B_1,2( r)     = -B_1,2(- r) = B_1,2(𝒞_3z r), 
    
        B_1(x,y)     = -B_2(-x,y) = B_2(x,-y),

such that we only need to consider one layer. To verify these symmetries, we show the PMFs for both layers obtained from the acoustic and optical displacement fields in Fig. <ref> for two twist angles.


Similar as for the displacement fields, we can write the PMF as a Fourier series

    B( r) = ∑_ g ≠ 0 B_ g e^i g· r,

where we restrict the sum since the uniform part vanishes for periodic displacement fields, giving a zero net pseudoflux. We show the Fourier components of the acoustic (B_ac) and total PMF (B_tot = B_ac + B_op) for the first three reciprocal stars in Fig. <ref> as a function of twist angle. We see that the PMF is well approximated by

    B_1,2(x,y) ≃± 2 B_0 ∑_i=1^3 sin( g_i · r),

where g_1,2,3 are the shortest nonzero moiré reciprocal vectors related by 120^∘ rotations. Here |B_0| ≈ 3  T for twist angles in the range 1.5^∘ < ϑ < 4^∘ and where we used β = 3.37. The independence of the magnitude of the PMF for twist angles above the magic angle (ϑ≈ 1^∘) was also found in other studies that only considered the acoustic part <cit.>. Note that the form in Eq. (<ref>) satisfies the symmetry constraints in Eq. (<ref>) since the first star is symmetric under x ↦ -x in the presence of 𝒞_3z symmetry. 


Finally, the root mean square is then given by

    RMS = √(< B^2 >) = √(∑_ g |B_ g|^2),

which is shown in Fig. <ref> of the main text. If we define the magnetic length from the RMS, one finds ℓ_RMS > 9  nm for the twist angles under consideration. We also show the ratio √(<B_ac^2> / < B_tot^2 >) in Fig. <ref> as a function of twist angle. This ratio gives an estimate of the reduction factor <cit.> due to optical displacements.




 §.§ Molecular dynamics
 

For the molecular dynamics simulation of twisted bilayer graphene we divide the interactions between interlayer and intralayer. For the interlayer interaction we take the dihedral-angle-corrected registry-dependent potential (DRIP) benchmarked with EXX-RPA DFT calculations <cit.>, which has proven to be an improvement in comparison to the usual Kolgomorov-Crespi <cit.> potential. For the intralayer interactions we use the usual REBO potential <cit.> with a 2 Å cutoff to avoid interactions between atoms in different layers. For the geometric optimization we enforce periodic boundary conditions and use the "fire" minimization style which uses damped dynamics. This ensures that the system does not get stuck in a local energy minimum. As a stopping criteria, we take the force tolerance to be 10^-4 eV/Å between subsequent steps. 




 §.§ Comparing different microscopic models


In Fig. <ref> we show the first-star transverse Fourier component of the pseudo vector potential of the first layer, ẑ·(  A_1, g_1× g_1 ), and the RMS of the PMF as a function of twist angle for different MD potentials. In the main text we showed results for the interlayer DRIP <cit.> and intralayer AIREBO <cit.> potential. Here we also show results using the Kolgomorov-Crespi (KC) <cit.> interlayer potential instead of DRIP. Moreover, we also show data <cit.> for the intralayer REBO <cit.> and the same interlayer DRIP as before. We finally also show data for DRIP with the intralayer machine-learning potential GAP20 <cit.>.

We see that the specific microscopic details, corresponding here to different MD potentials, not only yield different total PMFs, but also different acoustic and optical contributions. Most strikingly, the results for DRIP with GAP20 yield a noticeably smaller reduction as compared to the other cases which have similar reduction factors. Even though DRIP with GAP20 gives smaller acoustic contributions, i.e., less strain compared to using REBO or AIREBO for intralayer interactions, the optical contributions are much smaller such that in the end the total PMF is larger for GAP20.



§ CORRUGATED GRAPHENE
 

In this section, we consider the structural and electronic properties of monolayer graphene subjected to a periodic corrugation. Specifically, we consider a triangular height modulation with C_3v symmetry, commensurate with the graphene lattice and defined by <cit.>

    h_sub( r) = h_0 ∑_n=1^3 cos( _n · r + θ),

with amplitude h_0 and where θ controls the shape. The superlattice is defined by the reciprocal vectors

    _1 = 4π/√(3)L[ 0; 1 ],     _2,3 = 4π/√(3)L[ ∓√(3)/2;    -1/2 ],

and _3 = - _1 - _2 where L = Na is the lattice constant of the height modulation. The corresponding lattice vectors l_1,2 are chosen such that g_i · l_j = 2 πδ_ij. Here we have taken the coordinate system shown in Fig. <ref> with the zigzag direction along the x axis.

The height profile from Eq. (<ref>) preserves C_3v = < 𝒞_3z, ℳ_x > symmetry on the superlattice scale where 𝒞_3z is a rotation by 120^∘ about the z axis and ℳ_x is a mirror x ↦ -x. Note that these are not the microscopic symmetries of the graphene lattice, which are broken by the corrugation. In general, the corrugation breaks 𝒞_2z rotation symmetry since this operation is equivalent to θ↦ -θ. Moreover, 
in the long-wavelength limit L ≫ a we can restrict to θ∈ [0, π/3[. This follows from the fact that ±θ are 𝒞_2z partners, while θ↦θ + 2π/3 is equivalent to a translation y ↦ y + 2L/√(3). Hence, for the special case θ = 0   mod  π/3, the point group of the superlattice becomes C_6v = < 𝒞_6z, ℳ_x >.



 §.§ Symmetry constraints


To minimize the elastic energy, the corrugated graphene lattice will relax, giving rise to in-plane displacement fields. Since the corrugation is smooth and periodic on the atomic scale, the in-plane and out-of-plane displacement fields can be written in terms of Fourier series,

    u_σ( r)     = ∑_ u_σ, e^i· r, 
    
         h_σ( r)     = ∑_ h_σ, e^i· r,

respectively, where σ = A, B is the sublattice index,  are reciprocal lattice vectors of the corrugation, and h_ = h_-^* and u_ =  u_-^* are Fourier components. Here the uniform components (| g|=0) are set to zero as these only result in an overall translation of the graphene.


As described in the main text, we can then define acoustic and optical displacement fields. The displacement fields calculated with lammps for a corrugated substrate with h_0 = 5  Å and θ = 15^∘ are shown in Fig. <ref>. Since the relaxed structure is assumed to be adiabatically connected to the rigid corrugation, the displacements fields obey the same symmetries. In particular, due to 𝒞_3z symmetry and the reality of the fields, the displacements are characterized by three complex numbers α_m, β_m, and h_m for each star of reciprocal vectors, which is indexed by m. Here, we define a reciprocal star by six reciprocal vectors that are related by 𝒞_6z and we use a Helmholtz decomposition for the in-plane fields. 

We show the symmetry-allowed values for a corrugation with C_3v symmetry in Table <ref> for the first five stars. 
We also show the volumetric Fourier coefficients of the in-plane displacements, i.e., projected on the xy plane, in Fig. <ref> for θ = 15^∘ as a function of h_0/L for the first nine reciprocal stars. 





 §.§ Continuum elasticity


The long-wavelength acoustic displacements in graphene can be modeled using the continuum theory of elasticity <cit.>. Here one views the graphene as an elastically isotropic membrane. The elastic potential energy <cit.>
and substrate interaction are modeled as <cit.>

    H_elas    = 1/2∫ d^2  r [ λ u_ii u_ii + 2 μ u_ij u_ji + κ( ∇^2 h )^2 ], 
    
        H_sub    = γ/2∫ d^2  r [ h( r) - h_sub( r) ]^2,

where λ and μ are in-plane Lamé constants, κ is the out-of-plane bending rigidity, and γ controls the interaction with the substrate. Here summation over repeated indices is implied and

    u_ij( r) = 1/2( ∂ u_j/∂ r_i + ∂ u_i/∂ r_j + ∂ h/∂ r_i∂ h/∂ r_j),

is the strain tensor with i,j = x,y. We further assume that the graphene is pinned to the substrate such that h( r) = h_sub( r). This is justified in the limit L ≫( κ / γ)^1/4≈ 1  nm <cit.> where L is periodicity of the corrugation and the numerical value is for graphene on SiO_2 <cit.>. In this case, the interaction with the substrate dominates over the curvature term [last term of Eq. (<ref>)] since the latter scales as κ h^2 / L^4 while the former scales as γ h^2. In this limit, H_elas is a functional of the in-plane field only.

Under these assumptions, the elastic energy density


    ℋ_elas = 1/A_c∫_cell d^2  r [ ( λ/2 + μ) ( u_xx^2 + u_yy^2 ) + λ u_xx u_yy + 2 μ u_xy^2 ] + constant,


where A_c is the area of the supercell defined by the periodic height modulation. For the triangular height profile, we have A_c = √(3) L^2 / 2. For convenience, we define the tensor

    f_ij( r) ≡[ ∂_i h( r) ] [ ∂_j h( r) ] = ∑_ f_ij  e^i · r,

where

    f_ij = - ∑_' h_' h_ - ' g_i' ( g_j - g_j'),

with f_ij,- = f_ij^*. The strain tensor thus becomes

    u_ij( r) = 1/2∑_[ i ( g_i u_j + g_j u_i) + f_ij] e^i · r,

where we set u_i 0 = 0 since this amounts to a uniform translation. In terms of the Helmoltz decomposition, 

    u_xx g    = g_x ( g_x β_ g - g_y α_ g)/g^2 + f_xx g, 
    
        u_yy g    = g_y ( g_y β_ g + g_x α_ g)/g^2 + f_yy g, 
    
        u_xy g    = 2 g_x g_y β_ g + ( g_x^2 - g_y^2 ) α_ g/g^2 + f_xy g,

with g = | g| and where α_ g and β_ g are the rotational and volumetric components of the in-plane displacement field, respectively. Plugging the Fourier expansions into the energy density ℋ_elas gives <cit.>


    1/A_c∫_cell d^2  r   u_ii^2     = 1/A_c∑_,'∫ d^2  r ( i g_i u_i + f_ii /2) ( i g_i' u_i' + f_ii '/2) e^i (  + ' ) · r
        = ∑_| i g_i u_i + f_ii /2|^2, 
    1/A_c∫_cell d^2  r   u_xx u_yy    = ∑_( i g_x u_x + f_xx /2) ( -i g_y u_y^* + f_yy ^*/2) 
        = 1/2∑_[ ( i g_x u_x + f_xx /2) ( -i g_y u_y^* + f_yy ^*/2) + c.c.], 
    1/A_c∫_cell d^2  r   u_(xy)^2     = 1/4∑_( i g_x u_y + i g_y u_x + f_xy) ( -i g_x u_y^* - i g_y u_x^* + f_xy^* ).

We obtain

    ℋ_elas    = ( λ/2 + μ) ∑_( i g_x u_x + f_xx /2) ( -i g_x u_x^* + f_xx ^*/2) 
        + ( λ/2 + μ) ∑_( i g_y u_y + f_yy /2) ( -i g_y u_y^* + f_yy ^*/2) 
        + λ/2∑_[ ( i g_x u_x + f_xx/2) ( -i g_y u_y^* + f_yy^*/2) + c.c.] 
        + μ/2∑_( i g_x u_y + i g_y u_x + f_xy) ( -i g_x u_y^* - i g_y u_x^* + f_xy^* ) + κ/2∑_ | h_ |^2 g^4.

Minimizing with respect to u_i^* for nonzero g yields the solutions for the Fourier components u_i in terms of f_ij (and thus h_i). We find

    ∂ℋ_elas/∂ u_x^*    = -i g_x [ ( λ + 2μ) ( i g_x u_x + f_xx /2) + λ( i g_y u_y + f_yy /2) ] - i μ g_y ( i g_x u_y + i g_y u_x + f_xy), 
    ∂ℋ_elas/∂ u_y^*    = -i  g_y [ ( λ + 2 μ) ( i g_y u_y + f_yy /2) + λ( i g_x u_x + f_xx /2) ] - i μ g_x ( i g_x u_y + i g_y u_x + f_xy),

and setting these equations equal to zero, gives

    u_x    = i/2 ( λ + 2 μ) g^4{ f_xx g_x [ g_x^2 ( λ + 2 μ) + g_y^2 ( 3 λ + 4 μ) ] + ( f_yy g_x - 2 f_xy g_y ) [ g_x^2 λ - g_y^2 ( λ + 2 μ) ]  }, 
    
        u_y    = i/2 ( λ + 2 μ) g^4{ f_yy g_y [ g_y^2 ( λ + 2 μ) + g_x^2 ( 3 λ + 4 μ) ] + ( f_xx g_y - 2 f_xy g_x ) [ g_y^2 λ - g_x^2 ( λ + 2 μ) ]  },

and

    α_ g    = i  g × u_ g = ( f_xx g - f_yy g) g_x g_y + f_xy g( g_y^2 - g_x^2 )/g^2 = 1/g^2∑_ g' h_ g' h_ g- g'(  g' · g ) ( ẑ· g' × g ), 
    β_ g    = i  g · u_ g = μℱ_ g/λ + 2 μ - f_xx g + f_yy g/2,

with ℱ_ g = ( g_x^2 f_yy - 2 g_x g_y f_xy + g_y^2 f_xx) / g^2. We find that α_ g = 0 when h( r) is restricted to the first star. For completeness, we also give the explicit forms for f_ij g for the height profile from Eq. (<ref>):

    f_ g_1    = 2π^2/3h_0^2 e^-2iθ/L^2[  3  0;  0 -1 ], 
     f_ g_1+2 g 2    = 2π^2/3h_0^2/L^2[ -3  0;  0  1 ], 
     f_2 g_1    = 2π^2/3h_0^2 e^2iθ/L^2[  0  0;  0 -2 ],

which transform as f_𝒮 g = 𝒮 f_ g𝒮^-1.


Using the relations

    μ = E/2 ( 1 + ν ),     λ = ν E/1 - ν^2,

for isotropic linear elastic two-dimensional materials, with E the Young modulus and ν the Poisson ratio of graphene, we obtain <cit.>

    u_xx + u_yy    = 1-ν/2ℱ_ g, 
    
    
    
        u_xx - u_yy    = 1+ν/2g_y^2 - g_x^2/g^2ℱ_ g, 
    
        u_xy + u_yx    = -1+ν/22 g_x g_y/g^2ℱ_ g,

for nonzero g and u_ij 0 = f_ij 0/2. These are the volumetric and shear strains, respectively. In this work we use the value ν = 0.165 from experiment <cit.>. Importantly, since |ℱ_| ∼ h^2 / L^2 we expect the linear theory to be valid only for L ≫ h.

For the triangular corrugation with C_3v symmetry, continuum elasticity gives a relaxed in-plane acoustic displacement field that is irrotational such that α_ g = 0 even though it is symmetry allowed in some stars. The volumetric components are finite only in the first three stars with

    β_1 = β_ g_1    = ( 1 - 3 ν) π^2/3h_0^2/L^2  e^-2iθ, 
    β_2 = β_ g_1+2 g_2    = ( 3 - ν) π^2/3h_0^2/L^2, 
    β_3 = β_2 g_1    = 2π^2/3h_0^2/L^2  e^2iθ,

consistent with the symmetry analysis. Hence, continuum elasticity yields u = ∇ϕ with ϕ_ g = -β_ g/g^2.



 §.§ Electronic band structure


In the presence of a periodic scalar and gauge fields, the valley-projected electronic continuum Hamiltonian from Eq. (<ref>) can be diagonalized by Fourier transformation,

    ψ_τ( r) = 1/√(A)∑_ k∑_ e^i ( k -  ) · r  c_τ, k - ,

where A is the system size and the sum over k is restricted to the superlattice Brillouin zone (SBZ) and c_τ, k-^† (c_τ, k-) are two-component creation (annihilation) operators that create a fermion in valley τ with momentum k -. The Hamiltonian becomes


    H = 1/A∑_τ∑_ k,  k'∑_, '∫ d^2 r   c_τ, k'-'^† e^-i( k' - ') · r{ħ v_F [  k -  + τ e/ħ  A( r) ] ·( τσ_x, σ_y ) + V( r) σ_0 } e^i( k - ) · r c_τ, k-.


Next, for any function f( r) with the periodicity of the superlattice, we have

    ∫ d^2  r   e^-i( k' - ') · r  f( r)   e^i( k - ) · r = A δ_ k k' f_-',

where we used that k -  g with k in the SBZ and g a reciprocal vector of the superlattice, is a unique momentum decomposition. We further used ∫ d^2  r = ∑_ R∫_cell d^2 r with ∑_ R e^i  k · R = N δ_ k 0 where the sum runs over superlattice cells. We then obtain


    H = ∑_τ∑_ k∑_, ' c_τ, k-'^†{ħ v_F [ (  k - ) δ_' + τ e/ħ  A_-'] ·( τσ_x, σ_y ) + V_-'σ_0 } c_τ, k-,
  

which can be diagonalized numerically by taking a sufficient number of  vectors for convergence. All results shown in this work were obtained with a cutoff || < 12k_0 where k_0 = 4π / 3L.

For the triangular periodic height profile that is aligned with the zigzag direction, the electronic continuum theory has wallpaper group 14 (p3m1) with point group C_3v and superlattice translations. The symmetries of the valley-projected Hamiltonian are generated by 𝒞_3z and ℳ_x 𝒯 where ℳ_x is a mirror symmetry (x ↦ -x) and 𝒯 is spinless time-reversal symmetry with 𝒯^2 = 1. While ℳ_x and 𝒯 both interchange valleys, their combination leaves the valleys invariant. This yields the magnetic point group 3m' also denoted as C_3v(C_3) <cit.>. 

In Fig. <ref>(a) we show the band structure calculated with the continuum model in the presence of a perpendicular electric field V( r) = V_0 h( r)/h_0 where h( r) is the relaxed height profile. Parameters are h_0 = 6  Å, θ = 15^∘, and V_0 = 71  meV corresponding to Fig. <ref>(c) and (d). We also show the Berry curvature and quantum metric of the first valence band in valley K_+.




 §.§ Quantum geometry


We calculate the quantum geometry of an isolated Bloch band in a given valley with the gauge-invariant product. To this end, consider a square plaquette in the Brillouin zone of area δ^2 centered at k with corners: k_1 =  k + δ2(-1, -1), k_2 =  k + δ2(-1, 1), k_3 =  k + δ2(1, 1), and k_4 =  k + δ2(1, -1). The gauge-invariant product is then given by 
 
    ∏_m=1^4 ⟨ u_ k_m | u_ k_m+1⟩,

where k_5 =  k_1. It is straightfoward to show that

    tr( g_ k)     = lim_δ→ 0δ^-2 Re( 1 - ∏_m=1^4 ⟨ u_ k_m | u_ k_m+1⟩), 
    Ω_ k    = lim_δ→ 0δ^-2 ∏_m=1^4 ⟨ u_ k_m | u_ k_m+1⟩,

where

    g_ k^μν    = Re( ⟨∂^μ u_ k | ∂^ν u_ k⟩) + ⟨ u_ k | ∂^μ u_ k⟩⟨ u_ k | ∂^ν u_ k⟩, 
    Ω_ k    = -2  Im( ⟨∂_k_x u_ k | ∂_k_y u_ k⟩),

is the (single band) Fubini-Study quantum metric with tr( g_ k) = g_ k^xx + g_ k^yy and Berry curvature, respectively. They form the real and imaginary components of the (single band) quantum geometric tensor,

    Q_ k^μν    = ⟨∂^μ u_ k | ( 1 - | u_ k⟩⟨ u_ k | ) | ∂^ν u_ k⟩
        = g_ k^μν - i/2 ϵ^μνΩ_ k.

One can further show that <cit.>

    tr(g_ k) ≥ | Ω_ k |,

which is the trace inequality. The trace condition refers to the saturation of this bound and holds for Landau levels. One route of engineering flat bands that may host fractional Chern insulators at fractional filling of the band, is so-called Landau level mimicry <cit.>. For example, a flat Bloch band with unit Chern number that satisfies the trace condition emulates the lowest Landau level, whose exact ground state at one over odd integer filling in the presence of short-range interactions is given by a Laughlin-type wave function <cit.>.



 §.§ Valley polarization


The band structure from tight-binding calculations always features Kramers' pairs of minibands that approximately belong to different valleys since the long-wavelength corrugation does not couple the valleys. We quantify this from the expectation value of the valley polarization P_v. The simplest choice for P_v is given in terms of the next-nearest neighbor Haldane hopping,

    P_v = 1/i3√(3)∑_<< i,j > > e^3iϕ_ij c_i^† c_j,

where ϕ_ij is the angle of the next-nearest neighbor bond vector with the x axis. For pristine graphene, one can diagonalize P_v in momentum space:

    P_v = ∑_ k g( k) c_ k^† c_ k,

with g( K_± +  q) = ± 1 + 𝒪(q^2a^2). In the presence of potentials that vary slowly with respect to the graphene lattice constant a, the valleys remain approximately decoupled such that < P_v > ≈± 1 for valley K_±.



 §.§ Molecular dynamics


For the lammps simulations, we consider a graphene sheet spanned by l_1,2 with L = | l_1,2 | = 59a ≈ 14.5 nm. In order to model a generic substrate that induces a smooth van der Waals force on the graphene sheet, we use a honeycomb lattice with lattice constant a/3. This ensures commensurability with the graphene and an interaction potential that is smooth on the graphene lattice scale due to the higher density of the substrate. We then apply the out-of-plane displacement field given by Eq. <ref> to the substrate for h_0 ranging from 0.05 nm to 0.6 nm in steps of 0.05 nm and for θ = 15^∘. Energy minimization was performed with the lammps code <cit.> using the FIRE minimization procedure, using force tolerance of 10^-6 eV/Å as a stop criteria. To describe the potential energy landscape of the system we use a combination of the AIREBO interatomic potential <cit.> for short-range interactions between carbon atoms in the graphene sheet with a cutoff of 2 Å, and a 12-6 Lennard-Jones potential for the interaction between the graphene and the substrate. The Lennard-Jones interaction was parametrized with an energy constant of 10 meV, a zero-crossing distance of 3.7 Å, and a cut-off radius of 10 Å. 

