Analytic solution to pseudo-Landau levels in strongly bent graphene nanoribbons

Nonuniform elastic strain is known to induce pseudo-Landau levels in Dirac materials. But these pseudo-Landau levels are hardly resolvable in an analytic fashion when the strain is strong, because of the emerging complicated space dependence in both the strain-modulated Fermi velocity and the strain-induced pseudomagnetic field. We analytically characterize the solution to the pseudo-Landau levels in experimentally accessible strongly bent graphene nanoribbons, by treating the effects of the nonuniform Fermi velocity and pseudomagnetic field on equal footing. The analytic solution is detectable through the angle-resolved photoemission spectroscopy (ARPES) and allows quantitative comparison between theories and various transport experiments, such as the Shubnikov-de Haas oscillation in the complete absence of magnetic fields and the negative strain-resistivity resulting from the valley anomaly. The analytic solution can be generalized to twisted two-dimensional materials and topological materials and will shed a new light on the related experimental explorations and straintronics applications.

Unfortunately, the pseudo Landau levels induced by the aforementioned strain patterns [33][34][35][36][37][38][39][40] are dispersive and thus are not directly interpretable by the standard Dirac theory established for the ordinary dispersionless Landau levels. For weak strain, the pseudo Landau level dispersions are often overlooked for simplicity until a recent study [40] analytically and nonperturbatively solves such dispersions in a uniaxially strained graphene nanoribbon with a nonuniform Fermi velocity but a uniform pseudomagnetic field. Nevertheless, understanding how pseudo Landau levels disperse in the presence of strong strain is a much more complicated problem remaining largely unexplored. This is presumably * tliu@pks.mpg.de † luhz@sustech.edu.cn because the pseudo Landau levels are expected to occupy a large portion of the Brillouin zone with increased strain; and the standard procedure solving pseudo Landau levels using the linearized Hamiltonians [15, 16, 20-22, 34, 36, 39, 40] at the Brillouin zone corners consequently fails.
In this paper, we present an analytic approach to solve the pseudo Landau levels in bent zigzag graphene nanoribbons under strong strain. In Sec. II, we briefly review two commonly used and analytically solvable Dirac models for weakly bent graphene nanoribbons and demonstrate the applicability as well as the limitations of such models. In Sec. III, we show that the graphene nanoribbon unit cell [ Fig. 1(a)] is effectively a Su-Schrieffer-Heeger model [41] with strain-modulated bipartite hoppings, giving rise to a zero-energy topological domain wall mode [ Fig. 1(b)], which is actually the zeroth pseudo Landau level by nature. Linearizing the lattice model in the vicinity of the domain wall (i.e., the pseudo Landau level guiding center) into an analytically solvable Schrödinger differential equation, we obtain the pseudo Landau level dispersions in a wide range of the Brillouin zone. In Sec. IV, we elucidate that the superiority of the lattice model over the commonly used Dirac models lies in the real-space linearization, which treats the strain-modulated Fermi velocity and the strain-induced pseudomagnetic field on equal footing. In Sec. V, we derive the dispersions of the pseudo Landau levels for more realistic graphene models with the Semenoff mass, the intrinsic spin-orbit coupling, the electric fields, and the next nearest neighbor hoppings. The resolved analytic dispersions enable us to explore, in Sec. VI, the transport resulting from the pseudo Landau levels, exemplified by the Shubnikov-de Haas oscillation in the absence of magnetic fields and the negative strainresistivity arising from the valley anomaly. Section VII concludes the paper and addresses the potential generalization of our real-space approach to a various of Dirac is the x(y) direction spacing between two neighboring sites belonging to the same sublattice. (b) Schematic plot of a fan-shaped graphene nanoribbon circularly bent from a rectangular graphene nanoribbon of length L and width W . Note that the central arc of the bent nanoribbon coincides with the length L of the undeformed nanoribbon, while the width of the bent nanoribbon is identical to its counterpart in the absence of strain. The bend may create in the bulk a domain wall (dashed) at which the bipartite hoppings are identical. The localized domain wall state is the zeroth pseudo Landau level |ψ0(y) by nature. Inset: The circular bend is parameterized by the curvature of the central arc (orange curve), denoted as λ, such that the radius of curvature of the central arc is λ −1 .

