Unusual magnetotransport in twisted bilayer graphene from strain-induced open Fermi surfaces

Significance Because of its rich array of correlated phases, twisted bilayer graphene (TBG) near the magic angle has captivated the condensed matter physics world. The large moiré length scale not only promotes interaction-related effects but also allows for extrinsic factors such as strain to play a major role. In a previous work, we presented measurements of a TBG device with several unusual behaviors in magnetotransport and conjectured that uniaxial strain could explain our measurements. Here, we model magnetotransport in TBG by incorporating uniaxial heterostrain into the Bistritzer–MacDonald Hamiltonian. The theory not only reproduces the unusual phenomena from the previous work but also predicts additional features unnoticed before. Our work therefore demonstrates the crucial role of heterostrain in TBG devices.

Strain-especially differing lattice distortions in the two layers, termed heterostrain-is believed to play an important role in the phase diagram of TBG (37)(38)(39)(40)(41). Scanning probe measurements 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, as opposed to strain applied equally to both layers, the linear distortion of the moiré unit cell is amplified by a factor of ∼ 1/ relative to the linear distortion of the microscopic atomic lattice. For example, 0.2% uniaxial heterostrain causes an ∼ 8% change in the largest linear dimension of the moiré unit cell for a twist angle of 1.38 • . The moiré unit cell 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 underestimates of 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 a moiré unit cell area of 90 nm 2 (corresponding to = 1.38 • ) displayed several unusual phenomena 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 nearmagic-angle devices. Surprisingly, though, over a broad filling range near half-filling the longitudinal magnetoresistivity (MR) exhibited a B 2 increase up to ∼ 100-fold at ≈ 5 T, after which quantum oscillations set in. The authors found that a toy Hofstadter model with anisotropy showed multiple features similar to those in the data, and they conjectured that uniaxial strain might cause such anisotropy.
In this work, we present a systematic theoretical study of the impact of uniaxial heterostrain on the narrow-band dispersion of TBG above the magic angle, analyze its consequences for weak field magnetotransport, and compare it with experimental data from ref. 45. We base our theory on the Bistritzer-MacDonald (BM) continuum

Significance
Because of its rich array of correlated phases, twisted bilayer graphene (TBG) near the magic angle has captivated the condensed matter physics world. The large moiré length scale not only promotes interaction-related effects but also allows for extrinsic factors such as strain to play a major role. In a previous work, we presented measurements of a TBG device with several unusual behaviors in magnetotransport and conjectured that uniaxial strain could explain our measurements. Here, we model magnetotransport in TBG by incorporating uniaxial heterostrain into the Bistritzer-MacDonald Hamiltonian. The theory not only reproduces the unusual phenomena from the previous work but also predicts additional features unnoticed before. Our work therefore demonstrates the crucial role of heterostrain in TBG devices. model (19), incorporating heterostrain in the form of a deformation potential, a pseudomagnetic field (46,47), and a distortion of the moiré pattern.
Our key theoretical result is that heterostrain lifts the energetic degeneracy of the two Dirac points as well as that of the three van Hove points of a given band. The splitting of the two Dirac points leads to a semimetallic state with small Fermi pockets near the charge neutrality point (CNP). More interestingly, the splitting 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 Boltzmann equation, the open FSs generally lead to a nonsaturating B 2 MR (48), accounting for this previously unexplained feature in the experimental data of ref. 45.
This theory makes a number of falsifiable predictions. 1) If the direction of current flow in the lithographically patterned Hall bar is not aligned with the 1D-like principal axis of transport in the open FS regime, longitudinal magnetoresistance should mix substantially into the measured resistance at Hall contacts.
2) A subtle cusp should appear in resistivity as filling crosses the lowest-energy van Hove point. 3) A Lifshitz transition from two FS pockets to one should also coincide with crossing this lower van Hove point. We reanalyze experimental data from ref. 45 and find that these predictions are verified. The strained BM model studied here has electron-hole symmetry. We leave the discussion of electron-hole asymmetry in the experimental data to future works.
The theory also has implications that are not so far directly probed by the experiment, but are striking. First, on each side of the CNP, the divergence and sign change in the Hall number near half-filling of a 4-fold degenerate band does not coincide with any of the van Hove points but instead occurs within the filling range where FSs are open. Second, the transport principal axis continuously rotates by up to 90 • as density is tuned from the CNP to the open FS regime. Such rotation of the transport axes has been presented as evidence for interaction-induced nematic order (14), but here, we find that it can arise purely due to strain-induced band structure effects.
This work clearly demonstrates that the effects of even miniscule amounts of heterostrain in TBG cannot be neglected. Dramatic and unexpected phenomena occur in strained TBG even in the single-particle regime, without the strong correlation effects that arise near the magic angle. Given the amplifying effect of a small heterostrain on the moiré length scale, it is tantalizing to consider strain engineering of such devices to achieve effects that would be impossible in regular solids due to structural instabilities.
The paper is organized as follows: In Section I, we present the theoretical calculation of the changes in band structure resulting from uniaxial strain. In Section II, we use the Boltzmann transport equation to calculate the electron transport properties in magnetic fields resulting from the strain-induced alterations in the band structure. Section III is a comparison of our predictions with experiment, and it is sufficiently self-contained that a reader who does not need all the theoretical details can jump directly there. In Section IV, we summarize our results.

