Pseudo Landau levels, negative strain resistivity, and enhanced thermopower in twisted graphene nanoribbons

As a canonical response to the applied magnetic field, the electronic states of a metal are fundamentally reorganized into Landau levels. In Dirac metals, Landau levels can be expected without magnetic fields, provided that an inhomogeneous strain is applied to spatially modulate electron hoppings in a similar way as the Aharonov-Bohm phase. We here predict that a twisted zigzag nanoribbon of graphene exhibits strain-induced pseudo Landau levels of unexplored but analytically solvable dispersions at low energies. The presence of such dispersive pseudo Landau levels results in a negative strain resistivity characterizing the $(1+1)$-dimensional chiral anomaly if partially filled and can greatly enhance the thermopower when fully filled.

In this paper, we propose a general method based on the band theory to analytically derive the band structure of pLLs induced in a twisted zigzag GNR and then use the resolved pLL dispersions to interpret the transport signatures of the twisted GNR. In Sec. II, we show that the twisted GNR exhibits a bulk zero mode which is the strain-induced zeroth pLL (pLL 0 ) by nature. By linearizing the Hamiltonian of the twisted GNR in the vicinity of the bulk zero mode, i.e., the pLL guiding center, we derive the dispersions of the pLLs at low energies. In Sec. III, we study the low-energy transport of the twisted GNR in the framework of the semiclassical Boltzmann formalism and elucidate that the dispersive pLLs engender a negative strain resistivity if partially filled and can enhance the thermopower if fully filled. Section IV concludes the paper and envisages a few other venues to which our general method may be applied.