II. DIRAC MODELS IN THE WEAK STRAIN LIMIT
We begin by briefly reviewing the commonly used Dirac models of strained graphene [15, 16, 20-22, 34, 36, 39, 40] with a focus on their applicability and limitations. In the framework of nearest neighbor tight-binding theory, the graphene Hamiltonian reads where a R (b R+αi ) annihilates an electron on the A (B) sublattice at position R (R + α i ) with the nearest neighbor vectors (α 1 , α 2 , α 3 ) = (  Fig. 1(a)] measured by the lattice constant a = 1.42Å; and t i is the electron hopping parameter between the site located at R and its ith nearest neighboring site at R + α i . In the absence of strain and anisotropy, the nearest neighbor hopping parameters are set as t i = t = −2.8 eV [42]. External elastic strain alters the positions of lattice sites and thus spatially modulates the hopping parameters. In graphene, such a strain effect is incorporated through the empirical formula where u(r) is the displacement of the lattice site located at position r and g = 3.37 is the Grüneisen parameter [43]. In the weak strain limit, the displacement field u(r) varies slowly on the lattice scale, i.e., |α i ·∇u| |u|. As a common practice [13-15, 21, 29, 34, 39, 40], the empirical formula [Eq. (2)] of the strain-modulated hopping can then be approximated by expanding to the linear order of ∇u as where u µν = 1 2 (∂ µ u ν + ∂ ν u µ ) is the strain tensor. The strain tensor should take its value at the position R+ 1 2 α i such that the hoppings along α i and −α i are the same. For constant strain tensors, the hopping parameters determined by Eq. (3) incorporate no space dependence and the translational symmetry is preserved. By the Fourier transform (a r , b r ) where n uc is the number of unit cells, we obtain the Bloch Hamiltonian where the Pauli matrices σ x,y are defined. According to Eq. (3), t i → t for weak strain. Therefore, the low-energy theory of Eq. (4) can be obtained by linearizing H k in the vicinity of the Brillouin zone corners k η = (η 4π where η = ±1 is the valley index; (v η x , v η y ) = 3ta 2 (−η, 1) is the Fermi velocity; and q = k − k η is measured from the corners of the Brillouin zone. Since h η q is in a Peierls substitution form, we can define a strain-induced vector potential A η = η g 2ae (u yy − u xx , 2u xy ). Though h η q is obtained by assuming constant strain, we argue that it is in fact a legitimate theory even if the strain tensor incorporates space dependence, because the strain only varies slowly. In particular, for the weak circular bend, the displacement field reads u = λ(xy, − 1 2 x 2 ) [29], where λ is the curvature of the central arc of the bent nanoribbon [Inset, Fig. 1 where a strain-induced uniform pseudomagnetic field can be defined as B η = ∇ × A η = η g 2ea λẑ. The spectrum of h η q comprises of the dispersionless Dirac-Landau levels where the integer n is the Landau level index. Equation (7) is often referred to as the pseudo Landau levels in order to be distinguished from those Landau levels produced by ordinary magnetic fields. Unfortunately, Eq. (7) is only capable of capturing the numerical band structure, which is obtained by diagonalizing H [Eq. (1)] under the strain modulation t 1,2 = t(1− 3 4 gλy) and t 3 = t [Eq. (3)], right at the projected Brillouin zone corners 2(a) and 2(b)], because h η q merely encloses the terms linear in the momentum q and the strain tensor u xx = λy (note u xy = u yy = 0). To improve the match, we include additional higher order terms O(q y q x ), O(u xx q x ), and O(u xx q y q x ) to h η q and obtain a modified Dirac theory where the renormalized Fermi velocities areṽ η , but the pseudomagnetic field is intact [cf., Eq. (6)] to the lowest order of y. The diagonalization of h η q is analogous to the Sturm-Liouville problem analyzed in Ref. [40] and the spectrum of h η q can be analytically solved as which indeed better fits the numerical band structure in the vicinity of the projected Brillouin zone corners [Figs. 2(c) and 2(d)]. It is worth noting that E η n (q x ) only captures the bulk bands bounded between the two projected Dirac cones DC max = ± ṽ η x (q x − η g 2a λy)| y=±W/2 , while the dispersive energy bands inside the projected Dirac cones and the flat energy bands emerging from the projected Dirac points are clearly originated from the marginal regions as reflected by the average of the position operatorȳ = dy ψ * nkx (y) y ψ nkx (y), where ψ nkx (y) is the wave function, as illustrated in Figs. 2(c) and 2(d). The real-space position of these energy bands can also be resolved by the spectral function, which is detailed in Appendix A.
We mention that E η n (q x ) derived from the modified Dirac Hamiltonian h η q can gradually lose its validity when the bend curvature λ is increased. In fact, the acquisition of h η q relies on two important approximations: (i) A momentum space expansion (with respect to q) of the Bloch Hamiltonian [Eq. (4)] in the vicinity of the Brillouin zone corners. (ii) A real space linearization (with respect to y) of the exponentially varying strain-modulated hopping [Eq. (2)]. However, a strong strain inevitably extends the pseudo Landau levels in the momentum space and renders the momentum space expansion around the Brillouin zone corners inadequate. Moreover, the overlooked higher order terms by the linearization in the real space can become more important at strong strain. Consequently, a more sophisticated theory valid for strong strain would be desired and worthy of investigation.

III. LATTICE MODEL IN THE STRONG STRAIN LIMIT
In Sec. II, we have seen that the Dirac models are only applicable in the weak strain limit, but the expansion  around the Brillouin zone corners would lose its ground for strong strain and the additional higher order terms can transform the low-energy theories to non-Dirac models, where neither the Fermi velocity nor the pseudomagnetic field can be well defined. In the present section, we develop a real-space approach based on the band topology analysis to derive the dispersions of the pseudo Landau levels induced by strong (as well as weak) circular bend.
For the circular bend lattice deformation [ Fig. 1(b)], the length of the central arc coincides with the nanoribbon length before bending and the width of nanoribbon is unchanged. This implies that the azimuthal projection of a chemical bond alters linearly with the y coordinate, while the radial projection of the bond is unchanged. Specifically, along the bonds α 1,2 , the projections in the azimuthal direction becomex·α 1,2 (1+λy). But the bond α 3 remains intact. According to the empirical formula [Eq. (2)], the modulated hopping parameters are which preserve the x direction translational symmetry.
We are thus able to perform the partial Fourier transform (a r , b r ) T = N −1/2 uc kx e ikxx (a kx,y , b kx,y ) T , where N uc is the number of unit cells of the bent graphene nanoribbon [ Fig. 1(b)], to obtain a tight-binding Hamiltonian for the circularly bent graphene nanoribbon as [2t(y) cos( 1 2 k x δ x ) +tŝ δy ]a kx,y− δy 6 + H.c., (11) where δ x = √ 3a, δ y = 3 2 a, andŝ δy is a shift operator satisfyingŝ δy a kx,y = a kx,y+δy . At a given momentum k x , the nanoribbon tight-binding Hamiltonian [Eq. (11)] becomes a Su-Schrieffer-Heeger model [41] with intracell hopping 2t(y) cos( 1 2 k x δ x ) and intercell hopping t. Due to the y dependence of the hopping parameters [Eq. (10)], for momenta |k x | ≤ k c = 2 δx arccos( 1 2 e −g/2 ), a domain wall can possibly appear at where the two hoppings are the same, while no domain wall can exist if k c < |k x | ≤ π δx , in which case the intercell hopping is always overwhelmed.
The position of the domain wall has a profound influence on the band topology of the nanoribbon tightbinding Hamiltonian [Eq. (11)]. For an undeformed nanoribbon with λ = 0, the domain wall can only be located within the nanoribbon at the k x = ±k D . For |k x | > k D (|k x | < k D ), the intercell (intracell) hopping dominates and the unit cell becomes a topological (trivial) Su-Schrieffer-Heeger chain with (without) a pair of end modes. It is such end modes that constitute for the momenta k D ≤ |k x | ≤ π δx the well-known flat zigzag edge states [ Fig. 3(a)]. For a moderately bent graphene nanoribbon with 0 < λ < λ c , where For a given momentum k x is this range, the upper (lower) sector of the unit cell is topological (trivial), giving rise to an end mode and a domain wall mode at the charge neutrality point [ Fig. 3(b)]. The end modes at all allowed momenta, i.e., k − max ≤ |k x | ≤ k + max , constitute a dispersionless energy band located at the stretched zigzag edge, while the domain wall modes result in a flat bulk band, which must be interpreted as the zeroth pseudo Landau level, since no other bulk states are expected to be dispersionless. For the momenta |k x | > k + max (|k x | < k − max ), the unit cell realizes a purely topological (trivial) Su-Schrieffer-Heeger model [ Fig. 3(b)]. Therefore, a pair of flat edge states composed of Su-Schrieffer-Heeger end modes are expected at k + max < |k x | ≤ π δx , which corresponds to the momentum-space scope of the edge state located at the compressed edge. As for the stretched edge, the ranges of the edge state add up to k − max < |k x | ≤ π δx . For a critically bent nanoribbon with λ = λ c , the pseudo Landau levels from the left half and the right half of the which is related to the nanoribbon tight-binding Hamiltonian [Eq. (11)] through H = kx,y ψ † kx,y H kx,y ψ kx,y with the sublattice basis ψ kx,y = (a kx,y−δy/6 , b kx,y+δy/6 ) T . Note that we have taken the continuum limit in Eq. (11) such that the shift operator is written asŝ δy ≈ 1 + δ y d dy .
Because of the complicated space dependence of t(y), analytically solving the Schrödinger differential equation characterized by H kx,y is generally not feasible. But the band topology analysis has revealed the nature of the zeroth pseudo Landau level being the Su-Schrieffer-Heeger domain wall mode, and thus locates the common guiding center of all pseudo Landau levels in the real space, provided that there are no electric fields or next nearest neighbor hoppings, whose effects are detailed in Secs. V C and V D. Since the pseudo Landau levels are well localized states, their dispersions can be in principle accurately approximated by studying the nanoribbon Bloch Hamiltonian [Eq. (13)] in the vicinity of their common guiding center. We find it more convenient to work with the momenta k x ∈ [ π δx , 3π δx ] and then maps the resolved dispersions of the pseudo Landau levels back to the conventional first Brillouin zone. Such a manipulation introduces no artifacts because the legitimate energy bands must have a 2π δx period in k x , even though the nanoribbon Bloch Hamiltonian [Eq. (13)] seemingly has a 4π δx period due to the specific form of the Fourier transform we have chosen. For the momenta k x ∈ [ π δx , 3π δx ], the position of the domain wall should be rewritten as which can be reduced to Eq. (12) by setting k x → k x + 2π δx . In the vicinity of the domain wall, i.e., the common guiding center, the nanoribbon Bloch Hamiltonian [Eq. (13)] is restored to a standard Dirac Hamiltonian where Alternatively, such a Dirac Hamiltonian may be written as a matrix operator where B = 2|Ω 0 tδ y | is the energy scale. In Eq. (16), a andâ † are the ladder operators defined aŝ in which we have defined the dimensionless parameter ξ a = (y − 0 )/l B with l B = |tδ y /Ω 0 | being the magnetic length. To solve the spectrum of h kx,y , we adopt the trial solution |ψ n>0 = (ζ A,n |n , ζ B,n |n − 1 ) T and |ψ 0 = (ζ A,0 |0 , 0) T , where |n is defined to be an eigenstate of the bosonic number operatorâ †â , satisfyingâ †â |n = n |n . Explicitly, |n can be written as where H n (·) is the nth Hermite polynomial. It is straightforward to verify that |ψ n>0 (|ψ 0 ) is the eigenvector of h kx,y when ζ 2 A,n = ζ 2 B,n (ζ A,0 = 0). We here choose ζ A,n = ∓1/ √ 2, ζ B,n = 1/ √ 2, and ζ A,0 = 1. And the explicit eigenvectors are which correspond to the spectra n>0 = ± B √ n and 0 = 0, respectively. Mapping back to the first Brillouin zone through k x → k x + 2π δx , we obtain the explicit dispersions of the pseudo Landau levels (19) is our key result, whose validity is justified by the good match in a wide range of momenta to the numerical band structure resulting from directly diagonalizing the nanoribbon tight-binding Hamiltonian [Eq. (11)] for a maximally bent graphene nanoribbon [ Fig. 4(a)]. It is also worth noting that the derivation of n (k x ) does not depend on the specific value of the bend curvature λ. Therefore, Eq. (19) is in fact applicable for both strong and weak strain. Consistent with our aforementioned analysis, Eq. (19) is defined for k x ∈ [−k c , k c ], in which the domain wall l 0 can possibly exist, while the range of the pseudo Landau levels cannot exceed the subset [−k + max , k + max ] in order to confine the domain wall l 0 inside the nanoribbon. Comparing to Eqs. (7) and (9) derived from Dirac models [Eqs. (6) and (8)] in the weak strain limit, Eq. (19) is equally accurate at the projected Brillouin zone corners k x = ±k D but exhibits much lower discrepancy with respect to the numerical band structure elsewhere for |k

IV. SUPERIORITY OVER DIRAC MODELS IN THE WEAK STRAIN LIMIT
In Sec. III, we have elucidated that the dispersive pseudo Landau levels [Eq. (19)] are more accurate than those [Eqs. (7) and (9)] arising from the Dirac models [Eqs. (6) and (8)] in the strong strain limit. Such a finding may not be surprising because the Dirac models are only applicable in the weak strain limit. We are thus motivated to examine whether the superiority of Eq. (19) can retain in the weak strain limit. According to Sec. II, the modified Dirac model [Eq. (8)] is a more accurate low-energy theory for weak strain. We thus focus on the comparison between Eqs. (19) and (9) in the present section.
We intuitively expect Eqs. (19) and (9) to have similar performance in fitting the numerical band structure in the vicinity of the projected Brillouin zone corners k x = ±k D for weak strain. This is because the hopping modulation [Eq. (10)], which is the ground for Eq. (19), can be reduced in the weak strain limit to t 1,2 = (1 − 3 4 λgy) and t 3 = t, identical to the condition [i.e., Eq. (3) with the displacement field u = λ(xy, − 1 2 x 2 )] we use to derive Eq. (9). Surprisingly, we find Eq. (19) exhibits much smaller deviation to the numerics than Eq. (9) even for weak strain [Figs. 5(a) and 5(b)]. Since the only difference between the two analytic dispersions n (k x ) and E η n (q x ) lies in the hopping modulation, we thus attribute the difference to the higher order terms [e.g., O(λ 2 y 2 )] overlooked during the linearization of the  (19), (9), and (7)] are plotted as solid, dashed, and dot dashed curves, respectively. Left (right) inset enlarges the energy differences associated with Eqs. (19) and (7) [Eq (19) and (9)] in the vicinity of kx = −kD (dotted line).