Geometric and Energetic Effects of Uniaxial Heterostrain on TBG
In the limit of small deformations, both uniaxial heterostrain and a small twist angle are captured via a coordinate transformation: r l = r + u l (r), where l = t, b labels the Top (Bottom) graphene layers, and u l (r) ≈ E l r is the local deformation field. The symmetric and antisymmetric part of the 2×2 tensor E l describes strain and rotation, respectively. For twist angle ( ) and a uniaxial heterostrain of strength ( ) and direction ( ), we parameterize , and given by: Here, R is the two-dimensional rotation matrix, and ≈ 0.16 is the Poisson ratio (3). Physically, > 0 corresponds to compressing the Top layer while stretching the Bottom layer along the direction determined by , as illustrated in Fig. 1A. A relative deformation E between the graphene bilayers generates a moiré superlattice, with moiré reciprocal lattice vectors 2 , where G i are reciprocal lattice vectors of the undeformed monolayer graphene. The moiré lattice vectors L i=1,2 are uniquely defined through the relation L i · g j = 2 ij . It is important to note that only relative deformations generate the moiré superlattice. Deformations that act homogeneously on the two layers have little effect on the narrow band physics, and we accordingly neglect them in this work * .
Under rotation R , the strain tensor transforms as a headless vector that remains invariant under → + 180 • . Combined with the C 3z symmetry of the undeformed graphene lattice, the strained electronic dispersion within a given graphene valley simply rotates 60 • under → + 60 • . We therefore only report results for ∈ [0 • , 60 • ). For concreteness, we define the microscopic unit cell vectors a i=1,2 of an undeformed graphene lattice as a 1 Fig. 1 A and B illustrates the geometric effects of heterostrain for twist angle = 1.38 • . At = 0, AA-stacked regions of the moiré form an equilateral triangular superlattice. Introducing uniaxial heterostrain changes the spacings between neighboring AA-stacked regions (L i=1,2,3 ). For = 0.2%, typical in these systems (3,9,11), the variation in spacings can be as large as / ≈ 8%. Such dramatic amplification of the microscopic strain makes moiré materials uniquely suited to strain engineeringconventional materials become structurally unstable at distortions only 10% as large as those achieved in the moiré superlattice. Note that the effect of uniaxial heterostrain on the moiré unit cell area is small at 2 2 / 2 , as some spacings become larger while others become smaller (SI Appendix, section I).
We proceed to discuss the energetic effects in the context of the continuum BM model (19). We work in the limit where both E l and the wavevector k in the moiré Brillouin zone are small and consider only the leading order terms in both. This would mean, for instance, that terms such as Ek are omitted as higherorder terms. This treatment is generally justified away from the magic angle because higher-order terms can play an important role only close to the magic angle where the bandwidth becomes comparable to these terms (49,50). We explicitly checked that * We checked numerically that adding a small homogeneous strain in addition to a heterostrain of similar strength yields band and transport properties almost identical to those produced by adding a heterostrain alone. at ≈ 1.38 • , the effects of such higher-order terms are indeed negligibly small. To leading order, the strained BM Hamiltonian for a given valley is given by where = ±1 labels K (K ) valleys of monolayer graphene. The interlayer Hamiltonian is given by ,l (r) ≡ ( ,l,A (r), ,l,B (r)) T is a spinor in the sublattice basis for a given valley and layer. We have suppressed the spin index for simplicity. q j=1,2,3 are the three nearest neighbor bonds of the reciprocal honeycomb lattice, and [4] ( 0 , x , y ) are Pauli matrices acting on sublattice degrees of freedom.
The intralayer Hamiltonian is given by: ,l (r).