II. ELECTRONIC STRUCTURE OF THE TWIST-INDUCED PSEUDO-LANDAU LEVELS
We start from the strain-free tight-binding Hamiltonian of graphene with only nearest neighbor hopping arXiv:2006.09048v2 [cond-mat.mes-hall] 12 Aug 2021 Denoting the honeycomb lattice constant by a, we have δx = √ 3a and δy = 3 2 a. The nearest neighbor vectors (blue) are thus (α1, α2, α3) = ( 1 2 δxx + 1 3 δyŷ, − 1 2 δxx + 1 3 δyŷ, − 2 3 δyŷ), and the next nearest neighbor vectors (red) are (β1, β2, β3) = ( 1 2 δxx + δyŷ, − 1 2 δxx + δyŷ, δxx). In the presence of open boundary conditions along the y direction, the unit cell of the zigzag GNR is chosen as the shaded area (yellow) containing 2N sites on each sublattice. (b) The hexagonal Brillouin zone (green) of the honeycomb lattice in panel (a) contains two Dirac cones at the BZ corners labeled by K and K , respectively. In the presence of spatially uniform anisotropic hoppings (t1, t2, t3) = (t + δt, t + δt, t), where δt < 0, the two Dirac cones are oppositely displaced away from BZ corners as illustrated by the red and blue arrows. (c) A twisted GNR (blue) is obtained by applying a torsional strain that rotates the right (left) edge of the undeformed GNR (light blue) of length L and width W by a small angle of Φ/2 (−Φ/2).
where a r (b r ) is the electronic annihilation operator at site r = (x, y) belonging to the A(B) sublattice of the honeycomb lattice with lattice constant a = 0.142 nm [59][60][61], and α i is the ith nearest neighbor vector along which the hopping is t i independent of r. In the following, unless otherwise specified, we choose the x direction to be parallel to the zigzag edges [ Fig. 1(a)], so that . This tight-binding Hamiltonian encodes two energy bands ε(k) = ±| i t i e ik·αi |, which exhibit Dirac cones at the Brillouin zone (BZ) corners [62] for isotropic hoppings t i=1,2,3 = t. For anisotropic hoppings (t 1 , t 2 , t 3 ) = (t + δt, t + δt, t), the Dirac cones are translated from the BZ corners to k ± W = [± 2 √ 3a arccos(− 1 2 t t+δt ), 0], as illustrated in Fig. 1(b).
An elastic strain can deform the crystalline lattice of graphene, break the translational symmetry, and spatially modulate the hoppings through the empirical for-mula [63] where g = 3.37 is the Grüneisen parameter [63,64], α i (r) is the strain-modulated spacing between a chosen lattice site at r and its ith nearest neighbor, and α i = a is the strain-free counterpart ofα i (r) as illustrated in Fig. 1(a). For the simple twist lattice deformation [ Fig. 1(c)] characterized by the parameter λ = Φ/L that measures the rotational angle of the GNR unit cell [ Fig. 1(a)] per unit length along the x direction, we haveα i (r) ≈ (α 2 i + λ 2 α 2 i,x y 2 ) 1/2 for a sufficiently small twist λa 1 [65]. We note thatα i (r) preserves the modified x-direction translational symmetry Π(δ x ) = T (δ x )R(λδ x ), which should be defined as a regular translation by δ x along the x direction combined with a counter-clockwise rotation by an angle of λδ x around the x axis [ Fig. 1(c)]. Applying Fourier transform in the x direction, we find that the modulation [Eq. (2)] yields for the twisted zigzag GNR a tight-binding Hamiltonian whereŝ −δy b kx,y = b kx,y−δy and t(y) = t exp{g[1 − (1 + 3 4 λ 2 y 2 ) 1/2 ]} corresponding to the hoppings along α 1 and α 2 , while the hopping along α 3 is preserved as t. Therefore, the effect of the twist is similar to that of the aforementioned anisotropy (t 1 = t 2 = t 3 ), relocating the Dirac cones but in a space-dependent fashion.
To scrutinize this relocation, we take the continuum limit such that the shift operator can be estimated through linearization asŝ ±δy ≈ 1 ± δ y d dy , which leads to the Bloch Hamiltonian for the twisted GNR The nanoribbon tight-binding Hamiltonian [Eq. (3)] then becomes H = kx,y ψ † kx,y H kx,y ψ kx,y with Pauli matrices τ x and τ y acting on ψ kx,y = (a kx,y , b kx,y+δy/3 ) T . If t( W 2 ) < 1 2 t, for any given momentum 4π 3δx ≤ k x ≤ 8π 3δx , we can always find within the twisted GNR a pair of spatial coordinates at which the first term in Eq. (4) changes sign. Therefore, for each choice of sign in Eq. (5), there exists a bulk zero mode Ψ 0 (y), which is an even function of y−y 0 , satisfying H kx,y Ψ 0 (y)| y=y0 = 0. Such a bulk zero mode comprises the twist-displaced Dirac points associated with different values of y 0 . To better understand the nature of these bulk zero modes, we investigate the spectrum of the Bloch Hamiltonian [Eq. (4)] at low energies. Because of the exponentially decaying t(y), analytically resolving the eigenvalues of H kx,y is usually not feasible. However, if the twist is sufficiently small, t(y) varies slowly on the lattice scale and can be well estimated in the vicinity of the bulk zero mode through the linearization t(y) ≈ t(y 0 )[1 − Ω y0 (y − y 0 )/t], where Ω y0 = 3 4 λ 2 y 0 gt/[1 + 3 4 λ 2 y 2 0 ] 1/2 . This helps simplify the Bloch Hamiltonian [Eq. (4)] into a minimally coupled Dirac Hamiltonian whose eigenvalues are dispersive pLLs (see Appendix A for detailed derivations) where the momentum dependence is incorporated through Ω y0 by making use of Eq. (5) and mapping its domain of definition from [ 4π 3δx , 8π 3δx ] to [− 2π 3δx , 2π 3δx ], which is located in the first BZ of the twisted GNR. We mention that these pLLs are bounded from above: By noting that the pLL energies cannot exceed the merging points of the two Dirac cones, i.e. the Lifshitz transition points ε ± Lif = ±t, one can estimate from Eq. (7) the maximal pLL index as n max = 2 where · is the floor function; but we emphasize that Eq. (7) fails before n reaches n max because it is based on the linearization of the Bloch Hamiltonian [Eq. (4)]. It is worth noting that the pLLs characterized by Eq. (7) are doubly degenerate due to the multivaluedness of sgn(tδ y Ω y0 ) associated with the contributions from both the upper (y 0 > 0) and the lower (y 0 < 0) sectors of the GNR. The bulk zero modes of the nanoribbon tight-binding Hamiltonian [Eq. (3)] connecting Dirac points K and K are none other than the doubly degenerate pLL 0 and belong to the family of pLLs.
We have applied numerical simulations to the nanoribbon tight-binding Hamiltonian [Eq. (3)] and find that the dispersive pLLs are indeed well captured by Eq. (7) at low energies [ Fig. 2(a)]. We also notice that the pLLs are actually bounded between two Dirac cones in the momentum dimension. This is because when k x approaches the Dirac points, y 0 eventually becomes comparable to the spatial extent of the pLL wave functions; the wave functions then begin to overlap, lose their degeneracy and recombine into dispersive snake states which constitute the Dirac cones [67] [ Fig. 2(b)]. While the energy bands inside the Dirac cones are qualitatively analogous to those in strain-free Dirac cones [ Fig. 2(c)], the strained Dirac cones host far fewer bands, because the strain confines the snake states to the center of the GNR, thereby reducing the effective width of the system.
Besides the pLLs and the bulk states in the Dirac cones of the twisted GNR, we also find a pair of flat bands -π -π/2 0 π/2 π -0 Band structure of the twisted GNR with N = 1200.
(a) The low-energy spectrum (blue) is composed of dispersive pLLs and two Dirac cones from which these pLLs emerge.
The dispersions of these pLLs are accurately captured by Eq. (7) (red), which is overlaid on the energy bands. The twist parameter adopted is λ = 0.0005a −1 , at which the maximal C-C bond elongation appearing at the edges of the GNR is 27%. Such strain, though large, should be sustainable in graphene [66]. (b) A closer look of the left Dirac cone in panel (a). The dashed curve is the envelope of the cone ε ± D (kx) = ±[2 cos( 1 2 kxδx) − 1]t. Each pLL is doubly degenerate, formed by the confluence of two bulk bands at the right sector of the envelope. (c) Dirac cone located at K of an undeformed GNR of the same width harbors many more bands than that in panel (b) in the same energy window.
traversing the whole BZ in Fig. 2(a). We argue that such bands can only be the edge states of the GNR. To justify this argument, we analyze the nanoribbon tight-binding Hamiltonian [Eq. (3)] from another point of view by treating it as a Su-Schrieffer-Heeger (SSH) model [68] with bipartite hoppings 2t(y) cos( 1 2 k x δ x ) and t. Between the Dirac points, the GNR unit cell is divided into three segments by the two domain walls at y 0 such that the outer two segments (i.e., |y| > |y 0 |) are topological SSH chains with 2t(y) cos( 1 2 k x δ x ) < t, while the inner segment (i.e., |y| < |y 0 |) is a trivial SSH chain with 2t(y) cos( 1 2 k x δ x ) > t. Therefore, there are in total four "end states" associated with the outer two segments. Two of them are the doubly degenerate pLL 0 centered at y 0 , while the other two located at ± W 2 constitute the edge states. In the rest of the first BZ, recalling t(y) < t, we have 2t(y) cos( 1 , making the entire GNR unit cell a topological SSH chain possessing a doubly degenerate flat edge state. Therefore, the BZ-wide flat bands are the edge states of the twisted GNR. Before we leave this section, we briefly discuss the effect of electron-electron interactions which we have ignored so far. It is well known that a repulsive Coulomb interaction can create local ferromagnetic order at the two edges of a strain-free zigzag GNR such that the two edges have opposite magnetization [69][70][71][72][73]. The most important consequence on band structure is that a small band gap opens up at the Dirac points [71][72][73], which should in general deform the pLLs in Eq. (7). However, this band gap is inversely proportional to the ribbon width [71], and quickly becomes much smaller than the pLL spacing for realistic interaction strengths and wide ribbons we consider here. (An additional gap arises at k = π δx , but is unlikely to affect the pLLs located within [− 2π 3δx , 2π 3δx ].) Therefore, we would expect Eq. (7) to characterize the pLL dispersions even in the presence of interaction. Similar to the zigzag edge state, the pLL 0 may also carry local ferromagnetic order and form in the twisted zigzag GNR an "edge-compensated" antiferromagnetism [74].

III. TRANSPORT SIGNATURES OF THE TWISTED GRAPHENE NANORIBBON
In Sec. II, we have obtained the pLL dispersions [Eq. (7)] by studying the Bloch Hamiltonian [Eq. (4)] in the vicinity of y 0 , where the bulk zero mode, i.e., the pLL 0 , is centered. In the present section, we adopt these results to study the transport signatures of the twisted GNR at low energies.

A. Density of states
To derive the transport properties of the twisted GNR, it is instructive to first investigate the density of states (DOS) contributed by the doubly degenerate pLLs and the two Dirac cones from which these pLLs emerge. Without loss of generality, we will set the chemical potential µ > 0 in the following, while the µ < 0 case can be treated using the particle-hole transformation.
We first examine the pLLs whose dispersions have been found in Eq. (7). It is critically important to note that these pLLs are only well-defined between the Dirac cones. Therefore, the actual extent of the n-th pLL . The bound k n x can be determined by finding the intersec- Upon finding k n x , the DOS contributed by pLLs is where d + n dkx | µ is calculated in the left BZ here and below. In Eq. (8), we also define for the pLL n the occupancy pa-

Heaviside step function and
] marking the minimum (maximum) of the electron-like pLL n . The values of D n and Γ n for the lowest few pLLs are listed in Table I.
As for the bulk bands harbored by the Dirac cones, we may treat their contribution to the DOS as identical to that of Dirac cones in a strain-free GNR of renormalized width, as explained in Sec. II and Figs. 2(b) and 2(c). In the limit of large GNR width, the DOS associated with the bulk states in the Dirac cones can then be written as where ε(q) = ± 3 2 ta(q 2 x + q 2 y ) 1/2 are the dispersions characterizing the Dirac cones before projection into the onedimensional BZ, and the multiplier ξ(µ) = N λ (µ)/N 0 (µ) concerns the difference in energy band numbers with πt µ] being the number of bands in a single Dirac cone of the twisted (undeformed) GNR intersecting the chemical potential µ > 0. A more substantial derivation of ξ(µ) is provided in Appendix B. The total DOS is then the combination of Eqs. (8) and (9), which is substantiated by comparing to the DOS of the nanoribbon tight-binding Hamiltonian [Eq. (3)] numerically evaluated through the tetrahedron method [75] as illustrated in Fig. 3(a).