strain-modulated hopping t(y).
To substantiate this claim, we rewrite the modified Dirac Hamiltonian [Eq. (8)] as with the nonuniform velocity parameters where we temporarily do not specify the space dependence of the strain-modulated hopping t(y); and the strain-induced vector potential A η x = 2 η 3ea t(y)−t t(y) gives rise to a pseudomagnetic field  Because of the simultaneous spatial inhomogeneity in the velocity parameters and the pseudomagnetic field, the Schrödinger differential equation associated with h η q is generally not analytically solvable except for t(y) with simple (e.g., linear) space dependence.
For the purpose of deriving the spectrum of h η q , we shall follow the strategy established in Sec. III by studying h η q in the vicinity of the pseudo Landau level guiding center y 0 , which coincides with the domain wall 0 when t(y) adopts the form of Eq. (10). For a strain-modulated hopping t(y) of generic space dependence, according to Ref. [44], the guiding center of the pseudo Landau levels is determined by A η x (y 0 ) = − q x /e such that there exists a zero-energy mode in the spectrum of h η q to be interpreted as the zeroth pseudo Landau level. By expanding in the vicinity of the guiding center y 0 , it is straightforward to find whose spectrum is completely determined by the velocity parameters and the pseudomagnetic field at the guiding center y 0 . Making use of the condition A η which are independent of the specific space dependence of t(y). However, the pseudomagnetic field sensitively depends on the form of t(y) due to the appearance of ∂ y t(y). Explicitly, it reads where the coefficient reads f qx = 1 + 3 2 ηaq x for the linearized hopping modulation [Eq. (3)] and f qx = 1 + 1 2g ηaq x for the full empirical hopping modulation [Eq. (10)]. The resulting pseudo Landau level dispersions are where the former is simply the slightly dispersive pseudo Landau levels E η n (q x ) in Eq. (9); and the latter corresponds to n (k x ) in Eq. (19) expanded in the vicinity of the projected Brillouin zone corners ηk D . Indeed, the latter is much less dispersive than the former by a ratio of 3g, which confirms our observation in Figs. 5(a) and 5(b).
The finding that the higher order terms overlooked during the linearization of the strain-modulated hopping do affect the pseudomagnetic field [Eq. (25)] but do not impact the Fermi velocity [Eq. (24)] up to the linear order of q x suggests that the widely used strainmodulated hoppings with linear space dependence [13-15, 21, 29, 34, 39, 40] may be insufficient in characterizing the dispersions of the strain-induced pseudo Landau levels. To find the accurate dispersions, one would need to adopt the full space dependence of the hopping parameters without any approximation. But the complicated space dependence of such hopping parameters may hardly result in analytically solvable Schrödinger differential equations, which govern the dispersions of the pseudo Landau levels. In contrast, our analytic method is rooted in the band topology analysis; does not rely on the specific form of the space dependence of the hopping parameters; and thus can be transplanted to strain patterns beyond circular bend as long as such strain patterns are still characterized by t 1,2 → t(y) and t 3 → t.