[5]
Here, the first term is the deformation potential that couples to the electron density. Its value is not precisely known in the literature, with numbers ranging from −4.1 eV to 30 eV depending on the methodology (51-54). We use = −4.1 eV in this work based on first principles calculations (54), although for heterostrain ≈ 0.2% varying the deformation potential over the range proposed in the literature leads to only minor quantitative differences in band dispersions. A ,l is the pseudovector potential that comes from changes in the intersublattice hopping due to deformations and changes sign between graphene valleys. It is given 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 ℏv F /a = 2.68eV based on Fermi velocity in monolayer graphene v F ≈ 10 6 m/s (55), w 0 = 88meV, and w 1 = 110meV (19) in our calculations and also set ℏ = 1 in the remainder of the paper.
To leading order approximation, the strained BM Hamiltonian in a given valley (Eq. 2) has particle-hole symmetry under P l (r) = l i( y ) ll l (−r) (56), where y is a Pauli matrix acting on the layer degrees of freedom. This means that for every single-electron state at energy E and wavevector k, there is a state at energy −E and wavevector −k. This particle-hole symmetry has been investigated extensively for the unstrained BM model, e.g., refs. 26 and 57, and here, it is generalized to the strained case. The source of experimentally evident particle-hole asymmetry in the off-magic-angle device (45) could be higherorder gradient terms beyond what is captured in the BM model in 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-axis defined in Fig. 1A, as specified by ∈ [0 • , 60 • ). For simplicity, we show only contour maps of the upper band from valley K in the moiré Brillouin zone specified by k = k 1 g 1 + k 2 g 2 , where k 1,2 ∈ [0, 1). Heterostrain preserves C 2 T and valley U (1) (23), so the Lower and Upper bands remain connected via two Dirac points. The Upper band features six special points-two Dirac points (black stars), three van Hove points (colored dots), and one band maximum (black cross). The six special points of a given band are related to "critical points" in the context of the Morse theory, which states that where i is the index of the i-th critical point, and is the Euler characteristic of a manifold (62); vanishes for the Brillouin zone which is a torus. Although a Dirac point is strictly a point of nonanalyticity and is not directly covered by Morse theory, if we imagine adding a tiny gap term, it will become a legitimate band extremum allowing Morse theory to apply. Whereas the two band minima (Dirac points) and the band maximum have even , and so each contributes +1 to the sum, every conventional van Hove point (i.e., not a higher order) has an odd and contributes −1. The overall sum thus vanishes. Therefore, the van Hove points can only be annihilated/created by colliding with local minima/maxima. For a relatively small heterostrain as shown in Fig. 1, the number 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 aforementioned collisions and the appearance of tilted type II Dirac cones (63,64).
A key finding of the present work is that the respective energy degeneracies of the two Dirac points and the three van Hove points are lifted by uniaxial heterostrain, by amounts depending sensitively on . In the absence of strain, Fig. 1C, the three van Hove points are at equal energy, and separate closed contours of constant energy centered around the Dirac points from closed contours centered around the band maximum. As illustrated in Fig. 1 D-F, uniaxial heterostrain splits the energy degeneracy of the two Dirac points, leading to a semimetallic state with small Fermi pockets near the CNP (40). The three van Hove points also split in energy. The two outermost van Hove points (i.e., closer to the band maximum) bound a filling range of open FSs near = 2, whereas the innermost van Hove point moves closer to one of the Dirac points. If we continue increasing , a collision of the critical points occurs, the innermost van Hove disappears, the two Dirac points become type II tilted, and a new ordinary minimum is created. Note that a small mass added to type II tilted Dirac points will not introduce band extrema, and as a consequence, type II tilted Dirac points are not critical points of Morse theory. Therefore, after the collision, Eq. 6 still holds.
Interestingly, the elongation of the FSs shows a strong filling dependence. Close to the CNP, the bigger Fermi pocket that encloses a Dirac point is stretched along a direction perpendicular to that of the open FSs; see Figs. 1 D-F. As explained later, this leads to a dramatic rotation of the principal transport axis when the filling is tuned from the CNP to the open FS range.
The dependence of the energy and filling of the band structure special points on at a fixed is shown in Fig. 1 G and H. Of notable interest is the sensitivity of the filling range with open FSs to . This filling range must in fact vanish at some between 0 • and 60 • , when the energies of the two outermost van Hove points cross. As seen in Fig. 1 D-F, this also alters the elongation of the open FSs.