B. Negative strain-resistivity
We now present the analysis of the longitudinal electrical conductivity of the twisted GNR. At sufficiently low temperatures, we may apply the Sommerfeld expansion  to the Boltzmann formalism and keep only the lowest order contribution, which, for the pLLs, reads where v x n (k x ) = 1 d + n dkx is the band group velocity for the pLL n and we have used the energy-dependent relaxation time approximation [76] and assumed identical relaxation time τ (µ) for all energy bands. The contribution to the electrical conductivity from the bulk states in the Dirac cones may be written as where v x (q) = 1 ∂ε(q) ∂qx is the band velocity associated with the Dirac cone. Similar to the DOS contribution [Eq. (9)], which is modified from the strain-free Dirac cone contribution by a multiplier ξ(µ), the conductivity contribution [Eq. (12)] also requires a multiplier χ(µ) to incorporate the difference in band numbers. We argue that χ(µ) ≈ 1 2 ξ(µ) because the energy bands in Dirac cones of the twisted GNR are no longer "V-shaped" as those hosted by the strain-free Dirac cones, but are nearly half "V-shaped" as illustrated in Fig. 2(b), making the number of Dirac cones contributing to the electrical conductivity effectively one (see Appendix B).
A comprehensive understanding of the electrical conductivity also requires the knowledge of the relaxation time τ (µ), especially its energy dependence, which is sensitive to the details of scattering mechanism. For strainfree graphene, the screened Coulomb scattering dominates [77,78] and the relaxation time becomes τ (µ) ∝ µ due to the "Diracness" of the charge carriers [79,80].
In the presence of twist, the Dirac cones are broken into pLLs and the µ dependence of the relaxation time would change. Using Fermi's golden rule [81], the scattering time can be written as τ (µ) = C/g(µ), where C = 2π |V if | −2 encodes the detailed information of the scattering mechanism with V if being the scattering matrix element between the initial state |i and the final state |f . When the Drude contribution is the major source of scattering, the parameter C becomes a constant proportional to the density of the impurities [82] consistent with the first-order Born approximation prediction [83]. For other types of scatterers, such as charged impurities, the parameter C is expected to have a smooth dependence on µ. In terms of the parameter C, the total longitudinal electrical conductivity comprising Eqs. (11) and (12)  (13) Assuming Drude scattering, the electrical conductivity σ xx is calculated as a function of µ in Fig. 3(b); our theoretical prediction Eq. (13) is seen to capture the essential features of the numerical results. Though Fig. 3(b) is quantitatively correct only for C = const., the pLLs should generally manifest themselves through singularities in the DOS g(µ) and hence dips in σ xx , as long as C has a smooth µ dependence.
We now analyze the λ dependence of the electrical conductivity [Eq. (13)]. When the twist λ is allowed to vary in a narrow range in which only the dispersive pLL 1 is partially filled with the chemical potential µ 1 = 0.024t [lower orange lines, Figs. 4(a) to 4(d)], the values of the occupancy parameters ν n (µ) and the multipliers ξ(µ) and χ(µ) remain intact. Moreover, the band velocities ∼ d n dkx of the pLLs are increasing functions of λ (see Appendix A). Consequently, both the total electrical conductivity σ xx and the pLL contribution σ xx L increase with the twist λ, giving rise to a negative strainresistivity ρ xx = 1/σ xx [84] as illustrated by the orange curve in Fig. 4(e). Such an effect is analogous to the nega-tive magnetoresistivity [14][15][16][17][18][19][20] in topological semimetals with only chiral LLs (cf. pLL 1 ) partially filled (i.e. in the quantum limit) and may serve as a manifestation of the (1 + 1)-dimensional chiral anomaly [13], which coincides with the valley anomaly in graphene [82]. For a higher chemical potential µ 2 = 0.057t [upper orange line, Figs. 4(a) to 4(d)] outside the quantum limit, pLL 1 (pLL 2,3,4 ) is always fully (partially) filled during the variation of the twist; a negative strain resistivity arising from pLL 2,3,4 then emerges [ Fig. 4(f)]. We mention that the pLL-induced negative strain resistivity can only occur when λ is scanned in a narrow range; otherwise the bottoms of subbands in the Dirac cones [Figs. 4(c) and 4(d)] can in general cross the chemical potential, giving rise to quantum-oscillation-like signals in both the resistivity and the conductivity. Consequently, the negative strain resistivity may become obscured. Lastly, we note that the thermal conductivity of the GNR is related to the electrical conductivity through the Wiedemann-Franz law κ xx = π 2 k 2 B 3e 2 T σ xx . Therefore, it also increases with λ in the quantum limit, giving rise to a negative strain-thermal resistivity similar to that in Weyl superconductors [85].