V. DISPERSIONS OF PSEUDO LANDAU LEVELS IN REALISTIC GRAPHENE
In Sec. III, we derive the dispersions of the pseudo Landau levels using a simple nearest neighbor tight-binding model [Eq. (11)] of a bent graphene nanoribbon. In realistic graphene samples, there are several inevitable effects: (i) the Semenoff mass arising from the interplay with the substrate; (ii) the Haldane mass due to the intrinsic spin-orbit coupling; (iii) the electric fields; (iv) the next nearest neighbor hoppings. The deformation of pseudo Landau levels in the presence of such effects are respectively analyzed in this section.

A. Semenoff Mass
The interplay between the graphene and the substrate where it is hosted breaks the chiral symmetry by introducing a staggered potential characterized by a Semenoff mass [45]. The magnitude of the Semenoff mass m S closely relies on the details of the substrates. For hexagonal boron nitride (hBN) substrates [46], density functional calculations reveal m S = 27 meV, while m S in silicon carbide (SiC) [47,48] can be as large as m S = 135 meV. Due to the presence of the Semenoff mass, the linearized Bloch Hamiltonian [Eq. (15)] acquires an extra term and becomes which may be rewritten in terms of the ladder operators [Eq. (17)] as With the trial solution |ψ I n>0 = (ζ I A,n |n , ζ I B,n |n − 1 ) T and |ψ I 0 = (ζ I A,0 |0 , 0) T , we find that h I kx,y can be diagonalized when the parameters adopt the following val- Note that the zeroth pseudo Landau level is no longer located at the charge neutrality point but is pushed to m S in the energy dimension [Figs. 6(a) and 6(b)]. Analysis of y reveals that the two segments of the zeroth pseudo Landau level are still connected by the edge state originating from the compressed edge, which has the same sublattice support, while the other edge state located on the stretched edge and originally degenerate with the zeroth pseudo Landau level in the absence of m S is now separated from the zeroth pseudo Landau level by a band gap of 2m S .