Boltzmann Equation and Magnetoresistivity in TBG
Having understood the heterostrain effects on the bandstructure, we proceed to discuss the implications for magnetotransport. We begin by considering the general structure of the two-dimensional resistivity tensor subject to heterostrain. The resistivity tensor is defined via where E = (E x , E y ) T and j = (j x , j y ) T are electric field and current vectors, respectively. Under rotation by , the resistivity tensor transforms as: If the underlying system has a point group symmetry higher than C 2z (e.g., C 3z , C 6z ), then Since heterostrain breaks the point group symmetry down to C 2z , we generally expect xx = yy , xy = − yx . Nevertheless, it is always possible to define transport principal axes after a suitable rotation of the coordinate system, such that principal = Here, 1,2 are longitudinal resistivities along the principal transport directionsê 1,2 , respectively. The rotation angle is determined up to 180 • by requiring 1 < 2 .
Below, we first derive the MR tensor using a Boltzmann approach for a general noninteracting electronic system within the relaxation time approximation. Since there is currently limited understanding of the scattering mechanisms determining electrical transport in TBG, here, we follow ref. 65 and use the relaxation time approximation. We then present the results for heterostrained TBG, showing that in the open FS region, the low-resistivity principal axis (ê 1 ) is nearly perfectly aligned with the shortest moiré bond direction. However, there is a dramatic rotation of the principal axis as the filling moves toward the CNP. We further show that the open FSs lead to a B 2 nonsaturating MR alongê 2 and a saturating resistivity alongê 1 . For random orientation ( 0 ) of the principal axis to the electrical current axis in the Hall bar geometry, e.g., as in ref. 45, the longitudinal resistivity is given by: xx = 1 cos 2 0 + 2 sin 2 0 . This is dominated by the 2 ∼ B 2 component. As a result, the experimental measurements should observe the nonsaturating MR component if there is a misalignment of the Hall bar orientation with respect to the principal transport axis. Such misalignment is generically to be expected and indeed must occur over most of the relevant filling range since the principal axes rotate with filling while the device geometry remains fixed.