C. Enhanced thermopower
We now study the thermoelectric effect in the twisted GNR. Because of the preservation of the time-reversal symmetry, the twisted GNR exhibits no Nernst effects or thermal Hall effects in the presence of a temperature gradient ∂ x T . However, the Seebeck coefficient is compatible with the time-reversal symmetry and can be conveniently determined through the Mott relation [86] which arises from the Sommerfeld expansion [76] at low temperatures k B T µ. For the broadening δ = 0.0024t chosen for σ xx in Fig. 3(b), the temperature is k B T ∼ δ . If we are only interested in pLLs (outside the gray patches in Fig. 3) which are far from the charge neutrality point, k B T µ indeed stands and the application of the Mott relation is legitimate. On the other hand, for higher temperatures comparable to the Fermi temperature, a deviation from the Mott relation has been observed [87][88][89]. The thermopower becomes sensitively dependent on the scattering mechanisms [77,78,90,91] and is likely to be dominated by charged impurities [77,78], while electron-electron scattering and electron-optical-phonon scattering may also contribute [90,91].
We note that every time a pLL is fully filled, σ xx is greatly suppressed but experiences a sudden boost giving rise to large values of both 1 σ xx and | dσ xx dµ |, thus significantly enhancing the Seebeck coefficient as illustrated by the sequence of dips around Γ n in Fig. 3(c), which mimic those produced by the ordinary flat Dirac-Landau levels resulting from uniform real magnetic fields  [ 87-89, 92, 93]. The enhanced Seebeck coefficient is accompanied by a boosted Peltier coefficient through the Thomson relation Π xx = T S xx [76].