B. Spin-orbit coupling
The chiral symmetry can also be broken intrinsically by the spin-orbit coupling H SO ∼ s · (∇V × k), where s is the Pauli matrix in spin space [49,50]. Such a spin-orbit coupling term further breaks the time-reversal symmetry and is known to topologically gap out the Dirac cones of graphene by introducing a Haldane mass [51], which possesses opposite signs at the different projected Brillouin zone corners. The effect of the spin-orbit coupling can be modeled by the following imaginary next nearest ✏ neighbor hopping terms where r a (r b = r a + α 1 ) labels the lattice sites belonging to the A (B) sublattice; and (β 1 , β 2 , β 3 ) = ( ) are the next nearest neighbor vectors [red arrows, Fig. 1(a)]; and d i measures the strength of the spin-orbit coupling associated with β i in the presence of the circular bend. For simplicity, we assume d i to be exponentially varying, similar to the modulation of the nearest neighbor hoppings [Eq. (10)]. Explicitly, d i reads where d measures the spin-orbit coupling without strain. By applying the partial Fourier transform in the x direction, the nanoribbon tight-binding Hamiltonian [Eq. (11)] should be supplemented by where, for transparency, we have defined the parameter D kx,y = 2 sin( 1 2 k x δ x )[d 1 (y + 1 2 δ y )ŝ δy + d 1 (y − 1 2 δ y )ŝ −δy ] − 2d 3 (y) sin(k x δ x ) and set y a = y − 1 6 δ y and y b = y + 1 6 δ y such that we may write Eq. (32) in the sublattice basis ψ kx,y = (a kx,y−δy/6 , b kx,y+δy/6 ) T as H SOa + H SO b = kx,y ψ † kx,y H SO ψ kx,y with the correction to the nanoribbon Bloch Hamiltonian [Eq. (13)] being a purely diagonal matrix H SO = diag(D kx,y−δy/6 , −D kx,y+δy/6 ). For experimentally available bend with λa 1, it is straightforward to see from Eq. (31) that all d i are slowly varying on the lattice scale such that H SO can be estimated through linearization as where the first chiral symmetry breaking term is associated with the Haldane mass and opens up a band gap; and the second term emerges from the small separation of sublattices in the y direction and shifts the energy bands in a y dependent fashion. Although the parameter D kx,y explicitly encloses shift operatorsŝ ±δy , it can be approximated as a purely scalar function of y where we work in the continuum limitŝ ±δy ≈ 1 ± δ y d dy ; take the linearization of d 1 (y ± 1 2 δ y ); and only keep the lowest order terms. Since we are only interested in the low-energy pseudo Landau levels, which are localized around the domain wall 0 , it would be sufficient to study D kx,y exactly at this domain wall. The resulting momentum dependent D kx, 0 acts as the Haldane mass m H ≡ D kx, 0 . In such an approximation, the linearized Bloch Hamiltonian [Eq. (15)] should be rewritten as where we have neglected the second term in Eq. (33), because such a term only contributes at 0 a tiny shift to the pseudo Landau levels. The spectrum of h II kx,y can be directly written down by comparing to Eq. (29) as

C. Electric field
In the presence of a uniform electric field E = Eŷ along the y direction, each of the electrons on the lattice acquires a potential energy −eφ(y), where the electric potential is chosen as φ(y) = −Ey − φ 0 . The linearized Bloch Hamiltonian [Eq. (15)] is then rewritten as h III kx,y = h kx,y + eφ 0 σ 0 + eEyσ 0 − 1 6 eEδ y σ z , which has chiral symmetry preserving onsite terms e(φ 0 + Ey)σ 0 and a mass term − 1 6 eEδ y σ z due to the small separation of sublattices along the direction of the applied electric field. We write h III kx,y in a matrix form as where we define the parameters E = eEl B / √ 2 and m = − 1 6 eEδ y for transparency. To solve the eigenvalues of h III kx,y , we construct the following relation where III n is the eigenvalue of h III kx,y with respect to the eigenvector |ψ III n and we have defined the auxiliary ma- kx,y −m) 2 with no ladder operators in its off-diagonal entries. The dispersions of the pseudo Landau levels can then be obtained by resolving the eigenvalues of K. To diagonalize K, we apply a reversible (but not unitary) transformation to the eigenvector |ψ III n = P |ψ III where we have defined the parameter ω = 2 B − 4 2 E . After the transformation, Eq. (39) can be rewritten as where P −1 KP is a purely diagonal matrix operator and reads We now remove the terms linear inâ andâ † by translationâ where the shifted ladder operators arê with the dimensionless parameter ξ b = 2 √ 2 III n E /ω 2 + ξ a . In terms of these shifted ladder operators, Eq. (42) becomes We then remove the pairing ladder operators (i.e.,b †b † andbb) through the Bogoliubov transformation where the rotated ladder operators arê with the dimensionless parameter In terms of the these rotated ladder operators, Eq. (45) becomes We plug Eq. (48) into Eq. (41) and solve the dispersions of the pseudo Landau levels to be  Such pseudo Landau levels are symmetric with respect to the Brillouin zone center because of the time-reversal symmetry, and thus are fundamentally different from the ordinary Landau levels that produce quantum Hall effects [52,53].
We now briefly mention the effects of electric fields in the other two directions. A z direction electric field breaks the mirror symmetry and brings up an extrinsic spin-orbit coupling Rashba term H R ∼ẑ·(s×k) [54]. The Rashba spin-orbit coupling arising from experimentally available electric fields is typically small comparing to the nearest neighbor hoppings [55][56][57][58][59], and thus should not drastically alter the strain-induced pseudo Landau levels in principle. On the other hand, an x direction electric field can drive a current of electrons along the nanoribbon and lead to longitudinal transport, which will be detailed in Sec. VI C.