Boltzmann Equation and
Method of Characteristics. We begin with a brief description of the method of characteristics used to solve the Boltzmann equation perturbatively in electric field E but without a restriction on the strength of the perpendicular magnetic field B = Bẑ, as long as the semiclassical regime holds (66). Due to C 2z T symmetry of TBG at B = 0, there is no Berry curvature contribution to the semiclassical equations of motion. Then, within the relaxation time approximation, the Boltzmann equation for a given energy band becomes where qE + qv k × B is the total force on the Bloch electrons, with v k ≡ ∇ k k and charge q; n 0,k is the equilibrium Fermi-Dirac distribution and n k is the desired nonequilibrium distribution function.
We consider a stationary solution to the Boltzmann equation by parameterizing the distribution function as: n k = n 0,k + n 1,k . [11] As a result, the Boltzmann equation for the deviation of the distribution function from equilibrium is: Note that the magnetic field only couples to n 1 since (qv k × B) · ∇ k n 0,k = (qv k × B) · v k ∂ k n 0,k = 0. To solve the above partial differential equation (PDE), we seek a family of curves covering the k-space which we parameterize as k(s) with s ∈ [0, s 0 ), such that along these curves, the PDE becomes an ordinary differential equation (ODE). If a curve k(s) satisfies d k(s) ds = qv(s) × B, [13] then n 1,k(s) ≡ n 1 (s) satisfies the curve k(s) must coincide with the contour of constant energy. Thus, the Boltzmann equation becomes: The ODE is readily solved with where 0 is a constant determined by the following argument. Since k(s) describes a constant energy contour in a twodimensional Brillouin zone, it is either a closed contour or several open contours that terminate on boundaries of the Brillouin zone such that they form a closed loop on a torus. In either case, k(s) is periodic under s → s + s 0 modulo a moiré reciprocal lattice vector, where s 0 is the periodicity. The periodicity condition n 1 (s 0 ) = n 1 (0) leads to which determines the desired n 1 (s). In the low-temperature limit, the steady-state current from a given energy band is calculated as: where ( , ) = x, y, and we have defined the cyclotron frequency as: We have also made use of the periodicity of velocity under s → s + s 0 to write it in terms of Fourier series, v(s) = ∞ n=−∞ v n e −in c s . To show that the second line of Eq. 19 holds, note that at every k, we can define a local coordinate system (ê v ,ê s ) such that v ≡ vê v where v ≥ 0, andê s =ê v ×ẑ. The infinitesimal wavevector can be equivalently written as: Eq. 13 can then be written as dk/ds = qvBê s , or equivalently dk s = qvBds. As a result, The conductivity tensor is therefore given by the following expression: Eq. 21 gives the magnetoconductivity for a given FS contour. In the case of multiple FS contours and multiple bandsassociated for example with spin and valley degeneracy in TBGconductivities 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 labels respectively for a given energy level.
To better understand Eq. 21, consider an example of a parabolic dispersion with k = 1 2m 0 (k 2 x + k 2 y ), where m 0 is the bare electron mass. At a fixed energy , the contour is a circle parameterized as (k x , k y ) = √ 2m 0 (cos , sin ), ∈ [0, 2 ). Using the method of characteristics, we get d ds = − Note that the total number density of filled electrons is given by n = d 2 k (2 ) 2 Θ( − k ) = m 0 2 . We therefore reproduce the well-known magnetoconductivity tensor: In this simple example of a closed FS, the longitudinal resistivity is given by m 0 nq 2 , independent of the magnetic field. The average of the velocity field, v n=0 ≡ 1 s 0 s 0 0 dsv(s), vanishes. However, for an open FS, generally, v n=0 = 0, i.e., electrons have a finite drift velocity when a magnetic field causes them to traverse the contour (SI Appendix, Fig. S2). The impact of such a finite drift velocity on the magnetotransport can be qualitatively understood using the following example: In the expression for the conductivity tensor (Eq. 21), we consider v x n=0 = 0 but v y n=0 = 0. This corresponds to an open FS with a drift velocity along the x direction. In the high-field limit ( c ∝ B 1), we truncate the Fourier series at the leading order, yielding where we made use of the equality: v −n = v * n . Inverting the matrix, we obtain the MR tensor:

