Hofstadter butterfly and Floquet topological insulators in minimally twisted bilayer graphene

We theoretically study Hofstadter butterfly of a triangular network model in minimally twisted bilayer graphene. The band structure manifests periodicity in energy, mimicking that of Floquet systems. The butterfly diagrams provide fingerprints of the model parameters and reveal the hidden band topology. In a strong magnetic field, we establish that mTBLG realizes Floquet topological insulators (FTIs) carrying zero Chern number, while hosting chiral edge states in bulk gaps. We identify the FTIs by analyzing the nontrivial spectral flow in the Hofstadter butterfly, and by explicitly computing the chiral edge states. Our theory paves the way for an effective practical realization of FTIs in equilibrium solid state systems.

Introduction-Orientation misalignment in twisted bilayer graphene gives rise to a moiré superlattice greatly changing the electronic band structure. When the twist angle θ is around the magic angle (∼ 1.1 • ), the lowenergy moiré bands become nearly flat [1], and strongly enhanced many-body interactions induce superconducting and correlated insulating states [2][3][4][5][6][7][8][9][10]. Away from the magic angle, the system can behave quite differently. The focus of this work is on minimally twisted bilayer graphene (mTBLG) with θ 1 • . Although the flatband physics may not be significant here, mTBLG is a highly interesting system with intriguing properties worthy of serious experimental attention.
The mTBLG hosts a triangular network of onedimensional (1D) conducting channels when a large electric field is externally applied in theẑ direction (out of plane) [11][12][13][14][15][16][17][18][19][20][21][22]. Low-energy electronic states in AB and BA regions are gapped out by the electric field, but domain walls separating the two regions support gapless 1D channels, which are valley Hall kink states with valleydependent chiralities [23]. The AA regions act as junctions that connect domain-wall states along different directions [ Fig. 1(b)]. The low-energy electronic structure in this system can be captured by a triangular network model [15]. While energetics of the network model has been studied, a full characterization of the system, including its topological properties, remain an interesting open question.
In this Letter, we study the electronic structure of mT-BLG in a magnetic field Bẑ by calculating the Hofstadter butterfly of the network model. In the low B field regime, we show that Dirac points can reemerge at finite fields, which should provide fingerprints for determining the network model parameters. In the strong B field regime, the probabilities of an electron to deflect right and left at a scattering junction can be significantly different, which leads to gap opening at Dirac points. In such a case, we demonstrate that the system effectively realizes Floquet topological insulators (FTIs) [24,25], where Chern numbers of bulk bands are zero, but there are chiral edge states traversing bulk gaps. We identify the FTIs based on a nontrivial spectral flow in the Hofstadter butterfly as well as an explicit calculation of the chiral edge states. We also discuss experimental implications of our work, predicting the exciting possibility of realizing a solid state equilibrium FTI in mTBLG.
Model.-In the network model of mTBLG [ Fig. 1], each domain-wall link that connects adjacent junctions has a length d (i.e., the moiré periodicity) and hosts two 1D channels per spin and per valley. We neglect the spin index thereafter, and treat the electronic structures in the two valleys K and K separately since the system is smooth at the atomic scale. Noting that the two valleys are related by C 2z symmetry (twofold rotation aroundẑ axis), we mainly focus on states in K valley. We further assume that the two 1D channels in K valley are decoupled at the domain wall and at the junctions, which leads to a single-channel network model [15]. Our theory can be generalized to a two-channel network model, and we expect that our qualitative predictions still apply.
In the single-channel model, an electron moves chirally along domain-wall links but experiences scattering at junctions [ Fig. 1 represent the incoming and outgoing wavefunctions respectively, andŜ is a 3 × 3 unitary scattering matrix. The form ofŜ is dictated by symmetries, in particular, C 3z (threefold rotation aroundẑ) and C 2z T (C 2z combined with spinless time-reversal symmetry T ). Under the C 3z condition,Ŝ can be modeled as followŝ To satisfy the unitarityŜ †Ŝ = I , we use the follow- Here P aā is the probability for forward scattering, while P ab and P ac are those for deflected scatterings. In addition to P aā + P ab + P ac = 1, these probabilities are further constrained by the requirements that the phases Φ and Φ are real valued [ Fig. 1(c)]. The symmetry C 2z T requires thatŜ =Ŝ , and therefore, P ab = P ac . TheŜ matrix in Eq. (1) is equivalent to that in Ref. [15] when P ab = P ac , but describes a generalization allowing a broader range of parameters including P ab = P ac .
Band structure.-The network model in Fig. 1 describes a periodic array of junctions in a triangular lattice with lattice vectors R 1 = ( √ 3d/2, d/2) and R 2 = (− √ 3d/2, d/2). The electronic structure can be computed by combining the Bloch's theorem and the scattering matrixŜ as follows whereT = diag e ik·(R1+R2) , e −ik·R2 , e −ik·R1 is the translation operator with k being the Bloch wave vector, e i Ed v describes the phase accumulation in traveling along a link, E is the energy, and v is the domain-wall state velocity.
The energy spectrum E(k) is then simply given by the phase of the eigenvalues ofT −1Ŝ . Thus, E(k) has an infinitely repeated pattern with a period 2π v/d [15] due to the resemblance of Eq. (3) to the Floquet problem [26][27][28]. The occurrence of infinitely many bands is a consequence of the continuum approximation in each 1D domain-wall state. In realistic mTBLG systems, we expect a large but finite number of repeated bands within the energy gap opened in AB/BA regions.
We plot the band structures in Fig. 2 for representative values of P aā . In the presence of C 3z and C 2z T symmetries, the band structures show generic features that are independent of P aā value: (1) there are three bands within one energy period 2π v/d; (2) each pair of adjacent bands are connected via a Dirac cone located at one of the three high-symmetry momenta (Γ,K,K ) in the moiré Brillouin zone; (3) Dirac velocity is always v/2, precisely half of the domain-wall state velocity [15]. However, details of the band structures do vary with P aā . Particularly, particle-hole asymmetry of the Dirac cones becomes stronger as P aā approaches 1 (decoupled 1D chiral chains limit) [29]. Hofstadter butterfly and Dirac points.-In the presence of a magnetic field Bẑ, the phase accumulation of an electron traveling along a link becomes exp(i Ed v − i e A·dr), where −e < 0 is the electron charge and A is the vector potential. In the Landau gauge A = (0, Bx), the system is translationally invariant inŷ direction but generically not alongx. We focus on the commensurate cases such that eB √ 3d 2 /h = p/q where p and q are integers. The flux per triangle area ( Under such a flux, we choose a magnetic unit cell with lattice vectors R 1 = (q √ 3d, 0) and R 2 = (0, d). We calculate the Hofstadter butterfly (i.e., the energy spectra as a function of φ ) by choosing q to be prime numbers ≤ 37 and p = 0, 1, 2, .., 4q − 1 [29].
Another effect of the B field is to modify theŜ matrix, since C 2z T symmetry is explicitly broken by the B field. Therefore, the ratio γ ≡ P ab /P ac can deviate from 1. We first consider the case where the B field is small and γ can be approximated by its zero-magnetic-field value (i.e., γ = 1). In Fig. 3, we plot the Hofstadter butterfly for φ ∈ [0, 1] in one energy period (2π v/d) for different values of P aā with γ = 1. Independent of P aā , the electronic structure develops fractal energy gaps for generic φ except for φ = 1/2, where the bands are similar to that at zero field and also feature Dirac cones with velocity v/2. In addition, the energy spectra of 1 − φ and φ are identical given that γ = 1.
Details of the Hofstadter butterfly depend on the value of P aā . We first focus on φ 1 where Landau levels (LLs) emanating from zero-field Dirac cones can be identified. In particular, for P aā = 1/9, the energy of LLs can be accurately fitted by E = ±(v/2) √ 2 eBn for n = 0, 1, 2..., as expected from the LL spectrum of 2D massless Dirac fermions. For other P aā values, the LLs are complicated because of particle-hole asymmetry and other electron pockets coexisting with the Dirac cones in the zero-field band structure [ Fig. 2]. It is worth mentioning that the spectrum has no φ dependence when P aā = 1 (decoupled 1D chiral chains limit). Going beyond the low field limit, Dirac cones reemerge in the spectra at certain values of φ . In Fig. 4(a), we plot the energy levels versus the inverse flux 1/φ . The occurrence of the Dirac nodes sensitively depends on the value of P aā as we show in Fig. 4(b), which could be used to estimate the P aā value from transport experiments.
Floquet topological insulators.-We now turn to the case where the B field is strong enough such that γ has a noticeable deviation from 1. For definiteness, we consider a flux φ ∈ [N, N +1], where N represents a large integer. Within this field range, we assume that P aā and γ = 1 do not vary with B. The corresponding Hofstadter butterfly is shown in Fig. 5 for representative values of (P aā , γ).
Dirac cones at φ = N are gapped out because γ = 1, as demonstrated in Fig. 5. We focus on the parameter space where this Dirac gap opening leads to a full spectral gap, i.e., the blue regions in Fig. 1(d). Physically, an observable gap needs to exceed thermal or disorder broadening-explaining why a strong magnetic field is required in practice.
While the spectra are identical at φ = N and N + 1 under the approximation that γ is fixed, there is a nontrivial spectral flow from φ = N to N +1, which signals a nontrivial band topology [30]. In particular, the n 0 th gap at φ = N evolves to the n 1 th gap at φ = N + 1, where n 1 is n 0 − 1 and n 0 + 1, respectively, for γ > 1 and γ < 1. We note that the difference ν ≡ n 1 − n 0 is well defined, although n 0 and n 1 individually are not because the spectrum is unbounded in our continuum approximation. The spectral flow can be obtained analytically in the P ab = 1 (γ → ∞) or the P ac = 1 (γ = 0) limit. For P ab = 1 (P ac = 1), electrons move clockwise (counterclockwise) in localized chiral loops encircling isolated triangles [ Fig. 6(c)], and the entire spectrum is linearly shifted by −1/3 (1/3) of 2π v/d with the increase of one unit flux per triangle [29], which precisely leads to ν = −1 (ν = 1). Here ν serves as a topological invariant that counts n edge , the number of chiral edge modes, as proved in Ref. [30]. Remarkably, the invariant ν takes the same value for every gap at φ = N . This implies that the Chern number for each band at φ = N is zero, which we confirm through an explicit calculation as described in the Supplemental Material [29].
We establish the existence of chiral edge states by calculating the band structure on a stripe with an open (periodic) boundary condition alongx (ŷ) direction as illustrated in Fig. 6(c). The single-valley network model does not have an analog of a trivial insulator because it is based on 1D chiral states, "anomalous" by construction. Therefore, we take into account both valleys in the stripe-geometry calculation. We assume the two valleys, which are related by C 2z symmetry, to be decoupled in the bulk, but consider perfect inter-valley 180 • reflections along edges where an incoming K-valley (K -valley) state is back scattered at a junction [black dots in Fig. 6(c)] to an outgoing K -valley (K-valley) state along the same domain-wall link. The boundary then separates the bulk from a trivial insulator [29]. In Fig. 6(b), we plot the energy spectrum of the stripe and confirm the existence of chiral edge states. Both the number and chirality of edge states are the same for every bulk gap, which agrees with the bulk topological invariant ν. The existence of the edge state can be understood by the skipping orbits in the P ab = 1 or P ac = 1 limit, as shown in Fig. 6(c). The chirality of edge states depends on sgn(γ − 1), which is also consistent with the bulk butterfly diagram. Therefore, our system realizes FTI, of which n edge is determined by the Floquet winding number instead of Chern numbers [24][25][26][27][28]. Here we infer the topological invariant from the spectral flow of the Hofstadter butterfly, which provides a convenient way to diagnose the band topology. The invariant ν is also confirmed by explicitly computing the edge states. Because of the periodic energy spectrum, we predict a cascade of chiral edge states with the same chirality and mode number prevailing in bulk gaps.
Experimental signatures -Our predictions can be investigated by transport measurements in mTBLG. In mT-BLG with θ = 0.1 • , the moiré period is d ≈ 140 nm and the required magnetic field to generate one flux quantum per triangle is B 0 ≈ 0.49T. It clearly suggests that the mTBLG is more advantageous in studying Hofstadter butterfly than graphene on hBN [31,32] because the required magnetic field is even lower.
Our theory agrees with several existing experiments [13,14], which all demonstrate that the domain-wall network dominates the transport at low temperature and for low density. For example, the magnetoresistance shows oscillations in B which is consistent with the approximate periodicity in the magnetic flux [14]. The sign alternation in the Hall resistivity implies the existence of Dirac points and van Hove singularities in the network model [14,15]. The Dirac nodes can be identified by combining both the Hall resistivity (sign change) and the longitudinal resistivity (peaks). As we show in Fig. 4, the positions of Dirac nodes provide fingerprints to extract the parameters of the junction. Therefore, our theory establishes a bridge between theory and experiment by studying the position of Dirac nodes in a weak magnetic field. For the network model in a sufficiently large magnetic field so that (C 2z T breaking) bulk gaps exceed disorder/thermal broadening, we predict an effective realization of the FTI in an equilibrium system. In this topological state, chiral edge states with the same mode number and chirality can appear repeatedly in a cascade of bulk energy gaps. Experimental observation of a such phenomenon would provide smoking-gun evidence for FTI. We note that similar network models could also be realized in other Dirac systems with a periodic array of mass domains, for example, using the surface states of three-dimensional topological insulators.
We expect that the FTIs in our network model are robust against certain amount of disorder [33,34]. The mTBLG features a new type of disorder, i.e., the twist angle disorder [35], whose effect on the network model is an interesting future direction to pursue. In addition, many-body interaction can drive phase transitions [36][37][38]. Electron repulsion can generate the FTIs in the absence of Bẑ field by spontaneously breaking the C 2z T symmetry, and may possibly lead to fractionalized states in the network model.

CONSTRAINTS INŜ MATRIX
The scatteringŜ matrix with C 3z is characterized by P aā and γ = P ab /P ac as described in the main text. To ensure the unitary condition, two phase parameters, Φ and Φ , are required. These phases are defined modulo 2π. Changing Φ and Φ to −Φ and −Φ also satisfies the unitary condition, but the obtained energy spectrum changes by an overall minus sign. There are in total 6 ways to constructŜ for a given (P aā , γ). We make a particular choice (as described in the main text) without loss of generality. The findings and the conclusions in this work do not depend on this specific choice.

BAND STRUCTURE AND DENSITY OF STATES
We extend the discussion on band structures of the network model used in our theory. We first focus on a few special limits as follows: 1. P aā = 1: This limit realizes a system with decoupled chiral chains along three orientations of the domain-wall links,ŷ, 2x − 1 2ŷ . The band structure is simply a linear superposition of three chiral dispersions and is insensitive to an external magnetic field.
2. P ab = 1 (P ac = 1): The electrons move clockwise (counterclockwise) around localized isolated triangular loops for P ab = 1 (P ac = 1). The band structure features three flat bands in an energy period 2π v/d. The energies of those flat bands are E = ± 1 3 π v/d, π v/d within the first energy period for an integer magnetic flux per triangle.
3. P aā = 1/9 and P ab = P ac = 4/9: The band structure contains nearly particle-hole symmetric Dirac nodes. The low energy properties near the Dirac points agree well with the Dirac equation, e.g., the Landau level spectrum discussed in the main text of our paper. Besides the special limits, we plot the band structures along with density of states in Fig. S1 for a few representative parameters. Close to P aā = 1, there are electron pockets coexisting in energy with the Dirac cones, and the density of states is finite almost everywhere. Close to P aā ≈ 1/9, a clear Dirac nodal structure in the density of states is manifest. In the absence of C 2z T (γ = 1), direct band gaps open at the Dirac nodes, but the system might not develop a full band gap due to the complications of electron pockets.

NETWORK MODEL UNDER AN OUT-OF-PLANE MAGNETIC FIELD
We discuss the procedures for calculating the energy spectrum in the presence of an out-of-plane magnetic field. We first choose the Landau gauge, A = (0, Bx), such that the translational symmetry is generically broken along thex direction. The commensurate limit (eB √ 3d 2 /h = p/q with p,q being integers) is considered here. We choose a magnetic unit cell with and lattice vectors R 1 = (q √ 3d, 0) and R 2 = (0, d) as illustrated in Fig. S2. The unit cell can be split into q subcells that each contains 6 domain-wall links. The phase accumulations due to the Aharonov-Bohm effect (i.e., − e A · dr) are given by: where φ a,m is the ath link in the mth subcell. The flux pattern in the network model is presented in Fig. S2(b). Each clockwise (counterclockwise) loop encircles a flux φ = p 4q (−φ = − p 4q ). With the magnetic unit cell introduced in Fig. S2(a), the spectrum can be straightforwardly calculated by imposing the Bloch equation. We generate butterfly diagrams by choosing all the prime numbers q ≤ 37 and p = 1, 2, ..., 4q − 1. We consider a sufficiently fine grid in the Bloch wavevector space such that all the spectra converge. Practically, we need at least N x = 4 and N y = 12 where N x and N y are the number of wavevectors in k x and k y directions respectively. Now we discuss the spectral flows in the butterfly diagrams, which reveal the topological properties of the band. We plot the special limits, P ab = 1 and P ac = 1 in Fig. S3. The spectral shift under the increase of a unit flux, δφ = 1, is −1/3 (1/3) of 2π v/d for P ab = 1 (P ac = 1). Such energy shifts lead to ν = −1 and ν = 1 for P ab = 1 and P ac = 1 respectively. The amount of the spectral shift can be understood analytically in the following. For P ab = 1, each electron move clockwise and encloses an isolated triangle. Each state thus contains three domain-wall links, and the magnetic flux encircled is exactly φ [ Fig. S2(b)]. Since each link picks up exp(i Ed v − i e A · dr), we obtain 3 δEd v + 2πδφ = 0 where δE is the energy shift due to the change of φ . Similarly, we derive 3 δEd v − 2πδφ = 0 for P ac = 1. Under the increase of one flux per triangle, the n 0 th gap flows to the n 1 the gap, where ν ≡ n 1 − n 0 is well defined. Electrons that follow this flow have a density δρ = ν/A MUC , where A MUC = 2A is the area of one moiré unit cell (MUC) that encloses two triangles. Therefore, when the magnetic field increases by one flux per triangle, the number of electrons per triangle that follow the spectral flow is ν/2. This might indicate that ν/2 is the topological invariant. However, the above analysis only takes into account the K valley, whereas the full system includes both the K and K valleys. Because the two valleys are related by C 2z symmetry, they have the same spectrum pattern with identical spectral flow. It is crucial that the two valleys must be considered together in order to get the correct topological invariant. The reason is that the network model in each valley is based on 1D chiral states, and is therefore "anomalous" by construction. The topological invariant for the full system with both K and K valleys is given by ν, which is also consistent with our edge state calculation.

BERRY CURVATURES AND CHERN NUMBERS
We present an explicit calculation of Berry curvatures and Chern numbers for the network model. In the absence of an out-of-plane magnetic field, the C 2z T symmetry enforces the Berry curvature to be identically zero. A finite Berry curvature is allowed once the C 2z T symmetry is broken (i.e., γ = 1). For simplicity, we consider a magnetic field that generates an integer flux φ = N , where N ≥ 1 is an integer. The Aharonov-Bohm phases generated by this integer flux can be completely gauged away, and the moiré translational symmetry can be restored. In this case, the band energy and wave function are fully determined by the scatteringŜ matrix with γ = 1: whereT = diag e ik·(R1+R2) , e −ik·R2 , e −ik·R1 is the translation operator with k being the Bloch wave vector, and R 1,2 are lattice vectors for the moiré lattice. E k,n and |u k,n are, respectively, the energy and wave function of the nth band at momentum k. We calculate the Berry curverature Ω by numerically computing the Wilson loop. We discretize the momentum space by a rhombus grid, and each small grid spans a momentum-space area A 0 = A MBZ /N 2 , where A MBZ is the momentum-space area of the moiré Brillouin zone (MBZ) and N is a sufficiently large integer characterizing the linear dimension of the grid. The Berry curvature is approximated by ,n |u k2,n u k2,n |u k3,n u k3,n |u k4,n u k4,n |u k1,n ] /A 0 , where k 1 → k 2 → k 3 → k 4 → k 1 tracks in a counterclockwise manner a small rhombus grid with the area A 0 . The Chern number C of the nth band can then be calculated following the definition We show the calculated Berry curvature for a representative band structure of the network model with γ = 1 in Fig. S4. For each band in Fig. S4(a), the Berry curvature Ω reaches its maximum (positive) and minimum (negative) values around momentum points with gapped Dirac cones. The hot spots of Ω correspond to gap opening regions, as expected. The Chern number for each band shown in Fig. S4(a) is identically zero, because the two hot spots of Ω in each band take opposite values and cancel out.
Therefore, we confirm through explicit calculation that the bands of the network model at an integer flux have a zero Chern number.

EDGE STATE AND BOUNDARY CONDITION
FIG. S5. Network model with edges. We consider a strip with a finite width alongx and a periodic boundary alongŷ. η1, η2, η3, and η4 are the reflection phases on the edges.
We provide detailed procedures for calculating the edge states in the network model. The main obstacle is that the single-channel network model with only the K valley is composed of the chiral states. One cannot construct a topologically trivial state from the chiral states because the chiral states are always anomalous. The resolution is to introduce both the K and K valleys such that there are two counter-propagating modes in each domain-wall link, i.e., non-chiral state. We will discuss the construction in depth in the rest of the section.
First of all, we digress and discuss how to construct a trivial state from a model including both the K and K valleys. We consider intervalley scatterings such that perfect 180 • reflections take place at each junction. Incoming electrons in the K valley are scattered to outgoing electrons in the K valley and vice versa. Insulating states are therefore realized, and the wavefunction is localized within one segment. This clearly indicates that the system allows for a Wannier orbital representation, a hallmark of the topologically trivial states. The two-valley network model with perfect reflection is a topologically trivial insulator. The energy spectrum can be obtained analytically by solving the quantization condition 1 = e i 2Ed v e i(η+η ) , where η and η are the reflection phases at the end points. The energy spectrum is independent of momentum and the flatband energies are at v[−(η + η )/2]/d and v[π − (η + η )/2]/d within one energy period. Now, we discuss how to study the two-valley network model with boundaries. The setup is illustrated in Fig. S5. We consider both K and K valleys, which are connected by a C 2z operation in the bulk, where scattering at the junctions is characterized by a matrix as follows: whereŜ is given by Eq. (1) in the main text, The intervalley scattering in the bulk is completely ignored. In addition, the scattering matrix is the same for both K and K microscopic valleys because of C 2z symmetry. It is worthwhile to note again that the 1D domain-wall states associated with the two valleys have opposite chiralities, as also required by C 2z symmetry. Finally, we consider scattering on the boundaries. The perfect reflecting boundary conditions are imposed on the edge, which can be viewed as the interface between the bulk network model and topologically trivial insulators. The black dots in Fig. S5 indicate those perfect reflecting junctions. The scatterings at the rightmost junctions are given by similarly, the scatterings at the leftmost junctions are given by In the main text, we choose η 1 = η 2 = η 3 = η 4 = 0. The topological properties are independent of the choice of the reflection phases.

EDGE STATE
FIG. S6. Edge state with a different boundary condition. The plot corresponds to Pac = 1. Green arrows and brown arrows indicate K and K valley states respectively. The electron motions associated with the two valleys are related by C2z symmetry (rotation center at AA region), and both are counterclockwise in this plot.
The existence of the edge states can be understood by the skipping orbits in a special limit, P ac = 1 (or equivalently P ab = 1). Here, we provide an alternative way to construct the edge state in Fig. S6. The system is also in a strip geometry with a finite width inx and a periodic boundary condition inŷ. The electrons in the bulk move counterclockwise in localized triangular loops. On the other hand, the right and left boundaries host perfect forward moving edge states because each domain-wall link must contain a pair of counter propagating movers.