D. Next nearest neighbor hopping
In realistic graphene, electrons can also hop to the next nearest neighboring sites belonging to the same sublattice   and produce in the tight-binding Hamiltonian additional terms Unlike the Hamiltonian [Eq. (30)] used to model the spinorbit coupling, the hopping parameters in Eq. (51) are chosen to be purely real and exponentially varying as where t ∈ [0.02t, 0.2t] is the next nearest neighbor hopping in the absence of strain [60]. Following the procedure we have formulated in Sec. V B, it is straightforward to find out that the nanoribbon Bloch Hamiltonian [Eq. (13)] now approximately acquires an extra term where we have made use of the fact that t i are slowly varying on the lattice scale when λa 1 and defined parameter T kx,y = 2 cos( 1 2 k x δ x )[t 1 (y + 1 2 δ y )ŝ δy + t 1 (y − 1 2 δ y )ŝ −δy ]+2t 3 (y) cos(k x δ x ). We notice that H kx,y resembles the terms [cf. Eq. (37)] induced by an electric field E = Eŷ with T kx,y playing the role of the potential energy −eφ(y) = eφ 0 + eEy. Although the parameter T kx,y contains shift operatorsŝ ±δy and thus is different from the electrostatic energy, it is straightforward to show that such a parameter is approximately a purely scalar function of y as where the continuum limit ofŝ ±δy and linearization of t 1 (y ± 1 2 δ y ) are taken. For our purpose of finding the dispersions of low-energy pseudo Landau levels, it would be sufficient to study T kx,y in the vicinity of the domain wall 0 through the linearization T kx,y = T kx, 0 +T 0 (y − 0 ), where the derivative T 0 = dT kx ,y Consequently, the linearized Bloch Hamiltonian [Eq. (15)] should be rewritten as h IV kx,y = h kx,y + T kx, 0 σ 0 + T 0 (y − 0 )σ 0 − 1 6 δ y T 0 σ z , (55) which is analogous to Eq. (37) with T 0 in place of the force eE. By comparing to Eq. (50), we can immediately write down the pseudo Landau levels which well match the numerically calculated band structure [Figs. 7(c) and 7(d)]. Because of the similarity between the parameter T kx,y and the potential energy −eφ(y), the next nearest neighbor effect can be exactly cancelled at (and greatly suppressed around) the projected Brillouin zone corners by an electric potential

VI. TRANSPORT OF BENT GRAPHENE NANORIBBONS
We have performed a systematic study on the analytic dispersions of pseudo Landau levels in bent graphene nanoribbons in Secs. II-V. To allow comparison to experiments, analytic evaluations of transport signatures of bent graphene nanoribbons would be greatly favored. In the present section, we first justify the sufficiency of our nearest neighbor lattice model of bent graphene nanoribbons. We then phenomenologically find the analytic dispersions of the marginal energy bands spliced to the pseudo Landau levels. Ultimately, the transport signatures including the density of states (DOS), the longitudinal electrical conductivity, and the Seebeck coefficient are analytically evaluated and compared to their numerical counterparts.