[25]
We see that for an open FS, the longitudinal MR has nonsaturating B 2 behavior along the axis with a zero drift velocity (ŷ in the above example), and saturating behavior (constant in B in this simple model) along the other axis.

Comparison between Theory and Experiment
As we will see, the theory satisfactorily explains the weak-field magnetotransport measurements presented in ref. 45. We then present two predictions of the theory that we did not anticipate prior to starting this work: the dependence of the principal axis of transport on filling, and the behavior of magnetoresistance and quantum oscillations at densities between the CNP and the onset of quadratic MR. The former has important implications, but it cannot be confirmed with our present datasets because of limitations of the Hall bar geometry. The latter can be considered smoking gun evidence for the presence of the lowest-energy van Hove point and the energetic splitting of the Dirac cones.
We present calculations for = 1.38 • (independently measured for the region of the experimental device we focus on), and = 0.2% and = 0 • , parameters which are not measured in the experiment but are chosen to yield reasonable quantitative agreement between the theoretical and experimental results with respect to the filling range of open FSs and to the frequencies of magnetoresistance oscillations to be presented later. The general phenomena of open FSs and quadratic MR hold for a broad range of heterostrain parameters and . We do not perform fine-tuning of these input parameters for two reasons: 1. We do not expect our strained BM model in Eq. 2 to yield precise quantitative agreement with experiment. Specifically, the model has particle-hole symmetry, which is absent from experimental measurements. More sophisticated noninteracting model calculations (49,50) as well as interaction renormalizations (61) would likely be necessary to properly account for such details. 2. Increasing broadens the filling range that displays open FSs, but for a given , varying strongly tunes this range ( Fig. 1 G and H ), so there is some flexibility in assigning the two parameters to match the experimental measurements. This might be overcome by also seeking to match the position of the lowest van Hove singularity and the evolution of sizes of the two Fermi pockets near the CNP, but the simplifications in the BM model noted in (1) above caution us against drawing strong quantitative conclusions from the comparison.
In Fig. 2, we show the computed MR along the principal transport axes (A) and the Hall number (B). For comparison, we plot the experimentally measured longitudinal and transverse resistivities (C) and Hall number (D) for the TBG device studied in 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 by the two outermost van Hove points of the zero-field strained band structure. In experiment, we observe quadratic MR in longitudinal resistivity within a similar range of fillings. More strikingly, we observe quadratic MR in the transverse resistivity as well. In some cases, with increasing field, the symmetric part of the transverse resistivity becomes larger than that of the longitudinal resistivity. As discussed earlier, this degree of mixing can be attributed to the misalignment between the strain-induced but filling-dependent principal axis of transport and the direction of current flow in the Hall bar geometry.
At the first van Hove point ( ≈ ±0.6), the nonanalyticity in the density of states leads to a cusp in the first derivative of the zero-field resistivity with respect to filling (SI Appendix, Fig. S6). As shown in Fig. 2A, at B = 0 the longitudinal resistance as a function of filling develops a cusp at the first van Hove point. The cusp becomes more pronounced with increasing B. Experimentally, as shown in Fig. 2C, there is a cusp-like feature developing at | | ∼ 0.5 − 0.8 depending on the contact pair within the device used, consistent with theoretical predictions. In many contact pairs, this feature presents as a shoulder at B = 0, only developing into a cusp at B ∼ 0.1 T (SI Appendix, Fig. S7).
As depicted in Fig. 2B, the calculated filling dependence of the Hall number shows two singular sign changes near ≈ ±2, each falling within inside an open FS region. The filling at which each sign-changing singularity occurs is B-independent and is not directly associated with any van Hove point (see SI Appendix, Fig. S5 for a plot of H (B), which crosses zero at the same filling fraction inside the open FS filling range for varying field strength). Moreover, the filling dependence of the Hall number n H tracks the filling fraction in a broad filling range near the CNP, with the filling range being extended upon increasing B. Note also that the Hall number is generally field-dependent due to 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. Within the open FS filling range, however, the measured Hall number qualitatively deviates from the theoretical curves. We tentatively attribute this to a small constant offset in the magnetic field of order 10 to 20 mT, an amount typical for trapped flux in a NbTi/Nb 3  Our calculation finds a dramatic rotation of the principal axis with filling, as illustrated in Fig. 3. In the filling range with open FSs, the principal axis with saturating MR (ê 1 ) is aligned with direction of the shortest moiré triangular bond, suggesting that the electrons are hopping more efficiently along the shortest bond, which leads to a larger conductivity and therefore a smaller resistivity. Remarkably, when filling is changed from the second van Hove point ( ≈ ±1.3) to the vicinity of the CNP,ê 1 rotates by about 90°. This is likely associated with the larger Fermi pocket encircling a Dirac point being elongated in a direction nearly perpendicular to the open FS contours; see, for example, Fig. 1 D-F. Such rotation of the principal transport axis with filling has been attributed to interaction-induced nematicity (14), but we can now see that in some contexts, it could occur purely due to strain-induced bandstructure effects. Such a filling-dependent rotation of the principal transport axis was not possible to observe in ref. 45 using the Hall bar geometry, where only xx and yx are measured but not yy . Additional transport measurements are needed, where the filling dependence of the entire resistivity tensor can be mapped out.
Since this theory predicts a third van Hove point between the CNP and the filling range with open FSs, a direct experimental signature of this van Hove point is desired. In Fig. 4, we reanalyze quantum oscillation measurements of the TBG device discussed in ref. 45. The effective cyclotron mass m * is light in the filling range with two small closed Fermi pockets and dramatically heavier in the filling range with only one closed pocket ( Fig. 1  D-F and SI Appendix, Fig. S5). The large difference in masses on either side of the innermost van Hove singularity can account for the substantially lower-field onset of quantum oscillations close to the CNP than away from it, as shown in Fig. 4A. Fig. 4B is a Fourier transform of the quantum oscillation data with respect to 1/B. In the filling range of −0.7 ≤ ≤ 0.8, three distinct frequencies f i=1,2,3 are clearly observed in the data, with f 1 and f 2 corresponding to two small Fermi pockets, and f 3 = f 1 + f 2 to the breakdown orbit: By 1T, the inverse magnetic length is comparable to the momentum space distance between the two small Fermi pockets (67). Each edge of this filling range is marked by two features: 1. f 1 and f 2 disappear from the Fourier transform, leaving only f 3 at higher electron or hole filling. 2. A cusp-like feature occurs in longitudinal MR (Fig. 2 A and C). These both are predicted by our model as features of a Lifshitz transition on either side of the CNP, associated with crossing the lowestenergy van Hove points. This detailed match unambiguously demonstrates the existence of a third van Hove singularity to each side of the CNP, at electron or hole filling range (respectively) lower than the onset of B 2 MR. The Hall number does not show a sign-changing singularity at this van Hove point, as illustrated in Fig. 2 B and D. A B Fig. 3. (A) Rotation of the transport principal axisê 1 with respect to the global coordinate system for strained BM with = 0.2% and = 0 • . The three horizontal dashed lines are the bond directions. In the open FS region, the saturating MR axis is locked to the shortest bond (L 1 ) 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. Near the CNP,ê 1 is perpendicular to the shortest moiré bond direction. In the open FS filling range (e.g., ≈ 2.13), it is rotated to be parallel to the shortest bond direction.
The frequencies f 1,2 are a strong constraint on the amount of heterostrain in the TBG sample. Specifically, as illustrated in Fig. 4C, the frequency f 2 is roughly two times f 1 , showing that the two small Fermi pockets have an area ratio ∼ 2 : 1. Theoretically, as illustrated by the solid black lines in Fig. 4C, for a heterostrain strength = 0.2% and = 0 • , the areas A i=1,2 of the two small pockets, when converted to frequency via f −1 i ≡ (Δ 1 B ) i = 2 e ℏA i , are in good agreement with experiment. Furthermore, also note that the frequencies f 1,2 extrapolate to 0 at ≈ ±0.04, showing that the two Dirac points are shifted to finite (opposite) filling fractions by heterostrain.
We observe behavior qualitatively similar in all respects to that in Fig. 4A in all 3 longitudinal contact pairs for which we have dilution-fridge measurements (SI Appendix, Fig. S11).