IV. DISCUSSIONS AND CONCLUSIONS
During our derivation of the pLL dispersions and the associated transport properties, two assumptions have been made. First, we assume the continuum limit and strongly localized pLLs so that both the shift operatorŝ s ±δy and the hopping t(y) can be expanded only to the linear order. We argue that this is a legitimate estimate when the GNR width is not too small. In fact, as we will show in Appendix C, the quadratic and cubic order terms in the expansions ofŝ ±δy and t(y) can only contribute a correction of order of 10 −4 t when N = 1200, which should be safely negligible. Second, we assume that the effect of the next nearest neighbor hopping is negligible. This is because such an effect can be in princi-ple canceled by a fine-tuned electric field [82]. In realistic graphene, the strain-free next nearest neighbor hopping ranges from 0.02t to 0.2t [59] and thus should deform the twist-induced pLLs to some extent in the absence of the applied electric field. The response of pLLs to the next nearest neighbor hopping is discussed in detail in Appendix D.
To summarize, we have studied the strain-induced pLLs in a twisted GNR. By tracking the formation of the bulk zero mode located on the domain wall of the GNR unit cell, which is effectively an SSH model, we reduce the low-energy theory of the twisted GNR into an exactly solvable minimally coupled Dirac theory to derive the momentum dependence of the pLLs. Such dispersive pLLs produce a negative strain resistivity if partially filled and can enhance the thermopower if fully filled.
The method we have adopted deriving the pLL dispersions is by no means ad hoc and should be in principle transplantable to GNRs under other inhomogeneous strain patterns or magnetic fields with complicated or even arbitrary spatial profiles. It can also be conveniently generalized to other graphene-like Dirac materials made of photons [94], magnons [95], phonons [96,97], and Majorana particles [98,99], since strain-induced flat Dirac-Landau levels have been reported in these materials. Our work may also be adapted to superconducting Dirac matter such as Weyl superconductors [85,100,101] and dwave superconductors [102,103], where strain may be the only hope to Landau-quantize the Bogoliubov quasiparticles. In Sec. II of the main text, we have derived the dispersions [Eq. (7)] of the pLLs in a twisted GNR. We here present a detailed solution of the eigenvalue problem of the minimally coupled Dirac Hamiltonian H kx,y [Eq. (6)], from which the pLLs arise. Explicitly, the eigenvalue problem of H kx,y is written as , where f + (y) and f − (y) constitute the eigenvector Ψ(y) = e ikxx [f + (y), f − (y)] T corresponding to the eigenenergy . In a more compact form, Eq. (A1) can be rewritten as with s = ±. For transparency, we introduce the dimensionless coordinate which helps simplify Eq. (A2) as We define an auxiliary function g s (ξ y ) = e ξ 2 y /2 f s (ξ y ) that transforms Eq. (A4) into the exactly solvable Hermite's differential equation (A5) which possesses square-integrable solutions only when the eigenenergy adopts discrete values where n is a nonnegative integer. The solutions of Eq. (A5) are the Hermite polynomials g sgn(tδyΩy 0 ) (ξ y ) = H n (ξ y ) and g −sgn(tδyΩy 0 ) (ξ y ) = sgn(n)H n−1 (ξ y ), making the eigenvectors when sgn(tδ y Ω y0 ) > 0. When sgn(tδ y Ω y0 ) < 0, the eigenvectors should be written as τ x Ψ ± n>0 (y) and τ x Ψ 0 (y). The y-dependence of the eigenvectors is incorporated through Eq. (A3). By plugging these acquired eigenvectors into Eq. (A1), we may specify the sign of the eigenenergy n , which is not encoded in Eq. (A6). Specifically, the eigenenergies associated with the eigenvectors Ψ ± n>0 (y) [τ x Ψ ± n>0 (y)] and Ψ 0 (y) [τ x Ψ 0 (y)] are ± n>0 = ± 2n|tδ y Ω y0 | ( ± n>0 = ∓ 2n|tδ y Ω y0 |) and + 0 = 0 ( − 0 = 0), respectively. The pLL dispersions can then be written as Eq. (7), which possesses two-fold degeneracy.
As analyzed in Sec. III, the dispersive pLLs impact the transport properties, e.g., Eqs. (10) and (13), through the derivative of Eq. (7) calculated at the left BZ. Explicitly, it reads where we have defined the parameter It is straightforward to find that d + n dkx | µ is an increasing function of the twist λ because γ µ is a decreasing function of λ. Therefore, if λ is allowed to vary in a narrow range such that only pLL 1 is partially filled, the electrical conductivity Eq. (13) becomes an increasing function of λ, resulting in a negative strain resistivity.