A. Justification of the nearest neighbor lattice model of bent graphene nanoribbons
In Sec. V, we have elucidated that the pseudo Landau levels resulting from the nearest neighbor lattice model [Eq. (11)] are vulnerable to a variety of mechanisms such as the Semenoff mass, the spin-orbit coupling, the electric fields, and the next nearest neighbor hoppings. Despite appearing unavoidable at the first sight, these effects can actually be neglected in certain conditions. Specifically, we require a pristine graphene sample prepared on a proper substrate (e.g., hBN [46] would be superior over SiC [47,48]), where the Semenoff mass arising from the interplay with the sample is minimized; the Haldane mass (Rashba effect) resulting from the intrinsic (extrinsic) spin-orbit coupling is proved to be much smaller than the nearest neighbor hopping and can thus be neglected [49,50,[55][56][57][58][59]; and the ubiquitous next nearest neighbor hoppings can be compensated by a properly tuned uniform y direction electric field as discussed in Sec. V D. Under such conditions, it would be sufficient for the lattice model to only enclose the dominant nearest neighbor hopping terms.
Our nearest neighbor lattice model, for simplicity, only encloses the in-plane circular bend, which inhomogeneously stretches (compresses) the upper (lower) half of the nanoribbon [ Fig. 1(b)], while ignores the potential out-of-plane strain effects as a common practice [34,35,37]. In fact, the compressive strain, even as weak as 0.1% [61,62], can induce out-of-plane lattice deformation (e.g., bubbles and/or wrinkles) [61][62][63][64], which can further complicate the strain-modulated hoppings [Eq. (10)] by breaking the x direction translational invariance. To suppress such compression-induced buckling, graphene samples should be rigidly attached to the substrate or tightly sandwiched by two substrates such that the out-of-plane lattice deformation is constrained [61]. To avoid the buckling in the experimental implementation, a circular bend created only by tensile strain is preferred. Such a bend can still be modeled by Eq. (10) but the domain of definition of the coordinate should be adjusted to y ∈ [0, W ] from y ∈ [− W 2 , W 2 ]. Such a shift is analogous to a gauge transformation, which only relocates the guiding center but does not affect the dispersions of the pseudo Landau levels. Therefore, even for the more experimentally accessible bent graphene nanoribbons created by pure tensile strain, our key result [Eq. (19)] arising from the nearest neighbor lattice model [Eq. (11)] can still characterize the pseudo Landau levels.

B. Phenomenological analytics of marginal energy bands
A full analytic analysis of the transport of bent graphene nanoribbons requires the knowledge of all energy bands. The pseudo Landau levels [Eq. (19)] are clearly the bulk bands of the bent graphene nanoribbon because their common guiding center is constrained in the bulk through − W 2 ≤ l 0 ≤ W 2 . However, such pseudo Landau levels are distributed around l 0 with a characteristic width ∼ l B as reflected by their wave functions [Eq. (18)]. Therefore, when the guiding center approaches to the edges (i.e., within a few l B 's), the pseudo Landau levels begin to be affected by the edges and evolve into more dispersive energy bands in the marginal regions of the nanoribbon, consistent with our observation onȳ in Figs. 2(a)-2(d). Such energy bands are thus referred to as the "marginal energy bands" in order to be distinguished from the dispersionless topological edge bands. To investigate the transport of the bent graphene nanoribbon, we aspire to quantify such marginal energy bands on a phenomenological basis. For transparency, we only consider the energy bands in the left half of the Brillouin zone, while the energy bands in the right half can be obtained by time-reversal operation.
We first note that the width of a pseudo Landau level in the momentum space decreases with an increased Landau level index n as illustrated in Fig. 8(a). In fact, a higher pseudo Landau level has a more extensive wave function because of more nodes in the Hermite polynomial H n (·) in Eq. (18), making it easier to touch the edges of the nanoribbon and thus more confined in the momentum space due to the monotonic dependence of l 0 on k x [cf., Figs. 3(b)-3(d) and Eq. (12)]. Consequently, all the pseudo Landau levels are phenomenologically bounded between two projected Dirac cones [ Fig. 8 which are the projected spectra of the Bloch Hamiltonian H k with the parameters in Eq. (4) set to t 1,2 = t(± W 2 ), t 3 = t, and k y = ± 2π 3a as well as the strong strain counterparts of DC max . Denoting the ends of the nth pseudo Landau level as k l,r n , whose values are determined by finding the crossings of the projected Dirac cones [Eq. (57)] with the pseudo Landau levels [Eq. (19)], we find the real-space range of the nth pseudo Landau levels to be [l 0 (k l n ), l 0 (k r n )]. At l 0 (k l,r n ), the pseudo Landau levels begin to evolve into marginal energy bands, whose dispersions are governed bỹ which is obtained by linearizing the nanoribbon Bloch Hamiltonian [Eq. (13)] around l 0 (k l,r n ) = 0 (k l,r n + 2π δx ). The first two terms of Eq. (58) turn out to be a Dirac Hamiltonian [cf., Eq. (15)] characterizing pseudo Landau levels centered at l 0 (k l,r n ), while the last term can be understood as a shift to the guiding center. However, it is critically important to note that such a term must not be absorbed into the Dirac Hamiltonian, because the absorption would relocate the guiding center to somewhere outside the allowed scope of the nth pseudo Landau level. In the vicinity of k l,r n , where the linearization [Eq. (58)] of the nanoribbon Bloch Hamiltonian is legitimate, the last term in Eq. (58) can be treated as a perturbation. Performing the perturbation calculations forh 2 kx,y , we find the first order correction to the eigenvalue to be t 2 [1 − cos( 1 2 k x δ x )/ cos( 1 2 k l,r n δ x )] 2 , which is doubly degenerate due to the particle-hole symmetry. Therefore, the analytic dispersions of the marginal energy bands are We note that Eqs. (19) and (59) together with the flat topological edge bands at the charge neutrality point constitute an artificial band structure [ Fig. 8(b)] that satisfactorily mimics the numerical band structure [ Fig. 8(a)].
We thus expect such a band structure can phenomenologically capture the transport associated with the numerical energy bands.

C. Transport signatures
The phenomenologically derived artificial band structure allows analytic investigation of transport signatures and comparison to numerics as well as experimental observations. Without loss of generality, we here only consider the transport of electron-like energy bands (i.e., µ > 0) and conduct explicit calculations in the left half of the Brillouin zone (i.e., k x < 0), while the transport associated with the hole-like energy bands (i.e., µ < 0) and the right half of the Brillouin zone (i.e., k x > 0) can be found using the particle-hole symmetry and the time-reversal symmetry, respectively.
We first consider the DOS of the bent graphene nanoribbon. In the bulk, the energy bands are dispersive pseudo Landau levels [Eq. (19)]. The corresponding   19)]. The orange (green) scatters mark the connection of Eqs. (19) and (59a) [Eqs. (19) and (59b)] at the boundary of the projected Dirac cones (dotted). Quantum oscillations of (c) the DOS, (d) the electrical conductivity, (e) the Seebeck coefficient at a fixed bend curvature λ = 0.642 µm −1 . Quantum oscillations of (f) the DOS, (g) the electrical conductivity, (h) the Seebeck coefficient at a fixed chemical potential µ = 0.112 eV. In panels (c)-(h), the blue curves represent the quantities calculated from the numerical band structure in panel (a) using the tetrahedron method [65]; the red curves represent the quantities calculated from the analytic band structure in panel (b) using Eqs. (62), (65), and (66); the units g(µ, 0) = 4W µ 9πa 2 t 2 and σ xx (µ, 0) = 9e 2 8 2 Ca 2 t 2 are respectively the DOS and the electrical conductivity in the absence of strain; the parameter L = 2.45 × 10 −5 V/K 2 is closely related to the Lorenz number [66]; T is the temperature; and all the data are broadened by convolving in energy a Lorentzian of width δ = 5.6 meV to simulate the effects of disorder and finite temperature.

bulk DOS reads
where we define for the nth pseudo Landau level the occupancy parameter ν n (µ, λ) = θ[ n (k r n )−µ]−θ[ n (k r n )−µ] with θ(·) being the Heaviside function. The dependence on the bend curvature λ in ν n (µ, λ) is acquired from k l,r n . In the marginal regions, the energy bands are characterized by Eq. (59), whose contribution to the DOS reads The resulting total DOS of the bent graphene nanoribbon thus reads where the doubling is to include the contribution from the right half of the Brillouin zone. We find the calculated total DOS [Eq. (62)] satisfactorily fits the DOS numerically evaluated through the tetrahedron method [65] as illustrated in Figs. 8(c) and 8(f) for a bent graphene nanoribbon with varying chemical potential and bend curvature, respectively. Such good matches substantiate our claim on the dispersions of the marginal energy bands [Eq. (59)].
We then turn to calculate the longitudinal electrical conductivity of the bent graphene nanoribbon by the Boltzmann equation approach [66] at low temperatures (i.e., k B T t √ gλa). The bulk conductivity contributed by the pseudo Landau levels reads where, through change of variables, we can rewrite the relaxation time as τ a n (k x , λ) = τ a n ( a n , λ) with a n being the dispersion of the nth energy band in the artificial band structure [ Fig. 8(b)]. We further assume, for simplicity, an identical relaxation time τ a n (µ, λ) = τ (µ, λ) in the second line of Eq. (63). In the framework of the Fermi's golden rule, the relaxation time is inversely proportional to the DOS as τ (µ, λ) = C/g(µ, λ), where the proportionality coefficient C encodes the information of the scattering potential in the bent graphene nanoribbon.
It is worth noting that for a certain bend curvature λ that makes the nth pseudo Landau level partially occupied, i.e., ν n (µ, λ) = 1, the marginal energy bands always have little influence on the relaxation time because the DOS is mostly contributed by the nth pseudo Landau level. The bulk conductivity [Eq. (63)] is then reduced to σ xx b (µ, λ) = e 2 C 2 2 ( d n dkx ) 2 µ , which turns out to be a decreasing function of 1/λ [67] and implies a negative strain-resistivity analogous to the negative magnetoresistivity in the chiral magnetic effect of Weyl semimetals [4][5][6][7][8][9][10][11]. This negative strain-resistivity is originated from the dispersive pseudo Landau levels [Eq. (19)], which play the same role as the chiral zeroth Landau levels in Weyl semimetals. Moreover, it reflects the non-conservation of the valley charge η, i.e., the valley anomaly [40], which is a direct manifestation of the (1 + 1)-dimensional chiral anomaly [68].
Despite the interesting anomalous transport in the bulk conductivity σ xx b (µ, λ), the major source of contribution to the total longitudinal electrical conductivity is actually from the marginal regions as The total longitudinal electrical conductivity then reads With the µ dependence of σ xx (µ, λ) figured out, it is straightforward to calculate the Seebeck coefficient S xx (µ, λ) through the Mott relation [69] S xx (µ, λ) = − which is plotted in Figs. 8(e) and 8(h). The oscillatory behavior of the Seebeck coefficient is inherited from the longitudinal electrical conductivity [Eq. (65)].

VII. CONCLUSIONS
In conclusion, the dispersions and the transport of the pseudo Landau levels in a strongly bent graphene nanoribbon are analytically studied. Such a study is motivated by the fact that the widely used Dirac models [15, 16, 20-22, 34, 36, 39, 40] workable for comparatively weak strain become insufficient in the strong strain limit due to the oversimplification ignoring the nonlinear terms of the momentum and the strain tensor. Applying the band topology analysis based on the hidden chiral symmetry [70], we find that the unit cell of a bent graphene nanoribbon effectively maps to a Su-Schrieffer-Heeger model [41] with strain-modulated bipartite hoppings. A domain wall separating the topological and trivial sectors of the unit cell results from the strain modulation and carries a zero-energy mode, which is the zeroth pseudo Landau level by nature. In the vicinity of such a domain wall (i.e., the guiding center of the pseudo Landau levels), we restore the Schrödinger differential equation into an analytically solvable standard Dirac equation through linearizing the model Hamiltonian. In contrast to the standard linearization adopted when deriving the Dirac models around the Brillouin zone corners, our linear expansion is conducted in real space. It thus treats the strain-modulated Fermi velocity and the strain-induced pseudomagnetic field on equal footing to give an analytic solution to the pseudo Landau levels. The resolved pseudo Landau level dispersions are accurate in a wide range of the Brillouin zone for strong strain and are even superior over the Dirac models for weak strain.
Having acquired the dispersions of pseudo Landau levels using a nearest neighbor lattice model of bent graphene nanoribbons, we turn to consider more realistic models with chiral symmetry breaking masses, applied electric fields, and next-nearest-neighbor hoppings. The Semenoff (Haldane) mass arises from the interplay with the substrate (the intrinsic spin-orbit coupling) and opens up a trivial (topological) band gap to pseudo Landau levels. Nevertheless, comparing to the nearest neighbor hopping effect, the effect of the mass terms is comparatively small in graphene. On the other hand, the electric fields and the next nearest neighbor hoppings can be strong perturbations and thus drastically affect the electronic structure by suppressing and tilting the pseudo Landau levels. Fortunately, these two effects can cancel each other and the resulting bulk bands show no obvious difference from the pseudo Landau levels derived from the nearest neighbor lattice model. The analytically derived pseudo Landau levels and the phenomenologically approximated marginal energy bands constitute an artificial band structure allowing the analytic computation of the transport signatures (e.g., Shubnikov-de Haas oscillation in the absence of magnetic fields and the negative strain-resistivity resulting from the valley anomaly) and the comparison to numerics and experimental observations.
In Sec. II of the main text, we plot the spectrum [Figs. 2(a)-2(d)] of a weakly bent graphene nanoribbon and use the average value of the position operator [i.e., y = dy ψ * nkx (y) y ψ nkx (y), where ψ nkx (y) is the wave function] to mark the positions of the energy bands. We here show that such positions can also be resolved by the spectral function.
The spectral function of a bent graphene nanoribbon can be written as which represents the local density of states (LDOS) at the nth site with H nm (k x ) being the Hamiltonian matrix of the bent graphene nanoribbon. Equation (A1) allows us to study the spectral density in any part of the nanoribbon. For example, we may define the spectral function of the stretched (compressed) marginal region [the uppermost (lowermost) 5% of the bent graphene nanoribbon] by summing A n ( , k x ) over the sites belonging to that part of the nanoribbon. We first calculate the spectral function in the bulk of the bent graphene nanoribbon at low energies [Figs. 9(a) and 9(b)] and find a set of slightly dispersive bulk bands, which are the strain-induced pseudo Landau levels, bounded between the projected Dirac cones [white curves in Figs. 9(a) and 9(b)] characterized by DC max = ± ṽ η x (q x − η g 2a λy)| y=±W/2 [i.e., the weak strain limit of Eq. (57)]. On a phenomenological basis, the effect of the bend on a certain Dirac point is to relocate its position along the x direction in a y dependent fashion. The trace constituted by the displaced Dirac points at different y ∈ [− W 2 , W 2 ] is then a flat band spreading in the vicinity of the chosen Dirac point; and thus is the zeroth pseudo Landau level by nature. For a higher pseudo Landau level, the wave function has more nodes and consequently gets less localized in the real space. Its width in the momentum space becomes narrower. Eventually, all the pseudo Landau levels are bounded between the aforementioned two projected Dirac cones. We also notice that the energy bands outside the bounded area unambiguously belong to the stretched strongly dispersive and spliced to the pseudo Landau levels at the boundaries of the projected Dirac cones except for a pair of longer flat bands [Figs. 9(c) and 9(d)] degenerate with the zeroth pseudo Landau level and a pair of shorter flat bands [ Fig. 9(e) and 9(f)] connecting the two sectors of the zeroth pseudo Landau level across the Brillouin zone boundary. Further calculations of the spectral functions on the zigzag edges, i.e., the a 1 site and the b 2N site in Fig. 1(a), clarify that the longer (shorter) flat bands emerging from the projected Dirac points at ±k − max (±k + max ) are located on the stretched (compressed) edge as illustrated in Fig. 10(a) [Fig. 10(b)]. We thus refer to such flat bands as the edge states, which have a topological origin (cf., Sec III), while call those dispersive bands inside the projected Dirac cones the marginal energy bands.