Summary and Outlook
In summary, we have shown that due to the large size of the moiré unit cell at small twist angles, even a small amount of uniaxial heterostrain on the microscopic scale can lead to dramatic changes in the narrow bands of twisted bilayer graphene. A key feature of the strained bandstructure is the splitting of the respective energetic degeneracies of the two Dirac points and the three van Hove points. The splitting of the two Dirac points leads to a semimetallic state with two small Fermi pockets at the CNP. On the other hand, the two outermost van Hove points bound a broad filling range near = ±2 where the constant energy contours become open. Interestingly, the elongation of the larger Fermi pocket near the CNP is perpendicular to that of the open FSs, the latter being perpendicular to the direction of the shortest moiré triangular bond.
We have analyzed the resulting magnetotransport in strained TBG in the framework of the Boltzmann equation using the method of characteristics, treating the magnetic field nonperturbatively. We showed that a nonsaturating quadratic longitudinal magnetoresistance in a broad filling range near = ±2 naturally arises due to the heterostrain-induced open Fermi surfaces, therefore explaining in remarkable detail experimental results in off-magic-angle devices with lattice anisotropy (45). We have also shown that the sign-changing singularities in the Hall number occur in the open FS filling range and are not directly associated with any van Hove singularity as commonly assumed, e.g., in ref. 68. Furthermore, our theory reveals a dramatic rotation of the transport principal axis as the filling is tuned from the charge neutrality point to the filling range of open Fermi surfaces, without invoking interaction-induced electronic nematicity. Red dashed lines are frequencies from the experimental data. Solid black lines are predictions from the theory for = 0.2% and = 0 • . The two frequencies f 1 and f 2 sum to the one-pocket frequency f 3 that extends beyond the first van Hove point. They additionally account for the nontrivial relative strengths of the quantum oscillations within the bounds of the first van Hove points. As with other details of this work, the theory predicts electron-hole symmetry, while some asymmetry is observed in experiment.
Given the importance of energy-shifted van Hove points in the transport properties of TBG devices, we have analyzed previous quantum oscillation data, which has revealed a Lifshitz transition from two pockets to one pocket at a filling fraction where the innermost van Hove singularity is predicted to occur based on theoretical calculations, offering strong evidence of heterostrain effects on these devices. We have further proposed several additional signatures to look for in future experiments. These include a significant difference in cyclotron mass on either side of the innermost van Hove singularity (probed qualitatively but not quantitatively in the extant experiment) and a principal transport axis with saturating magnetoresistance in the open Fermi surface filling range.
Finally, given the amplifying effect of a small strain at the underlying carbon lattice scale on the moiré lattice scale, the latter of which controls the electronic behavior within the narrow bands, it is tantalizing to consider strain engineering of such devices to achieve effects which in regular solids would require applying strain magnitudes incompatible with structural stability. Data, Materials, and Software Availability. The data for Fig. 4 and SI Appendix, Fig. S11