Appendix B: Derivation of multipliers
The fact that Dirac cones of a twisted GNR [ Fig. 2(b)] contain far fewer bands than those of an undeformed GNR [ Fig. 2(c)] requires multipliers ξ(µ) and χ(µ) to be introduced when calculating the DOS [Eq. (9)] and the electrical conductivity [Eq. (12)], respectively. To figure out the values of these multipliers, we now analyze the bulk bands hosted by the left Dirac cone (i.e., the one located at k D = − 2π 3δx ) of twisted and/or undeformed GNR with close attention paid to the morphology and transport properties associated with these bulk bands.
For the bulk bands inside the Dirac cone of a twisted GNR, as illustrated in Fig. 5(a), when k x increases, an energy band first goes downhill along the left boundary of the Dirac cone and reaches its minimum in the vicinity of the right boundary of the cone. It then climbs uphill toward the right boundary, at which it merges with another energy band and they evolve into a doubly degenerate pLL. The transport associated with this Dirac cone is only contributed by the partially filled energy bands, labeled by ε t n (k x ), which intersect the chemical potential µ > 0. The number of such bands, denoted as N λ (µ), can be estimated by counting the number of pLLs emerging from the Dirac cone below the chemical potential µ. Explicitly, we have . (c) An artificial band structure is constructed by adding to (a) energy bands [red and denoted as ε a n (kx)], which are (almost) parallel to the right boundary of the Dirac cone, making the modified energy bands ε i n (kx) = {ε t n (kx), ε a n (kx)} also V shaped, resembling those in panel (b).
where θ(·) is the Heaviside step function. The first term in Eq. (B1) relates to the electron-like copy of the doubly degenerate pLL 0 . The DOS at chemical potential µ is then where the primed summation only runs over the N λ (µ) partially filled energy bands. It is worth noting that these energy bands are generally not monotonic and each of them may intersect the chemical potential µ at most twice. In that case, we may divide ε t n (k x ) at its minimum into two monotonic segments and let the summation run over both segments. Nevertheless, we conclude from Eq. (B2) that the DOS is dominated by the energy bands with small group velocities. In contrast, the electrical conductivity obtained through the Boltzmann formalism, is mainly contributed by the energy bands with large group velocities. Inside the Dirac cone of an undeformed GNR, we see that most of the energy bands are V shaped and go downhill (uphill) along the left (right) boundary of the cone, leaving the band minima in the center around k D . In addition, there is a state whose energy is a monotonously increasing function of k x . It evolves into the zigzag edge state away from k D on the left [ Fig. 5(b)]. The proximity of band minima to the Dirac point k D indicates a convenient way to estimate the number of partially filled en-ergy bands, denoted as N 0 (µ). Specifically, the nanoribbon tight-binding Hamiltonian [Eq. (3)] is reduced to a Toeplitz tridiagonal matrix at k D , whose eigenvalues are exactly solvable as ε u n (k D ) = −2t cos( nπ 4N +1 ) with n = 1, 2, · · · , 4N . By confining these band minima below the chemical potential µ but above the Dirac point, we obtain a requirement for band index n as 2N 4N +1 π). For µ < −2t cos( 2N +1 4N +1 π), the only partially filled energy band is the one evolving into the zigzag edge state, with index n = 2N + 1. This implies the number of energy bands intersecting the chemical potential µ may be estimated as (B4) We note that N 0 (µ) can be greatly simplified if the chemical potential µ is close to the Dirac point, imposing requirement ε u n (k D ) t, or, equivalently, n − 2N N , on the partially filled energy bands. Such a requirement leads to a uniform subband gap ∆ε u n = ε u n+1 (k D ) − ε u n (k D ) ≈ πt 2N , which gives an alternative estimate N 0 (µ) ≈ µ ∆ε u n = 2N πt µ. In the limit of large GNR width, these N 0 (µ) energy bands contribute to the DOS is the Fermi velocity of the Dirac cone. As for the electrical conductivity, we have which also depends on v F and N 0 (µ).
To figure out the values of the multipliers, a connection between g t D (µ) [σ t D (µ)] and g u D (µ) [σ u D (µ)] is useful. Such a connection can be established by introducing an artificial band structure, in which we require the energy bands ε t n (k x ) to merge with the right boundary of the Dirac cone rather than penetrating the cone and becoming pLLs. Equivalently, each energy band ε t n (k x ) is spliced with an artificial segment ε a n (k x ), whose group velocity is ∼ 3at 2 , as illustrated in Fig. 5(c). The reshaped energy bands, labeled as ε i n (k x ) = {ε t n (k x ), ε a n (k x )} are also V-shaped, resembling those in Fig. 5(b). Therefore, we may treat the artificial band structure in Fig. 5(c) as a Dirac cone in an undeformed GNR whose effective width is smaller than the one considered in Fig. 5(b), because the subband gap is larger. Since the Fermi velocity associated with this Dirac cone is approximately v F = 3at 2 , by comparing to Eqs. (B5) and (B6), the associated DOS g i D (µ) and electrical conductivity σ i D (µ) may be written down as Alternatively, by comparing to Eqs. (B2) and (B3), we write g i D (µ) and σ i D (µ) as where we have noticed that the group velocities (∼ 3at 2 ) of the artificially added energy bands ε a n (k x ) are much larger than those associated with the band bottoms of ε t n (k x ), whose contributions are enclosed in g t D (µ), so that the dominating contribution to g i D (µ) is from g t D (µ). On the other hand, the large band velocities indicate that the artificially added energy bands ε a n (k x ) should have a significant contribution to the electrical conductivity σ i D (µ). The contribution from ε a n (k x ) should be roughly the same as that from the energy bands ε t n (k x ), because the band bottoms of ε t n (k x ) barely affect the transport in the Boltzmann formalism.
We make use of Eqs. (B7) and (B8) and obtain the multipliers as which result in a transport theory [Eqs. (10) and (13)] that can actually capture the most important features of the numerical results.
Appendix C: Higher order terms in the continuum theory The pLL dispersions in the main text [Eq. (7)], while in good agreement with the numerics [ Fig. 2(a)], are obtained by expanding both the modified hopping t(y) and the shift operatorsŝ ±δy to the linear order in the nanoribbon tight-binding Hamiltonian [Eq. (3)]. This approximation becomes accurate under the following conditions: (i) the wave functions vary slowly on the lattice scale δ y , such that the linearization ofŝ ±δy is legitimate; and (ii) the bulk zero mode is strongly localized at the domain wall y 0 such that the widths of the pLL wave functions are sufficiently small, making the linear expansion of t(y) valid. Such a requirement also guarantees minimal overlap between the pLLs and the edge states. In more realistic conditions, we also need to consider in the expansions the higher order terms, which turn out to account for the slight discrepancy between Eq. (7) and numerics [ Fig. 2(a)]. In the rest of this section, we will show that the leading order discrepancy scales as n 3/2 , where n is the pLL index.
Starting from Eq. (3), we now expand both t(y) and s ±δy up to the cubic order and treat the quadratic order term δH D, (2) kx,y and cubic order term δH D, (3) kx,y as perturbations to the Dirac Hamiltonian [Eq. (6)]. Explicitly, these terms are where the k x dependence is acquired from y 0 through Eq. (5). The matrix elements of such perturbations are easily evaluated using ladder operators, and the perturbation theory calculations are lengthy but straightforward. Specifically, the quadratic order terms in Eq. (C1a) contribute a correction , (C3) at the first-order perturbation level. Therefore, one cannot naively expect that the quadratic order terms [Eq. (C1a)] always have a more profound influence than the cubic order terms [Eq. (C1b)]. In fact, the relative importance between Eqs. (C2) and (C3) is sensitive to the value of k x as well as the type of the lattice deformation, though they both scale as n 3/2 . For the twist deformation characterized by Eq. (2), the sizes of these contributions at n = 1 (relative to the pLL 1 dispersion + 1 ) are plotted in Fig. 6(a). Both diverge at the Dirac points, but δ ±,(2) n diverges more quickly; in contrast, δ ±,(3) n dominates around |k x δ x | = π/2, and the two corrections become comparable near the Γ point. Figure 6(b) compares δ +,(2) n + δ +,(3) n with +,exact n − + n , the difference between the numerical pLL dispersions and the Dirac theory prediction [Eq. (7)]. The discrepancy can be traced to subleading corrections arising from even higher-order terms in the expansions of t(y) andŝ δy . Since the linear order expansions ofŝ ±δy and t(y) already capture the numerical energy bands to a high accuracy of 10 −4 t, which can be further increased by including the quadratic and the cubic order terms, we choose not to include more terms in the expansions ofŝ ±δy and t(y).

