# Fileset

[pnas.2307151120.pdf](https://mdr.nims.go.jp/filesets/823cd107-84fd-4ab7-bbd0-a8d1daad78a6/download)

## Creator

Xiaoyu Wang, Joe Finney, Aaron L. Sharpe, Linsey K. Rodenbach, Connie L. Hsueh, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), M. A. Kastner, Oskar Vafek, David Goldhaber-Gordon

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Unusual magnetotransport in twisted bilayer graphene from strain-induced open Fermi surfaces](https://mdr.nims.go.jp/datasets/a379ab4d-c607-44bf-9bfe-7d69e5b8885b)

## Fulltext

Unusual magnetotransport in twisted bilayer graphene from strain-induced open Fermi surfacesRESEARCH ARTICLE PHYSICS OPEN ACCESSUnusual magnetotransport in twisted bilayer graphene fromstrain-induced open Fermi surfacesXiaoyu Wanga,1 ID , Joe Finneyb,c,1 ID , Aaron L. Sharped , Linsey K. Rodenbachb,c , Connie L. Hsuehc,e , Kenji Watanabef ID , Takashi Taniguchig ,M. A. Kastnerb,c,h,2 ID , Oskar Vafeka,i,2 , and David Goldhaber-Gordonb,c,2Contributed by Marc Kastner; received May 8, 2023; accepted July 15, 2023; reviewed by Leon Balents and Emanuel TutucAnisotropic hopping in a toy Hofstadter model was recently invoked to explain a richand surprising Landau spectrum measured in twisted bilayer graphene away from themagic angle. Suspecting that such anisotropy could arise from unintended uniaxialstrain, we extend the Bistritzer–MacDonald model to include uniaxial heterostrainand present a detailed analysis of its impact on band structure and magnetotransport.We find that such strain strongly influences band structure, shifting the threeotherwise-degenerate van Hove points to different energies. Coupled to a Boltzmannmagnetotransport calculation, this reproduces previously unexplained nonsaturatingB2 magnetoresistance over broad ranges of density near filling � = ±2 and predictssubtler features that had not been noticed in the experimental data. In contrast to thesedistinctive signatures in longitudinal resistivity, the Hall coefficient is barely influencedby strain, to the extent that it still shows a single sign change on each side of the chargeneutrality point—surprisingly, this sign change no longer occurs at a van Hove point.The theory also predicts a marked rotation of the electrical transport principal axes asa function of filling even for fixed strain and for rigid bands. More careful examinationof interaction-induced nematic order versus strain effects in twisted bilayer graphenecould thus be in order.moiré materials | heterostrain | magnetotransport | bilayer grapheneThe discovery of superconductivity and correlated insulating states in magic-angle twistedbilayer graphene (TBG) (1, 2) placed the material at the forefront of condensed matterphysics research (3–18). The moiré superlattice potential of TBG, resulting from a smallrelative twist angle � between the graphene layers, can induce nearly flat, topologicallynontrivial, isolated bands, consisting of electronic states near the Dirac points of eachmonolayer of graphene (19). As a result, TBG is an exceptional platform for studyingthe interplay of electron correlations and band topology (20–38).Strain—especially differing lattice distortions in the two layers, termed heterostrain—isbelieved to play an important role in the phase diagram of TBG (37–41). Scanning probemeasurements typically find uniaxial heterostrain in the range of 0.1 to 0.7% (3, 9, 11, 42)in samples fabricated with the tear-and-stack method (43, 44). For heterostrain, asopposed to strain applied equally to both layers, the linear distortion of the moiré unitcell is amplified by a factor of ∼ 1/� relative to the linear distortion of the microscopicatomic lattice. For example, 0.2% uniaxial heterostrain causes an ∼ 8% change in thelargest linear dimension of the moiré unit cell for a twist angle of 1.38◦. The moiré unitcell area changes by much less, only∼ 0.1%, but because we infer twist angle from moiréunit cell area in transport, this dependence of area on strain still leads to underestimatesof the uncertainty in twist angles presented in the transport literature, as noted in ref. 3.In a recent report by some of the authors [Finney et al. (45)], a TBG sample with amoiré unit cell area of 90 nm2 (corresponding to � = 1.38◦) displayed several unusualphenomena in magnetotransport. As anticipated for a twist well above the magic angle,the sample did not exhibit the strong interaction-driven effects typically observed in near-magic-angle devices. Surprisingly, though, over a broad filling range near half-filling thelongitudinal magnetoresistivity (MR) exhibited a B2 increase up to ∼ 100-fold at ≈ 5T, after which quantum oscillations set in. The authors found that a toy Hofstadtermodel with anisotropy showed multiple features similar to those in the data, and theyconjectured that uniaxial strain might cause such anisotropy.In this work, we present a systematic theoretical study of the impact of uniaxialheterostrain on the narrow-band dispersion of TBG above the magic angle, analyzeits consequences for weak field magnetotransport, and compare it with experimentaldata from ref. 45. We base our theory on the Bistritzer–MacDonald (BM) continuumSignificanceBecause of its rich array ofcorrelated phases, twisted bilayergraphene (TBG) near the magicangle has captivated thecondensed matter physics world.The large moiré length scale notonly promotes interaction-relatedeffects but also allows forextrinsic factors such as strain toplay a major role. In a previouswork, we presentedmeasurements of a TBG devicewith several unusual behaviorsin magnetotransport andconjectured that uniaxial straincould explain our measurements.Here, we modelmagnetotransport in TBGby incorporating uniaxialheterostrain into theBistritzer–MacDonaldHamiltonian. The theory not onlyreproduces the unusualphenomena from the previouswork but also predicts additionalfeatures unnoticed before. Ourwork therefore demonstratesthe crucial role of heterostrain inTBG devices.Competing interest statement: K.W. and T.T. are onmultiple papers with both reviewers due to their positionas materials sources in this field. Because they providehexagonal boron nitride to so many workers in the field,it is almost impossible to find referees who are notcoauthors with them.Copyright © 2023 the Author(s). Published by PNAS.This open access article is distributed under CreativeCommons Attribution License 4.0 (CC BY).1X.W. and J.F. contributed equally to this work.2To whom correspondence may be addressed.Email: mkastner@mit.edu, vafek@magnet.fsu.edu, orgoldhaber-gordon@stanford.edu.This article contains supporting information onlineat https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2307151120/-/DCSupplemental.Published August 14, 2023.PNAS 2023 Vol. 120 No. 34 e2307151120 https://doi.org/10.1073/pnas.2307151120 1 of 10Downloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.http://crossmark.crossref.org/dialog/?doi=10.1073/pnas.2307151120&domain=pdf&date_stamp=2023-08-10https://orcid.org/0000-0002-5035-5761https://orcid.org/0000-0001-7166-9754https://orcid.org/0000-0003-3701-8119https://orcid.org/0000-0001-7641-5438https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/mailto:mkastner@mit.edumailto:vafek@magnet.fsu.edumailto:goldhaber-gordon@stanford.eduhttps://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2307151120/-/DCSupplementalhttps://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2307151120/-/DCSupplementalmodel (19), incorporating heterostrain in the form of a deforma-tion potential, a pseudomagnetic field (46, 47), and a distortionof the moiré pattern.Our key theoretical result is that heterostrain lifts the energeticdegeneracy of the two Dirac points as well as that of the threevan Hove points of a given band. The splitting of the two Diracpoints leads to a semimetallic state with small Fermi pocketsnear the charge neutrality point (CNP). More interestingly, thesplitting of the van Hove points leads to open Fermi surfaces(FSs) in the filling range bounded by two of the van Hove points.In the weak field semiclassical regime governed by the Boltzmannequation, the open FSs generally lead to a nonsaturating B2 MR(48), accounting for this previously unexplained feature in theexperimental data of ref. 45.This theory makes a number of falsifiable predictions. 1) If thedirection of current flow in the lithographically patterned Hallbar is not aligned with the 1D-like principal axis of transportin the open FS regime, longitudinal magnetoresistance shouldmix substantially into the measured resistance at Hall contacts.2) A subtle cusp should appear in resistivity as filling crosses thelowest-energy van Hove point. 3) A Lifshitz transition from twoFS pockets to one should also coincide with crossing this lowervan Hove point. We reanalyze experimental data from ref. 45and find that these predictions are verified. The strained BMmodel studied here has electron–hole symmetry. We leave thediscussion of electron–hole asymmetry in the experimental datato future works.The theory also has implications that are not so far directlyprobed by the experiment, but are striking. First, on each sideof the CNP, the divergence and sign change in the Hall numbernear half-filling of a 4-fold degenerate band does not coincidewith any of the van Hove points but instead occurs within thefilling range where FSs are open. Second, the transport principalaxis continuously rotates by up to 90◦ as density is tuned from theCNP to the open FS regime. Such rotation of the transport axeshas been presented as evidence for interaction-induced nematicorder (14), but here, we find that it can arise purely due tostrain-induced band structure effects.This work clearly demonstrates that the effects of evenminiscule amounts of heterostrain in TBG cannot be neglected.Dramatic and unexpected phenomena occur in strained TBGeven in the single-particle regime, without the strong correlationeffects that arise near the magic angle. Given the amplifying effectof a small heterostrain on the moiré length scale, it is tantalizingto consider strain engineering of such devices to achieve effectsthat would be impossible in regular solids due to structuralinstabilities.The paper is organized as follows: In Section I, we presentthe theoretical calculation of the changes in band structureresulting from uniaxial strain. In Section II, we use the Boltzmanntransport equation to calculate the electron transport propertiesin magnetic fields resulting from the strain-induced alterations inthe band structure. Section III is a comparison of our predictionswith experiment, and it is sufficiently self-contained that a readerwho does not need all the theoretical details can jump directlythere. In Section IV, we summarize our results.1. Geometric and Energetic Effects of UniaxialHeterostrain on TBGIn the limit of small deformations, both uniaxial heterostrain anda small twist angle are captured via a coordinate transformation:r′l = r + ul (r), where l = t, b labels the Top (Bottom) graphenelayers, and ul (r) ≈ El r is the local deformation field. Thesymmetric and antisymmetric part of the 2×2 tensor El describesstrain and rotation, respectively. For twist angle (�) and a uniaxialheterostrain of strength (�) and direction ('), we parameterizeEt = −Eb ≡ E/2, where E ≡ T (�) + S(�,'), and given by:T (�) =(0 −�� 0), S(�,') = RT'(−� 00 ��)R'. [1]Here, R' is the two-dimensional rotation matrix, and � ≈0.16 is the Poisson ratio (3). Physically, � > 0 correspondsto compressing the Top layer while stretching the Bottom layeralong the direction determined by ', as illustrated in Fig. 1A.A relative deformation E between the graphene bilayers gener-ates a moiré superlattice, with moiré reciprocal lattice vectorsgi=1,2 = ETGi=1,2, where Gi are reciprocal lattice vectors ofthe undeformed monolayer graphene. The moiré lattice vectorsLi=1,2 are uniquely defined through the relation Li · gj = 2��ij.It is important to note that only relative deformations generatethe moiré superlattice. Deformations that act homogeneously onthe two layers have little effect on the narrow band physics, andwe accordingly neglect them in this work*.Under rotation R', the strain tensor transforms as a headlessvector that remains invariant under '→ '+ 180◦. Combinedwith the C3z symmetry of the undeformed graphene lattice, thestrained electronic dispersion within a given graphene valleysimply rotates 60◦ under ' → ' + 60◦. We therefore onlyreport results for ' ∈ [0◦, 60◦). For concreteness, we define themicroscopic unit cell vectors ai=1,2 of an undeformed graphenelattice as a1 = a( 12 ,−√32 ), a2 = a(1, 0), where a ≈ 2.46Å is the lattice constant. The positions of the sublattice A, Bwithin a unit cell are chosen as E�A = (0, 0) and E�B = a√3(0, 1).The reciprocal lattice vectors are G1 = 4�√3a(0,−1) and G2 =4�√3a(√32 , 12 ). Different conventions lead to different definitionsof the Dirac Hamiltonian (see for instance ref. 40), but the physicsis consistent.Fig. 1 A and B illustrates the geometric effects of heterostrainfor twist angle � = 1.38◦. At � = 0, AA-stacked regions of themoiré form an equilateral triangular superlattice. Introducinguniaxial heterostrain changes the spacings between neighboringAA-stacked regions (Li=1,2,3). For � = 0.2%, typical in thesesystems (3, 9, 11), the variation in spacings can be as large as�/� ≈ 8%. Such dramatic amplification of the microscopic strainmakes moiré materials uniquely suited to strain engineering—conventional materials become structurally unstable at distor-tions only 10% as large as those achieved in the moiré superlattice.Note that the effect of uniaxial heterostrain on the moiré unit cellarea is small at �2�2/�2, as some spacings become larger whileothers become smaller (SI Appendix, section I).We proceed to discuss the energetic effects in the context ofthe continuum BM model (19). We work in the limit where bothEl and the wavevector k in the moiré Brillouin zone are smalland consider only the leading order terms in both. This wouldmean, for instance, that terms such as Ek are omitted as higher-order terms. This treatment is generally justified away from themagic angle because higher-order terms can play an importantrole only close to the magic angle where the bandwidth becomescomparable to these terms (49, 50). We explicitly checked that*We checked numerically that adding a small homogeneous strain in addition to aheterostrain of similar strength yields band and transport properties almost identicalto those produced by adding a heterostrain alone.2 of 10 https://doi.org/10.1073/pnas.2307151120 pnas.orgDownloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialsϵ = 0.2 %ϵ = 0.2 % , φ = 0∘ ϵ = 0.2 % , φ = 40∘ ϵ = 0.2 % , φ = 60∘ϵ = 0.2 %L1L2L3xyϵ = 0.2 %A B G HC D E FFig. 1. (A) Schematics of applying a uniaxial heterostrain on the pair of microscopic unit cells of monolayer graphene making up TBG. (Upper sketch) Orange(blue) color corresponds to the Top (Bottom) layer. The uniaxial strain of strength +(−)�/2 and direction ' on the Top (Bottom) layer are represented as coloredarrows. (Lower sketch) Deformation of the moiré superlattice for twist angle 1.38◦ due to a uniaxial heterostrain of � = 0.2% and ' = 0◦. Unstrained (gray,dashed) and strained (black, solid) triangular lattice sites of AA stacking regions of the moiré superlattice are depicted. (B) Dependence of the three moirétriangular bond lengths on ' for a fixed strength. (C–F ) Energy maps of the Upper band of the BM Hamiltonian in valley K, plotted in the moiré Brillouin zonespecified by k = k1g1 + k2g2, where k1,2 ∈ [0,1). There are six special points of the band structure, i.e., two Dirac points (black stars), three van Hove points(colored dots), and one band maximum (black cross). The contour lines intersecting the van Hove points are plotted and labeled by their respective fillingfractions. In the unstrained case (C), the two Dirac points and three van Hove points are respectively at equal energies. The energy degeneracies are lifted inthe presence of uniaxial heterostrain, as illustrated in (D–F ). This leads to semimetallic behavior at the CNP, and a '-dependent filling range near � = 2 withopen FSs. (G and H) '-dependence of the energies and filling fractions of the band structure special points for a fixed heterostrain strength. The backgroundcolormap is the calculated density of states, with a broadening of � = 1meV. Green (blue) color represents high (low) density of states. The energetic minimumand maximum of the narrow bands are shown with horizontal dashed gray lines.at � ≈ 1.38◦, the effects of such higher-order terms are indeednegligibly small. To leading order, the strained BM Hamiltonianfor a given valley is given byH� = (∑l=t,bH intra�,l ) + H inter� , [2]where � = ±1 labels K (K′) valleys of monolayer graphene. Theinterlayer Hamiltonian is given byH inter� ≈∫d2r †�,t(r) ∑j=1,2,3T�,je−i�qj ·r �,b(r) + h.c.,[3]where  �,l (r) ≡ ( �,l,A(r), �,l,B(r))T is a spinor in thesublattice basis for a given valley and layer. We have suppressedthe spin index for simplicity. qj=1,2,3 are the three nearestneighbor bonds of the reciprocal honeycomb lattice, andT�,j = w0�0 + w1(cos2�(j − 1)3�x + � sin2�(j − 1)3�y).[4](�0, �x , �y) are Pauli matrices acting on sublattice degrees offreedom.The intralayer Hamiltonian is given by:H intra�,l = �∫d2r †�,l (r)(tr[El ]�0) �,l (r)−ℏvFa∫d2r †�,l (r)[(−i∇ − A�,l ) · (��x , �y)] �,l (r).[5]Here, the first term is the deformation potential that couplesto the electron density. Its value is not precisely known inthe literature, with numbers ranging from −4.1 eV to 30 eVdepending on the methodology (51–54). We use � = −4.1 eV inthis work based on first principles calculations (54), although forheterostrain � ≈ 0.2% varying the deformation potential over therange proposed in the literature leads to only minor quantitativedifferences in band dispersions. A�,l is the pseudovector potentialthat comes from changes in the intersublattice hopping due todeformations and changes sign between graphene valleys. It isgiven as refs. 46 and 47: A�,l =√3�2a �(�l,xx − �l,yy,−2�l,xy),where we choose � ≈ 3.14 from refs. 3 and 40. We further fixℏvF/a = 2.68eV based on Fermi velocity in monolayer graphenevF ≈ 106m/s (55), w0 = 88meV, and w1 = 110meV (19) inour calculations and also set ℏ = 1 in the remainder of the paper.To leading order approximation, the strained BM Hamilto-nian in a given valley (Eq. 2) has particle–hole symmetry underP l (r) =∑l ′ i(�y)ll ′ l ′(−r) (56), where �y is a Pauli matrixacting on the layer degrees of freedom. This means that forevery single-electron state at energy E and wavevector k, thereis a state at energy −E and wavevector −k. This particle–holesymmetry has been investigated extensively for the unstrainedBM model, e.g., refs. 26 and 57, and here, it is generalized to thestrained case. The source of experimentally evident particle–holeasymmetry in the off-magic-angle device (45) could be higher-order gradient terms beyond what is captured in the BM modelin Eq. 2, interaction effects (58–61), or their combination.Fig. 1 D–F displays the effects on the band structure of � =0.2% heterostrain applied in select directions relative to the x-axisdefined in Fig. 1A, as specified by ' ∈ [0◦, 60◦). For simplicity,we show only contour maps of the upper band from valley K inthe moiré Brillouin zone specified by k = k1g1 + k2g2, wherek1,2 ∈ [0, 1). Heterostrain preserves C2T and valley U (1) (23),so the Lower and Upper bands remain connected via two Diracpoints. The Upper band features six special points—two DiracPNAS 2023 Vol. 120 No. 34 e2307151120 https://doi.org/10.1073/pnas.2307151120 3 of 10Downloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.points (black stars), three van Hove points (colored dots), andone band maximum (black cross). The six special points of agiven band are related to “critical points” in the context of theMorse theory, which states that∑i(−1)
i = � , [6]where 
i is the index of the i-th critical point, and � is theEuler characteristic of a manifold (62); � vanishes for theBrillouin zone which is a torus. Although a Dirac point isstrictly a point of nonanalyticity and is not directly covered byMorse theory, if we imagine adding a tiny gap term, it willbecome a legitimate band extremum allowing Morse theoryto apply. Whereas the two band minima (Dirac points) andthe band maximum have even 
 , and so each contributes +1to the sum, every conventional van Hove point (i.e., not ahigher order) has an odd 
 and contributes −1. The overallsum thus vanishes. Therefore, the van Hove points can onlybe annihilated/created by colliding with local minima/maxima.For a relatively small heterostrain as shown in Fig. 1, thenumber of special points per band is the same as at � = 0.However, for larger heterostrain (e.g., � = 0.5%, see SI Appendix,Fig. S1), more striking behavior of the special points can occur,such as a change in their total number via aforementionedcollisions and the appearance of tilted type II Dirac cones(63, 64).A key finding of the present work is that the respective energydegeneracies of the two Dirac points and the three van Hovepoints are lifted by uniaxial heterostrain, by amounts dependingsensitively on '. In the absence of strain, Fig. 1C, the three vanHove points are at equal energy, and separate closed contours ofconstant energy centered around the Dirac points from closedcontours centered around the band maximum. As illustrated inFig. 1 D–F, uniaxial heterostrain splits the energy degeneracyof the two Dirac points, leading to a semimetallic state withsmall Fermi pockets near the CNP (40). The three van Hovepoints also split in energy. The two outermost van Hove points(i.e., closer to the band maximum) bound a filling range of openFSs near � = 2, whereas the innermost van Hove point movescloser to one of the Dirac points. If we continue increasing �,a collision of the critical points occurs, the innermost van Hovedisappears, the two Dirac points become type II tilted, and anew ordinary minimum is created. Note that a small mass addedto type II tilted Dirac points will not introduce band extrema,and as a consequence, type II tilted Dirac points are not criticalpoints of Morse theory. Therefore, after the collision, Eq. 6still holds.Interestingly, the elongation of the FSs shows a strong fillingdependence. Close to the CNP, the bigger Fermi pocket thatencloses a Dirac point is stretched along a direction perpendicularto that of the open FSs; see Figs. 1 D–F. As explained later,this leads to a dramatic rotation of the principal transportaxis when the filling is tuned from the CNP to the openFS range.The dependence of the energy and filling of the band structurespecial points on ' at a fixed � is shown in Fig. 1 G and H. Ofnotable interest is the sensitivity of the filling range with open FSsto '. This filling range must in fact vanish at some ' between0◦ and 60◦, when the energies of the two outermost van Hovepoints cross. As seen in Fig. 1 D–F, this also alters the elongationof the open FSs.2. Boltzmann Equation and Magnetoresistivityin TBGHaving understood the heterostrain effects on the bandstructure,we proceed to discuss the implications for magnetotransport. Webegin by considering the general structure of the two-dimensionalresistivity tensor � subject to heterostrain. The resistivity tensoris defined via (ExEy)=(�xx �xy�yx �yy)(jxjy), [7]where E = (Ex , Ey)T and j = (jx , jy)T are electric field andcurrent vectors, respectively. Under rotation by ��, the resistivitytensor transforms as:�′ = RT���R�� , R�� =(cos �� − sin ��sin �� cos ��). [8]If the underlying system has a point group symmetryhigher than C2z (e.g., C3z , C6z), then � = �0I − i�H �y isthe most general form of � invariant under such rotations.Here, �y is the Pauli matrix acting in the two-dimensionalcoordinate basis, �0(−B) = �0(B) is the longitudinal resistivity,and �H (−B) =−�H (B) is the Hall resistivity. The even/oddparity under time reversal is guaranteed by the Onsager reciprocalrelations.Since heterostrain breaks the point group symmetry down toC2z , we generally expect �xx 6= �yy, �xy 6= −�yx . Nevertheless, itis always possible to define transport principal axes after a suitablerotation �� of the coordinate system, such that�principal =12(�1 + �2)I +12(�1 − �2)�z + �H i�y. [9]Here, �1,2 are longitudinal resistivities along the principaltransport directions ê1,2, respectively. The rotation angle �� isdetermined up to 180◦ by requiring �1 < �2.Below, we first derive the MR tensor using a Boltzmannapproach for a general noninteracting electronic system withinthe relaxation time approximation. Since there is currentlylimited understanding of the scattering mechanisms determiningelectrical transport in TBG, here, we follow ref. 65 and use therelaxation time approximation. We then present the resultsfor heterostrained TBG, showing that in the open FS region,the low-resistivity principal axis (ê1) is nearly perfectly alignedwith the shortest moiré bond direction. However, there is adramatic rotation of the principal axis as the filling moves towardthe CNP. We further show that the open FSs lead to a B2nonsaturating MR along ê2 and a saturating resistivity along ê1.For random orientation (�0) of the principal axis to the electricalcurrent axis in the Hall bar geometry, e.g., as in ref. 45, thelongitudinal resistivity is given by: �xx = �1 cos2 �0 + �2 sin2 �0.This is dominated by the �2 ∼ B2 component. As a result, theexperimental measurements should observe the nonsaturatingMR component if there is a misalignment of the Hall barorientation with respect to the principal transport axis. Suchmisalignment is generically to be expected and indeed mustoccur over most of the relevant filling range since the principalaxes rotate with filling while the device geometry remainsfixed.Boltzmann Equation and Method of Characteristics. We beginwith a brief description of the method of characteristics used to4 of 10 https://doi.org/10.1073/pnas.2307151120 pnas.orgDownloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialssolve the Boltzmann equation perturbatively in electric field Ebut without a restriction on the strength of the perpendicularmagnetic field B = Bẑ, as long as the semiclassical regime holds(66). Due to C2zT symmetry of TBG at B = 0, there is no Berrycurvature contribution to the semiclassical equations of motion.Then, within the relaxation time approximation, the Boltzmannequation for a given energy band becomes∂nk∂t+ (qE + qvk × B) ·∂nk∂k= −nk − n0,k�, [10]where qE+qvk×B is the total force on the Bloch electrons, withvk ≡ ∇k"k and charge q; n0,k is the equilibrium Fermi–Diracdistribution and nk is the desired nonequilibrium distributionfunction.We consider a stationary solution to the Boltzmann equationby parameterizing the distribution function as:nk = n0,k + n1,k. [11]As a result, the Boltzmann equation for the deviation of thedistribution function from equilibrium is:(qE · vk)∂n0,k∂"k+ (qvk × B) ·∂n1,k∂k= −n1,k�. [12]Note that the magnetic field only couples to n1 since (qvk×B) ·∇kn0,k = (qvk × B) · vk∂"kn0,k = 0.To solve the above partial differential equation (PDE), we seeka family of curves covering the k-space which we parameterizeas k(s) with s ∈ [0, s0), such that along these curves, the PDEbecomes an ordinary differential equation (ODE). If a curve k(s)satisfiesdk(s)ds= qv(s)× B, [13]then n1,k(s) ≡ n1(s) satisfies(qE · vk)∂n0,k∂"k|k=k(s) +dn1(s)ds= −n1(s)�. [14]Becaused"(s)ds= v(s) ·dk(s)ds= 0, [15]the curve k(s) must coincide with the contour of constant energy.Thus, the Boltzmann equation becomes:[qE · v(s)]∂n0(s)∂"(s)+dn1(s)ds= −n1(s)�. [16]The ODE is readily solved withn1(s) = �0e−s/� − e−s/�∫ s0ds′es′/� [qE · v(s′)]∂n0(s′)∂"(s′). [17]where �0 is a constant determined by the following argument.Since k(s) describes a constant energy contour in a two-dimensional Brillouin zone, it is either a closed contour or severalopen contours that terminate on boundaries of the Brillouin zonesuch that they form a closed loop on a torus. In either case, k(s)is periodic under s → s + s0 modulo a moiré reciprocal latticevector, where s0 is the periodicity. The periodicity conditionn1(s0) = n1(0) leads to�0 =11− es0/�∫ s00ds′es′/�(qE · v(s′))∂n0(s′)∂"(s′), [18]which determines the desired n1(s).In the low-temperature limit, the steady-state current from agiven energy band is calculated as:j� = q∫d2k(2�)2 v�kn1,k =q2B(2�)2∫d"∫ s00dsv�(s)n1(s)=q3B(2�)�!c∞∑n=−∞v�n v�−n1 + in!c�E� ,[19]where (�, �) = x, y, and we have defined the cyclotron frequencyas:!c ≡ 2�/s0. [20]We have also made use of the periodicity of velocity unders → s + s0 to write it in terms of Fourier series, v(s) =∑∞n=−∞ vne−in!c s.To show that the second line of Eq. 19 holds, note that atevery k, we can define a local coordinate system (êv, ês) suchthat v ≡ vêv where v ≥ 0, and ês = êv × ẑ. The infinitesimalwavevector can be equivalently written as:dk = dkx êx + dky êy = dks ês + dkv êv.Eq. 13 can then be written as dk/ds = qvBês, or equivalentlydks = qvBds. As a result,∫dkxdky =∫dksdkv = qB∫d"ds.The conductivity tensor is therefore given by the followingexpression:��� =q3B2��!c∞∑n=−∞v(�)n v(�)−n1 + in!c�. [21]Eq. 21 gives the magnetoconductivity for a given FS contour.In the case of multiple FS contours and multiple bands—associated for example with spin and valley degeneracy in TBG—conductivities from different FS contours and bands add. Finally,the MR tensor is obtained by inverting the conductivity tensor,i.e., � =(∑n,i �n,i)−1, where n, i are band and contour labelsrespectively for a given energy level.To better understand Eq. 21, consider an example of aparabolic dispersion with "k = 12m0(k2x + k2y ), where m0 is thebare electron mass. At a fixed energy �, the contour is a circleparameterized as (kx , ky) =√2m0�(cos �, sin �), � ∈ [0, 2�).Using the method of characteristics, we get d�ds = − qBm0, or� = �0 − !0s, where !0 ≡qBm0is the cyclotron frequencyof bare electrons. This leads to the periodicity in s to bes0 = 2�/!0, where we have chosen the clockwise trajectory suchthat s0 > 0. The Fourier series of the velocity along the constantenergy contour is given by vx(s) =√�2m(e−i!0s + ei!0s), andvy(s) =√�2m1i(e−i!0s − ei!0s). Substituting into Eq. 21, weobtain the conductivity tensor:� = q2��2�11 + !20�2(1 −!0�!0� 1). [22]PNAS 2023 Vol. 120 No. 34 e2307151120 https://doi.org/10.1073/pnas.2307151120 5 of 10Downloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.Note that the total number density of filled electrons is givenby n =∫ d2k(2�)2 Θ(� − "k) = m0�2� . We therefore reproduce thewell-known magnetoconductivity tensor:� =nq2�m011 + !20�2(1 −!0�!0� 1). [23]In this simple example of a closed FS, the longitudinalresistivity is given by m0nq2� , independent of the magnetic field.The average of the velocity field, vn=0 ≡1s0∫ s00 dsv(s), vanishes.However, for an open FS, generally, vn=0 6= 0, i.e., electronshave a finite drift velocity when a magnetic field causes them totraverse the contour (SI Appendix, Fig. S2). The impact of such afinite drift velocity on the magnetotransport can be qualitativelyunderstood using the following example: In the expression forthe conductivity tensor (Eq. 21), we consider vxn=0 6= 0 butvyn=0 = 0. This corresponds to an open FS with a drift velocityalong the x direction. In the high-field limit ( !c� ∝ B � 1),we truncate the Fourier series at the leading order, yielding�open FS ≈q3B2��!c (vx0)2−2Im(vx−1vy1)!c�2Im(vx−1vy1)!c�|vy1|2!2c �2 , [24]where we made use of the equality: v−n = v∗n . Inverting thematrix, we obtain the MR tensor:�open FS ≈(2�)!cq3B�14Im(vx−1vy1)2 + (vx0)2|vy1|2×(|vy1|2 2Im(vx−1vy1)!c�−2Im(vx−1vy1)!c� (v(x)0 )2 (!c�)2).[25]We see that for an open FS, the longitudinal MR hasnonsaturating B2 behavior along the axis with a zero drift velocity(ŷ in the above example), and saturating behavior (constant in Bin this simple model) along the other axis.3. Comparison between Theory andExperimentAs we will see, the theory satisfactorily explains the weak-fieldmagnetotransport measurements presented in ref. 45. We thenpresent two predictions of the theory that we did not anticipateprior to starting this work: the dependence of the principal axisof transport on filling, and the behavior of magnetoresistanceand quantum oscillations at densities between the CNP and theonset of quadratic MR. The former has important implications,but it cannot be confirmed with our present datasets because oflimitations of the Hall bar geometry. The latter can be consideredsmoking gun evidence for the presence of the lowest-energy vanHove point and the energetic splitting of the Dirac cones.We present calculations for � = 1.38◦ (independentlymeasured for the region of the experimental device we focus on),and � = 0.2% and ' = 0◦, parameters which are not measuredin the experiment but are chosen to yield reasonable quantitativeagreement between the theoretical and experimental results withrespect to the filling range of open FSs and to the frequenciesof magnetoresistance oscillations to be presented later. Thegeneral phenomena of open FSs and quadratic MR hold fora broad range of heterostrain parameters � and '. We do notperform fine-tuning of these input parameters for two reasons: 1.We do not expect our strained BM model in Eq. 2 to yieldprecise quantitative agreement with experiment. Specifically,the model has particle–hole symmetry, which is absent fromexperimental measurements. More sophisticated noninteractingmodel calculations (49, 50) as well as interaction renormalizations(61) would likely be necessary to properly account for suchdetails. 2. Increasing � broadens the filling range that displaysopen FSs, but for a given �, varying ' strongly tunes this range(Fig. 1 G and H ), so there is some flexibility in assigning the twoparameters to match the experimental measurements. This mightbe overcome by also seeking to match the position of the lowestvan Hove singularity and the evolution of sizes of the two Fermipockets near the CNP, but the simplifications in the BM modelnoted in (1) above caution us against drawing strong quantitativeconclusions from the comparison.In Fig. 2, we show the computed MR along the principaltransport axes (A) and the Hall number (B). For comparison,we plot the experimentally measured longitudinal and transverseresistivities (C ) and Hall number (D) for the TBG device studiedin ref. 45.In the filling ranges with open FSs, the calculated �2(B)exhibits quadratic nonsaturating MR, whereas �1(B) saturates.The filling range for which quadratic MR occurs is bounded bythe two outermost van Hove points of the zero-field strainedband structure. In experiment, we observe quadratic MR inlongitudinal resistivity within a similar range of fillings. Morestrikingly, we observe quadratic MR in the transverse resistivityas well. In some cases, with increasing field, the symmetricpart of the transverse resistivity becomes larger than that of thelongitudinal resistivity. As discussed earlier, this degree of mixingcan be attributed to the misalignment between the strain-inducedbut filling-dependent principal axis of transport and the directionof current flow in the Hall bar geometry.At the first van Hove point (� ≈ ±0.6), the nonanalyticityin the density of states leads to a cusp in the first derivativeof the zero-field resistivity with respect to filling (SI Appendix,Fig. S6). As shown in Fig. 2A, atB 6= 0 the longitudinal resistanceas a function of filling develops a cusp at the first van Hovepoint. The cusp becomes more pronounced with increasing B.Experimentally, as shown in Fig. 2C, there is a cusp-like featuredeveloping at |�| ∼ 0.5 − 0.8 depending on the contact pairwithin the device used, consistent with theoretical predictions.In many contact pairs, this feature presents as a shoulder atB = 0,only developing into a cusp at B ∼ 0.1 T (SI Appendix, Fig. S7).As depicted in Fig. 2B, the calculated filling dependence ofthe Hall number shows two singular sign changes near � ≈ ±2,each falling within inside an open FS region. The filling at whicheach sign-changing singularity occurs is B-independent and isnot directly associated with any van Hove point (see SI Appendix,Fig. S5 for a plot of �H (B), which crosses zero at the samefilling fraction inside the open FS filling range for varying fieldstrength). Moreover, the filling dependence of the Hall numbernH tracks the filling fraction in a broad filling range near theCNP, with the filling range being extended upon increasing B.Note also that the Hall number is generally field-dependent dueto the impact of crystalline symmetries on the band dispersions(See SI Appendix, section IV.A for a detailed analysis). In Fig. 2D,we observe the same general shape of the Hall number. Withinthe open FS filling range, however, the measured Hall numberqualitatively deviates from the theoretical curves. We tentativelyattribute this to a small constant offset in the magnetic field oforder 10 to 20 mT, an amount typical for trapped flux in aNbTi/Nb3Sn superconducting magnet like ours. To explicate,in the open FS regime, we expect (and observe) a large quadraticsymmetric component of the transverse resistivity together witha vanishing antisymmetric component. Hence, quantification6 of 10 https://doi.org/10.1073/pnas.2307151120 pnas.orgDownloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialsA BC DFig. 2. Magnetotransport properties of strained TBG. (A and B) Theoretical calculations of transport properties as a function of magnetic field strengths !0�for � = 1.38◦, � = 0.2% and ' = 0◦. The cyclotron frequency !c defined in Eq. 20 is filling dependent, hence our choice to use the bare cyclotron frequency!0� = eB�/m0. The vertical dashed lines mark the calculated van Hove points, with yellow regions indicating open FSs. (A) Longitudinal MR along the principalaxes ê1 (dashed) and ê2 (solid) in units of �Q Γ�M, where �Q ≡ h/e2 is the quantum of resistance, Γ ≡ ℏ/� is the transport decay rate, and �M ≡ ℏvF |K|� isthe characteristic energy scale for moiré electrons. For a transport rate Γ = 0.1meV, �Q Γ�M≈ 9.6Ω, and !0� ≈ 0.13 is equivalent to a magnetic field strengthB ≈ 0.11T . (B) Hall number nH ≡ B/e�H plotted in units of ns , where 2ns is the total electron density for the narrow bands. (C and D) Experimental measurementsof longitudinal MR (contact pair 14-15) and transverse MR (contact pair 15 - 5) for the TBG sample in ref. 45, which has 20 contacts, at 1.6 K. Vertical dashed linesmark the densities that we ascribe to van Hove points based on the cusp near � ∼ 0.8 and the onset of quadratic MR (shaded yellow). Finite-field resistivitiesin panel (C) are symmetrized: � = (�(B) + �(−B))/2. Panel (D) is calculated from the antisymmetrized transverse resistivity. In the experimental figures, n/nslabels the electron filling fractions per moiré unit cell.of the antisymmetric component involves subtracting two largenumbers. An offset of only a few mT in the assigned field willlead to a small part of the symmetric component mixing intothe antisymmetric component, leading to these deviations fromtheory (SI Appendix, Fig. S8).Our calculation finds a dramatic rotation of the principal axiswith filling, as illustrated in Fig. 3. In the filling range with openFSs, the principal axis with saturating MR (ê1) is aligned withdirection of the shortest moiré triangular bond, suggesting thatthe electrons are hopping more efficiently along the shortestbond, which leads to a larger conductivity and therefore a smallerresistivity. Remarkably, when filling is changed from the secondvan Hove point (� ≈ ±1.3) to the vicinity of the CNP, ê1 rotatesby about 90°. This is likely associated with the larger Fermi pocketencircling a Dirac point being elongated in a direction nearlyperpendicular to the open FS contours; see, for example, Fig. 1D–F. Such rotation of the principal transport axis with fillinghas been attributed to interaction-induced nematicity (14), butwe can now see that in some contexts, it could occur purely dueto strain-induced bandstructure effects. Such a filling-dependentrotation of the principal transport axis was not possible to observein ref. 45 using the Hall bar geometry, where only �xx and �yxare measured but not �yy. Additional transport measurementsare needed, where the filling dependence of the entire resistivitytensor can be mapped out.Since this theory predicts a third van Hove point between theCNP and the filling range with open FSs, a direct experimentalsignature of this van Hove point is desired. In Fig. 4, we reanalyzequantum oscillation measurements of the TBG device discussedin ref. 45. The effective cyclotron mass m∗ is light in the fillingrange with two small closed Fermi pockets and dramaticallyheavier in the filling range with only one closed pocket (Fig. 1D–F and SI Appendix, Fig. S5). The large difference in masses oneither side of the innermost van Hove singularity can account forthe substantially lower-field onset of quantum oscillations closeto the CNP than away from it, as shown in Fig. 4A. Fig. 4B isa Fourier transform of the quantum oscillation data with respectto 1/B. In the filling range of −0.7 ≤ � ≤ 0.8, three distinctfrequencies fi=1,2,3 are clearly observed in the data, with f1 andf2 corresponding to two small Fermi pockets, and f3 = f1 + f2to the breakdown orbit: By 1T, the inverse magnetic length iscomparable to the momentum space distance between the twosmall Fermi pockets (67). Each edge of this filling range is markedby two features: 1. f1 and f2 disappear from the Fourier transform,leaving only f3 at higher electron or hole filling. 2. A cusp-likefeature occurs in longitudinal MR (Fig. 2 A and C ). These bothare predicted by our model as features of a Lifshitz transitionon either side of the CNP, associated with crossing the lowest-energy van Hove points. This detailed match unambiguouslydemonstrates the existence of a third van Hove singularity to eachside of the CNP, at electron or hole filling range (respectively)lower than the onset of B2 MR. The Hall number does not showa sign-changing singularity at this van Hove point, as illustratedin Fig. 2 B and D.PNAS 2023 Vol. 120 No. 34 e2307151120 https://doi.org/10.1073/pnas.2307151120 7 of 10Downloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialsABFig. 3. (A) Rotation of the transport principal axis ê1 with respect to theglobal coordinate system for strained BM with � = 0.2% and ' = 0◦. Thethree horizontal dashed lines are the bond directions. In the open FS region,the saturating MR axis is locked to the shortest bond (L1) direction. However,it rapidly rotates in the closed FS region upon approaching the CNP. (B)Principal transport axes ê1 (red) and ê2 (blue) for a few filling fractions. Nearthe CNP, ê1 is perpendicular to the shortest moiré bond direction. In the openFS filling range (e.g., � ≈ 2.13), it is rotated to be parallel to the shortest bonddirection.The frequencies f1,2 are a strong constraint on the amountof heterostrain in the TBG sample. Specifically, as illustrated inFig. 4C, the frequency f2 is roughly two times f1, showing that thetwo small Fermi pockets have an area ratio∼ 2 : 1. Theoretically,as illustrated by the solid black lines in Fig. 4C, for a heterostrainstrength � = 0.2% and ' = 0◦, the areas Ai=1,2 of the two smallpockets, when converted to frequency via f −1i ≡ (Δ 1B )i = 2�eℏAi ,are in good agreement with experiment. Furthermore, also notethat the frequencies f1,2 extrapolate to 0 at � ≈ ±0.04, showingthat the two Dirac points are shifted to finite (opposite) fillingfractions by heterostrain.We observe behavior qualitatively similar in all respects to thatin Fig. 4A in all 3 longitudinal contact pairs for which we havedilution-fridge measurements (SI Appendix, Fig. S11).4. Summary and OutlookIn summary, we have shown that due to the large size of the moiréunit cell at small twist angles, even a small amount of uniaxialheterostrain on the microscopic scale can lead to dramatic changesin the narrow bands of twisted bilayer graphene. A key featureof the strained bandstructure is the splitting of the respectiveenergetic degeneracies of the two Dirac points and the three vanHove points. The splitting of the two Dirac points leads to asemimetallic state with two small Fermi pockets at the CNP.On the other hand, the two outermost van Hove points bounda broad filling range near � = ±2 where the constant energycontours become open. Interestingly, the elongation of the largerFermi pocket near the CNP is perpendicular to that of the openFSs, the latter being perpendicular to the direction of the shortestmoiré triangular bond.We have analyzed the resulting magnetotransport in strainedTBG in the framework of the Boltzmann equation using themethod of characteristics, treating the magnetic field nonpertur-batively. We showed that a nonsaturating quadratic longitudinalmagnetoresistance in a broad filling range near � = ±2 naturallyarises due to the heterostrain-induced open Fermi surfaces,therefore explaining in remarkable detail experimental results inoff-magic-angle devices with lattice anisotropy (45). We have alsoshown that the sign-changing singularities in the Hall numberoccur in the open FS filling range and are not directly associatedwith any van Hove singularity as commonly assumed, e.g., inref. 68. Furthermore, our theory reveals a dramatic rotationof the transport principal axis as the filling is tuned from thecharge neutrality point to the filling range of open Fermi surfaces,without invoking interaction-induced electronic nematicity.A BCFig. 4. (A) Line cuts of MR near the CNP taken at 26 mK in contact pair 4 to 5 at the indicated field strengths, in Tesla. Vertical dashed lines indicate ourestimated location of the lowest-energy van Hove points, based on the cusps in resistivity at low field. Within the region bounded by these points, the quantumoscillations show up before 0.4 T, and their relative strengths do not follow a simple pattern. Outside of this region, the quantum oscillations onset at higherfield, and every multiple of 4 quantum Hall filling fraction is observed relatively equally. (B) Fourier transform of the quantum oscillation data with respect to 1/B.It reveals a transition from two pockets to one pocket at the lowest-energy van Hove points. (C) Schematic description of the frequencies observed in panel (B).Red dashed lines are frequencies from the experimental data. Solid black lines are predictions from the theory for � = 0.2% and ' = 0◦. The two frequenciesf1 and f2 sum to the one-pocket frequency f3 that extends beyond the first van Hove point. They additionally account for the nontrivial relative strengths of thequantum oscillations within the bounds of the first van Hove points. As with other details of this work, the theory predicts electron–hole symmetry, while someasymmetry is observed in experiment.8 of 10 https://doi.org/10.1073/pnas.2307151120 pnas.orgDownloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialsGiven the importance of energy-shifted van Hove pointsin the transport properties of TBG devices, we have analyzedprevious quantum oscillation data, which has revealed a Lifshitztransition from two pockets to one pocket at a filling fractionwhere the innermost van Hove singularity is predicted to occurbased on theoretical calculations, offering strong evidence ofheterostrain effects on these devices. We have further proposedseveral additional signatures to look for in future experiments.These include a significant difference in cyclotron mass on eitherside of the innermost van Hove singularity (probed qualitativelybut not quantitatively in the extant experiment) and a principaltransport axis with saturating magnetoresistance in the openFermi surface filling range.Finally, given the amplifying effect of a small strain at theunderlying carbon lattice scale on the moiré lattice scale, thelatter of which controls the electronic behavior within the narrowbands, it is tantalizing to consider strain engineering of suchdevices to achieve effects which in regular solids would requireapplying strain magnitudes incompatible with structural stability.Data, Materials, and Software Availability. The data for Fig. 4 and SIAppendix, Fig. S11 were previously published with ref. 45 in ref. 69. Allother experimental data are available at ref. 70. The code for calculatingmagnetotransport of TBG under uniaxial heterostrain is available at ref. 71.Tabular data with analysis notebook data have been deposited in StanfordDigital Repository (https://doi.org/10.25740/zs335dw3715).ACKNOWLEDGMENTS. X.W. acknowledges financial support from NationalMagLab through Dirac fellowship, which is funded by the NSF (Grant No.DMR-1644779) and the state of Florida. O.V. was supported by NSF GrantNo. DMR-1916958 and is partially funded by the Gordon and Betty MooreFoundation’s EPiQS Initiative Grant GBMF11070, National High Magnetic FieldLaboratory through NSF Grant No. DMR-1157490 and the State of Florida. Devicemeasurements and analysis were supported by the US Department of Energy,Office of Science, Basic Energy Sciences, Materials Sciences and EngineeringDivision, under contract DE-AC02-76SF00515. Measurement infrastructurewas funded in part by the Gordon and Betty Moore Foundation’s EPiQSInitiative through grant GBMF3429 and grant GBMF9460. D.G.-G. gratefullyacknowledges support from the Ross M. Brown Family Foundation. SandiaNational Laboratories is a multimission laboratory managed and operated byNational Technology & Engineering Solutions of Sandia, LLC, a wholly ownedsubsidiary of Honeywell International Inc., for the US Department of Energy’sNational Nuclear Security Administration under contract DE-NA0003525.K.W. and T.T. acknowledge support from JSPS KAKENHI (Grant Numbers19H05790, 20H00354, and 21H05233). Part of this work was performedat the Stanford Nano Shared Facilities, supported by the NSF under awardECCS-2026822.Author affiliations: aNational High Magnetic Field Laboratory, Tallahassee, FL 32310;bDepartment of Physics, Stanford University, Stanford, CA 94305; cStanford Institutefor Materials and Energy Sciences, SLAC National Accelerator Laboratory, Menlo Park,CA 94025; dMaterials Physics Department, Sandia National Laboratories, Livermore,CA 94550; eDepartment of Applied Physics, Stanford University, Stanford, CA 94305;fResearch Center for Functional Materials, National Institute for Materials Science,Tsukuba 305-0044, Japan; gInternational Center for Materials Nanoarchitectonics, Na-tional Institute for Materials Science, Tsukuba 305-0044, Japan; hDepartment of Physics,Massachusetts Institute of Technology, Cambridge, MA 02139; and iDepartment ofPhysics, Florida State University, Tallahassee, FL 32306Author contributions: X.W., J.F., M.A.K., O.V., and D.G.-G. designed research; X.W., J.F.,A.L.S., L.K.R., C.L.H., and O.V. performed research; K.W. and T.T. provided hexagonal boronnitride; M.A.K., O.V., and D.G.-G. supervised the experiments and analysis; X.W., J.F., A.L.S.,and O.V. analyzed data; and X.W., J.F., A.L.S., and O.V. wrote the paper.Reviewers: L.B., University of California, Santa Barbara; and E.T., The University of Texasat Austin.1. Y. Cao et al., Correlated insulator behaviour at half-filling in magic-angle graphene superlattices.Nature 556, 80–84 (2018).2. Y. Cao et al., Unconventional superconductivity in magic-angle graphene superlattices. Nature556, 43–50 (2018).3. A. Kerelsky et al., Maximized electron interactions at the magic angle in twisted bilayer graphene.Nature 572, 95–100 (2019).4. X. Lu et al., Superconductors, orbital magnets and correlated states in magic-angle bilayergraphene. Nature 574, 653–657 (2019).5. Y. Jiang et al., Charge order and broken rotational symmetry in magic-angle twisted bilayergraphene. Nature 573, 91–95 (2019).6. M. Yankowitz et al., Tuning superconductivity in twisted bilayer graphene. Science 363,1059–1064 (2019).7. Y. Choi et al., Electronic correlations in twisted bilayer graphene near the magic angle. Nat. Phys.15, 1174–1180 (2019).8. A. L. Sharpe et al., Emergent ferromagnetism near three-quarters filling in twisted bilayergraphene. Science 365, 605–608 (2019).9. Y. Xie et al., Spectroscopic signatures of many-body correlations in magic-angle twisted bilayergraphene. Nature 572, 101–105 (2019).10. U. Zondiner et al., Cascade of phase transitions and dirac revivals in magic-angle graphene. Nature582, 203–208 (2020).11. D. Wong et al., Cascade of electronic transitions in magic-angle twisted bilayer graphene. Nature582, 198–202 (2020).12. M. Serlin et al., Intrinsic quantized anomalous hall effect in a moiré heterostructure. Science 367,900–903 (2020).13. P. Stepanov et al., Untying the insulating and superconducting orders in magic-angle graphene.Nature 583, 375–378 (2020).14. Y. Cao et al., Nematicity and competing orders in superconducting magic-angle graphene. Science372, 264–271 (2021).15. X. Liu et al., Tuning electron correlation in magic-angle twisted bilayer graphene using coulombscreening. Science 371, 1261–1265 (2021).16. A. T. Pierce et al., Unconventional sequence of correlated Chern insulators in magic-angle twistedbilayer graphene. Nat. Phys. 17, 1210–1215 (2021).17. S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, E. Y. Andrei, Chern insulators, van Hove singularitiesand topological flat bands in magic-angle twisted bilayer graphene. Nat. Mater. 20, 488–494(2021).18. J. Yu et al., Correlated Hofstadter spectrum and flavour phase diagram in magic-angle twistedbilayer graphene. Nat. Phys. 18, 825–831 (2022).19. R. Bistritzer, A. H. MacDonald, Moiré bands in twisted double-layer graphene. Proc. Natl. Acad. Sci.U.S.A. 108, 12233–12237 (2011).20. E. Y. Andrei, A. H. MacDonald, Graphene bilayers with a twist. Nat. Mater. 19, 1265–1275 (2020).21. L. Balents, C. R. Dean, D. K. Efetov, A. F. Young, Superconductivity and strong correlations in moiréflat bands. Nat. Phys. 16, 725–733 (2020).22. M. Koshino et al., Maximally localized Wannier orbitals and the extended Hubbard model fortwisted bilayer graphene. Phys. Rev. X 8, 031087 (2018).23. H. C. Po, L. Zou, A. Vishwanath, T. Senthil, Origin of Mott insulating behavior and superconductivityin twisted bilayer graphene. Phys. Rev. X 8, 031089 (2018).24. J. Kang, O. Vafek, Symmetry, maximally localized Wannier states, and a low-energy model fortwisted bilayer graphene narrow bands. Phys. Rev. X 8, 031088 (2018).25. J. Kang, O. Vafek, Strong coupling phases of partially filled twisted bilayer graphene narrow bands.Phys. Rev. Lett. 122, 246401 (2019).26. Z. Song et al., All magic angles in twisted bilayer graphene are topological. Phys. Rev. Lett. 123,036401 (2019).27. M. Xie, A. H. MacDonald, Nature of the correlated insulator states in twisted bilayer graphene.Phys. Rev. Lett. 124, 097601 (2020).28. N. Bultinck et al., Ground state and hidden symmetry of magic-angle graphene at even integerfilling. Phys. Rev. X 10, 031034 (2020).29. Y. Zhang, K. Jiang, Z. Wang, F. Zhang, Correlated insulating phases of twisted bilayer graphene atcommensurate filling fractions: A Hartree-Fock study. Phys. Rev. B 102, 035136 (2020).30. T. Cea, F. Guinea, Band structure and insulating states driven by coulomb interaction in twistedbilayer graphene. Phys. Rev. B 102, 045107 (2020).31. J. Kang, O. Vafek, Non-abelian dirac node braiding and near-degeneracy of correlated phasesat odd integer filling in magic-angle twisted bilayer graphene. Phys. Rev. B 102, 035161(2020).32. O. Vafek, J. Kang, Renormalization group study of hidden symmetry in twisted bilayer graphenewith coulomb interactions. Phys. Rev. Lett. 125, 257602 (2020).33. B. Lian et al., Twisted bilayer graphene. IV. Exact insulator ground states and phase diagram. Phys.Rev. B 103, 205414 (2021).34. B. A. Bernevig et al., Twisted bilayer graphene. V. Exact analytic many-body excitations in coulombHamiltonians: Charge gap, goldstone modes, and absence of cooper pairing. Phys. Rev. B 103,205415 (2021).35. F. Xie et al., Twisted bilayer graphene. VI. An exact diagonalization study at nonzero integer filling.Phys. Rev. B 103, 205416 (2021).36. P. Potasz, M. Xie, A. H. MacDonald, Exact diagonalization for magic-angle twisted bilayer graphene.Phys. Rev. Lett. 127, 147203 (2021).37. Y. H. Kwan et al., Kekulé spiral order at all nonzero integer fillings in twisted bilayer graphene.Phys. Rev. X 11, 041063 (2021).38. D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, N. Bultinck, Strain-induced quantum phasetransitions in magic-angle graphene. Phys. Rev. Lett. 127, 027601 (2021).39. L. Huder et al., Electronic spectrum of twisted graphene layers under heterostrain. Phys. Rev. Lett.120, 156405 (2018).40. Z. Bi, N. F. Q. Yuan, L. Fu, Designing flat bands by strain. Phys. Rev. B 100, 035448 (2019).41. F. Mesple et al., Heterostrain determines flat bands in magic-angle twisted graphene layers. Phys.Rev. Lett. 127, 126405 (2021).PNAS 2023 Vol. 120 No. 34 e2307151120 https://doi.org/10.1073/pnas.2307151120 9 of 10Downloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://www.pnas.org/lookup/doi/10.1073/pnas.2307151120#supplementary-materialshttps://doi.org/10.25740/zs335dw371542. N. P. Kazmierczak et al., Strain fields in twisted bilayer graphene. Nat. Mater. 20, 956–963 (2021).43. K. Kim et al., van der Waals heterostructures with high accuracy rotational alignment. Nano Lett.16, 1989–1995 (2016).44. Y. Cao et al., Superlattice-induced insulating states and valley-protected orbits in twisted bilayergraphene. Phys. Rev. Lett. 117 (2016).45. J. Finney et al., Unusual magnetotransport in twisted bilayer graphene. Proc. Natl. Acad. Sci. U.S.A.119, e2118482119 (2022).46. H. Suzuura, T. Ando, Phonons and electron-phonon scattering in carbon nanotubes. Phys. Rev. B65, 235412 (2002).47. N. N. T. Nam, M. Koshino, Lattice relaxation and energy band modulation in twisted bilayergraphene. Phys. Rev. B 96, 075311 (2017).48. C. Kittel, Quantum Theory of Solids (Wiley, New York, 1963).49. O. Vafek, J. Kang, Continuum effective Hamiltonian for graphene bilayers for an arbitrary smoothlattice deformation from microscopic theories (2022).50. J. Kang, O. Vafek, Pseudo-magnetic fields, particle-hole asymmetry, and microscopic effectivecontinuum Hamitonians of twisted bilayer graphene (2022).51. E. H. Hwang, S. Das Sarma, Acoustic phonon scattering limited carrier mobility in two-dimensionalextrinsic graphene. Phys. Rev. B 77, 115449 (2008).52. D. K. Efetov, P. Kim, Controlling electron-phonon interactions in graphene at ultrahigh carrierdensities. Phys. Rev. Lett. 105, 256805 (2010).53. K. Kaasbjerg, K. S. Thygesen, K. W. Jacobsen, Unraveling the acoustic electron-phonon interactionin graphene. Phys. Rev. B 85, 165440 (2012).54. D. Grassano et al., Work function, deformation potential, and collapse of Landau levels in strainedgraphene and silicene. Phys. Rev. B 101, 245115 (2020).55. J. M. B. Lopes, N. M. R. dos Santos, Peres, A. H. Castro Neto, Graphene bilayer with a twist:Electronic structure. Phys. Rev. Lett. 99, 256802 (2007).56. B. A. Bernevig, Z. D. Song, N. Regnault, B. Lian, Twisted bilayer graphene. I. Matrix elements,approximations, perturbation theory, and a k · p two-band model. Phys. Rev. B 103, 205411(2021).57. J. Kang, B. A. Bernevig, O. Vafek, Cascades between light and heavy fermions in thenormal state of magic-angle twisted bilayer graphene. Phys. Rev. Lett. 127, 266402(2021).58. F. Guinea, N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twistedgraphene bilayers. Proc. Natl. Acad. Sci. U.S.A. 115, 13174–13179 (2018).59. L. Rademaker, D. A. Abanin, P. Mellado, Charge smoothening and band flattening due to Hartreecorrections in twisted bilayer graphene. Phys. Rev. B 100, 205114 (2019).60. Z. A. H. Goodwin, V. Vitale, X. Liang, A. A. Mostofi, J. Lischner, Hartree theory calculations ofquasiparticle properties in twisted bilayer graphene. Elect. Struct. 2, 034001 (2020).61. Y. Choi et al., Interaction-driven band flattening and correlated phases in twisted bilayer graphene.Nat. Phys. 17, 1375–1381 (2021).62. M. Audin, M. Damian, Morse Theory and Floer Homology (Springer, 2014).63. M. O. Goerbig, J. N. Fuchs, G. Montambaux, F. Piéchon, Tilted anisotropic dirac cones in quinoid-type graphene and � -(BEDT-TTF)2i3 . Phys. Rev. B 78, 045415 (2008).64. A. A. Soluyanov et al., Type-ii Weyl semimetals. Nature 527, 495–498 (2015).65. M. Xie, A. H. MacDonald, Weak-field hall resistivity and spin-valley flavor symmetry breaking inmagic-angle twisted bilayer graphene. Phys. Rev. Lett. 127, 196401 (2021).66. I. Lifshits, M. Azbel, M. Kaganov, Electron Theory of Metals (Springer, 1973).67. D. Schoenberg, Magnetic Oscillations in Metals (2009).68. J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, P. Jarillo-Herrero, Tunable strongly coupledsuperconductivity in magic-angle twisted trilayer graphene. Nature 590, 249–255 (2021).69. J. Finney et al., Data for: Unusual magnetotransport in twisted bilayer graphene, Stanford DigitalRepository (2021). https://doi.org/10.25740/tm725vs8229.70. X. Wang et al., Data for: Unusual magnetotransport in twisted bilayer graphene from strain-induced open fermi surfaces, Stanford Digital Repository (2023). https://doi.org/10.25740/zs335dw3715.71. X. Wang et al., Code for: Unusual magnetotransport in twisted bilayer graphene from strain-induced open fermi surfaces (2033). https://github.com/xywang2017/TBG_MR.10 of 10 https://doi.org/10.1073/pnas.2307151120 pnas.orgDownloaded from https://www.pnas.org by BUSSHITSU ZAIRYO KENKYUKIKO LIBRARY on August 18, 2023 from IP address 144.213.253.16.https://doi.org/10.25740/tm725vs8229https://doi.org/10.25740/zs335dw3715https://doi.org/10.25740/zs335dw3715https://github.com/xywang2017/TBG_MR 1. Geometric and Energetic Effects of Uniaxial Heterostrain on TBG 2. Boltzmann Equation and Magnetoresistivity in TBG 3. Comparison between Theory and Experiment 4. Summary and Outlook