Appendix D: Next nearest neighbor effects
Our discussion of the twisted GNR in the main text is based on a tight-binding model with only nearest neighbor hoppings, while the next nearest neighbor effect is in general not negligible in realistic graphene. In this section, we investigate the response of pLLs to the next nearest neighbor hoppings.
In the lowest order approximation, we may writê s −δy +ŝ δy ≈ 2 so that only onsite terms appear in Eq. (D3). This observation suggests that the effect of δH may be greatly suppressed if we apply an appropriate electric field E = −∂ y φ(y), which also induces an onsite potential in the GNR Hamiltonian. Since the pLLs sit at the domain wall y 0 , we require the external electric potential energy to cancel the onsite energy resulting from the next nearest neighbor hoppings. Therefore, the electric potential at y 0 is (D4) Making use of Eq. (5), we can remove the k x dependence in Eq. (D4) so that the applied electric potential reads With this electric field, both the electronic structure and the transport of the twisted GNR are governed by the nearest neighbor effect, which is the concern of the main text. Although the next nearest neighbor effect can be suppressed by an external electric field, the realization of such an electric field is unfortunately not easy due to the complicated space dependence associated with the electric potential φ(y) in Eq. (D5). In the following, we show that our prediction on the negative strain resistivity and enhanced thermopower should still be qualitatively correct even without the applied external electric field. By expanding the shift operatorsŝ −δy +ŝ δy ≈ 2+δ 2 y d 2 dy 2 and next nearest neighbor hopping t 1,3 (y) at y = y 0 to the quadratic order of y − y 0 , we can calculate for the pLL dispersions [Eq. (7)] the perturbative correction as which is strongly k x -dependent but only weakly depends on the pLL index n. Although Eq. (D6) is valid for 4π 3δx ≤ k x ≤ 8π 3δx , one can again shift the domain of k x into the first BZ so that the cos( 1 2 k x δ x ) terms change sign. The Band structure and DOS of a twisted GNR with next nearest neighbor hopping. As in the main text, N = 1200, and λ = 0.0005a −1 such that the maximal C-C bond elongation at the edges of the GNR is 27%. The strain-free nextneatest-neighbor hopping is t = 0.1t. (a) Low-energy bands from numerical simulations of the tight-binding model (blue) with the analytical predictions (red) of Eqs. (7) and (D6) for the pLLs and Eq. (D9) for the edge states overlaid. The particle-hole symmetry is broken by the next nearest neighbor hopping, and both the pLL0 and the edge states are now dispersive. (b) The corresponding numerical DOS, broadened by convolution with a Lorentzian of width δε = 0.0024t.
previously flat pLL 0 now becomes dispersive, and the particle-hole symmetry is manifestly broken. A comparison between our theory for the pLLs and numerical spectrum is given in Fig. 7(a). We find Eqs. (7) and (D6) yield reasonably good agreement for the lowlying pLLs comparing to the existing DFT calculations in Ref. [58]. Even better agreement could be achieved by keeping higher-order gradients in Eq. (D3), as we have done in Appendix C. As illustrated in Fig. 7(a), the next nearest neighbor hoppings do not destroy the dispersive pLLs. We thus expect the negative strain resistivity, which results from the intervalley transport mediated by partially filled dispersive pLLs, should persist. We also plot the numerical DOS in Fig. 7(b); as in the t = 0 case, because of the flatness of the pLLs near the Γ point, the DOS is sharply peaked and expected to produce an enhanced thermopower there.
Before we leave this section, we mention that the enhanced thermopower may also arise from the edge states, which are now dispersive and may produce features in the DOS similar to those resulting from the pLLs. It is therefore useful to find the edge state dispersion explicitly. In the framework of the perturbation theory, we first consider the zero mode localized on the A sublattice at the lower edge [see Fig. 1(a)] by writing the Schrödinger equation 2t(y j ) cos( 1 2 k x δ x )φ j + tφ j+1 = 0, where y j = − W 2 + (j − 1)δ y marks the position of the a j site and φ j is the associated probability amplitude. Note the wave function associated with the B sublattice vanishes near the lower edge. Since we have assumed strong deformation at the edge, φ n decays rapidly in the bulk for any k x ; thus we can approximate the slow-varying function t(y) by t(− W 2 ). The normalized wave function then reads where ≡ t(− W 2 )/t < 1/2. Further approximating t 1 (y) and t 3 (y) by their values at y = −W/2, we obtain the edge state dispersion as the expectation value of δH in Eq. (D3), (D9) This is degenerate with another edge state localized at the upper edge. As shown in Fig. 7(a), Eq. (D9) is a good approximation to the numerical edge state dispersion, which again traverses the entire BZ and is now well separated from the Dirac point. The edge band notably produces a shoulder on the pLL peak at ≈ 0.17t.