Replica Higher-Order Topology of Hofstadter Butterflies in Twisted Bilayer Graphene

The Hofstadter energy spectrum of twisted bilayer graphene is found to have recursive higher-order topological properties. We demonstrate that higher-order topological insulator (HOTI) phases, characterized by localized corner states, occur as replicas of the original HOTIs to fulfill the self-similarity of the Hofstadter spectrum. We show the existence of the exact flux translational symmetry of twisted bilayer graphene at all commensurate angles. Based on this result, we carefully identify that the original HOTI phase at zero flux is re-entrant at a half-flux periodicity, where the effective twofold rotation is preserved. In addition, numerous replicas of the original HOTIs are found for fluxes without protecting symmetries. Similar to the original HOTIs, replica HOTIs feature both localized corner states and edge-localized real-space topological markers. The replica HOTIs originate from the different interaction scales, namely, intralayer and interlayer couplings, in twisted bilayer graphene. The topological aspect of Hofstadter butterflies revealed in our results highlights symmetry-protected topology in quantum fractals.

Notably, the link with the magnetic translational symmetry and symmetry-protected topological phases of matter has been revealed recently [28][29][30][31][32][33][34][35][36][37].In a general lattice model with multiple sites per unit cell, the Hofstadter energy spectrum becomes approximately replicative under the addition of the flux periodicity, E(ϕ) ≈ E(ϕ+Φ), which constitutes the additional flux translational symmetry via the unitary transformation of the Hamiltonian, H(ϕ), as, where ) is a creation (annihilation) operator of an electron at R in the real space, A is the vector potential, and S is the unit cell area.Remarkably, the effective time-reversal symmetry UT is restored at a half-flux periodicity ϕ = 1  2 Φ [30], allowing for the existence of diverse topological states of matter protected by symmetries [28][29][30][31][32].
In this work, we study the higher-order topological insulator (HOTI) phases of Hofstadter butterflies in TBG.Archetypal HOTIs have been studied with respect to symmetry protection [28,30].By contrast, the replica HOTIs that we find here recur in the form of quasiperiodic replicas without explicit symmetry protection.Instead, they rely on the self-similar nature of Hofstadter butterflies.We prove that the full lattice model of TBG possesses the exact flux periodicity at all commensurate angles, which rigorously characterizes the band topological protection in the presence of the magnetic field.Two original HOTIs exist at time-reversal invariant fluxes (TRIFs) ϕ = 0 and ϕ = 1  2 Φ, where ϕ = −ϕ (mod Φ).In addition, replicas HOTIs recur at the specific fluxes ϕ = p 2Nrep Φ (p ∈ Z; p ̸ = N rep Z) [Fig.1b; See Eq. ( 3) for the definition of N rep ].To quantitatively diagnose HOTIs, we extend the concept of real-space topological markers [38][39][40][41][42][43]

Lattice model and symmetries
We use the Moon-Koshino tight-binding model for TBG [44] where the hopping integral t ij (R i − R j ) is modelled as an exponentially decaying function of [44] (see Methods).Magnetic flux ϕ is introduced using the Peierls substitution [45]: R j A 0 • dr , where A 0 is the vector potential in the Landau gauge (see Methods).We consider the atomic structure of TBG in the hexagonal space group # 177, generated by twisting the AA-stacked bilayer graphene about the hexagonal center with the twist angle θ m,n = arccos 1 2 1a).This construction of TBG preserves C 6z , C 2x , and C 2y rotational symmetries.The twist lowers the discrete translational symmetry, leading to the translational symmetry of the moiré lattice with the enlarged unit cell area by N L = m 2 +n 2 +mn times.In the presence of a uniform perpendicular magnetic field, a flux translational symmetry emerges, which locally restores crystalline symmetries for specific fluxes.

Hofstadter butterflies
For the nearest-neighbor tight-binding model of graphene, the flux periodicity is given as the magnetic field strength B = Φ 0 /S G , where Φ 0 = h e is the flux quantum and S G is the graphene unit cell area [11,44,46].However, when next neighbor hoppings are introduced, the minimal loop along the allowed hoppings, namely, the minimal Peierls path, has decreased inner area S G min = 1 6 S G (Fig. 1a), leading to an increased flux periodicity.A stronger magnetic field of B = Φ 0 /S G min = 6Φ 0 /S G is required to implement the full flux quantum into the decreased inner area S G min of the minimal Peierls path.Consequently, the entire cycle is completed by repeating six times modulated quasiperiodic replicas of the nearest-neighbor graphene spectrum [30].
For TBG, we show the existence of the exact flux periodicity at the twist angle θ m,n , dictated by, where gcd indicates the greatest common divisor and z (see Supplementary Note 2).For θ = 21.8 • (m = 1,n = 2), N rep = 7 corresponds to the area of the minimal Peierls path, S TBG min ≡ S = 1 7 S G min = 1 42 S G (see Fig. 1a).As a result, a self-similar pattern is rendered by 42 replicas of the original graphene spectrum, only having the nearest-neighbor hopping term.
Figure 2 shows the calculated Hofstadter butterflies for both graphene and TBG by using the kernel polynomial method (see Methods).Quasi-periodicity is exhibited, as our tight-binding model includes electron hopping beyond the nearest neighbors.For example, in the graphene spectrum (Fig. 2a), the quasi-periodicity of 1  6 Φ is displayed by having similar patterns recurring at every integer multiple of 1  6 Φ.Similarly, for TBG spectrum (Fig. 2b), a quasi-periodicity of 1 42 Φ occurs as expected.Moreover, the energy spectrum that resembles the graphene spectrum in Fig. 2a recurs at every integer multiple of 1  7 Φ.This modulation of the graphene spectrum by 1  7 Φ is weaker than 14 Φ, and j ϕ = 3 7 Φ, where ε = 1 1680 Φ.
that of 1 42 Φ because the interlayer hopping is relatively weaker in TBG compared to next-nearestneighbor intralayer hopping.Therefore, the quasi-periodicity of 1  7 Φ is more prominent than that of 1  42 Φ in TBG spectrum.The computed spectrum exhibits symmetries of Hofstadter butterflies (Fig. 2b).Translational flux symmetry is displayed in the recurring patterns at ϕ and ϕ = ϕ + Φ.Moreover, the C 2x symmetry that is broken under the flux gives rise to the mirror-symmetric spectrum about TRIFs (both ϕ = 0 and ϕ = 1 2 Φ).The Hamiltonian is transformed under the C 2x operator as Combined with the unitary matrix U, we obtain Therefore, the energy eigenvalues for ϕ and −ϕ about TRIFs are equivalent.

Exact HOTIs
The proposed tight-binding model reproduces the HOTI phase of TBG well at zero flux, showing good agreement with previous studies [47,48].Consequently, the system harbors localized states at the corner of a diamond-shaped flake under an open boundary condition (OBC) (Fig. 3b).In energy space, two corner states reside inside the spectral gap of the bulk (Fig. 3a).In general, these two (in-gap) corner states can have different energies owing to the finite-size effect, in which they spatially overlap and cause hybridization [48].
We suggest a HOTI marker given by χ ± (r) = − ⟨r| C ± 2x P ± XQ ± |r⟩, where C ± 2x = P ± C 2x P ± is a projected symmetry operator, and X is a position operator (see Methods).Here, is the projection operator to the occupied (unoccupied) c 2x = ±1 subspaces.In OBC, χ ± (r) successfully diagnoses the rotation-winding number in real space: χ ± (r) dictates the nontrivial rotation-winding number by being localized along the edge of the flake (Fig. 3d), whereas in the trivial case, it is delocalized over the entire geometry (see Supplementary Figure 3).Interestingly, the corner state appears at the boundary between the opposite signs of each HOTI marker χ ± (r).
The sum of the opposite HOTI markers χ + (r) + χ − (r) is zero, which indicate a trivial winding number.The HOTI marker can be applied to symmetry-breaking perturbations, as demonstrated in TBG under the uniform magnetic field.
To study the effect of the magnetic field on the HOTI states, we track the corner states by investigating their spectral flow at fixed filling [30] (see red and green lines in Fig. 2b).At zero flux, the highest occupied (HO) and lowest unoccupied (LU) states are identified as corner-localized states (Fig. 3b).They adiabatically evolve as a flux function and undergo a series of discontinuity transitions at specific fluxes.This discontinuity is indicative of a topological change due to bulk gap change [30,31].Indeed, we reveal that HO and LU states at the discontinuity transitions are quantum Hall chiral edge states (see Supplementary Note 4).
Remarkably, we find a reentrance of the HOTI phase at ϕ = 1 2 Φ, characterized by edge-localized marker χ + (r) (Fig. 4a).χ + (r) decays exponentially along the bulk as exp[−α(x − x 0 )] with α = 0.52, which is identical to that of the exact HOTI state at zero flux (Fig. 4b) (see also Supplementary Note 3 for the detailed quantitative analysis).The re-entrant HOTI phase relies on composite symmetry UC 2x exactly preserved at ϕ Note that the corner boundary modes of the re-entrant HOTI phase are localized at the corner, but the node appears slightly more concentrated off the corner (Fig. 4a) (see also Supplementary Figure 4 for the reason of the nodal structure of the corner states).

Replica HOTIs
In addition to the exact HOTIs at TRIFs (ϕ = 0 and ϕ = 1 2 Φ), replicas of the original HOTIs are found at the 1  7 Φ quasi-periodic counterparts of TRIFs.We employ HO and LU states as indicators of a replica of the original HOTI.We find that they are positioned within the spectral gap at the specific fluxes of quasi-periodicity ϕ = p 14 Φ (p ∈ Z; p ̸ = 7Z) (Figs.2e-j).A close inspection reveals that HO and LU states show oscillatory behavior of HO and LU energies as a function of flux, which originates from the Aharonov-Bohm tunneling in the presence of an external magnetic flux.
Notably, the oscillation is a finite-size effect rather than a characteristic behavior of corner states, as is evident in the oscillations of other states near HO and LU states.
To demonstrate the characteristics of the replica HOTIs, we plot the HO states in the left panels in Figs.4d-i.The real-space distribution arguably shows the corner-localized states, supporting the HOTI phases.Nonetheless, these states exhibit stark contrast to the corner states of the exact HOTI at zero flux in that they are complex-valued functions, while the exact HOTI hosts real-valued corner states (Fig. 3b).The complex-valued wave functions manifest the broken reality condition [(C 2z T ) 2 = 1] at finite fluxes, implying that the Stiefel-Whitney characterization is inapplicable.Furthermore, these quasiperiodic fluxes also break the C 2x and UC 2x symmetries, which were utilized to characterize the exact HOTIs at TRIFs.
Remarkably, the HOTI marker can be defined without the protecting symmetries, enabling the evaluation of rotation-winding numbers.We find that the HOTI marker can quantitatively characterize the corner states in the presence of flux, that is, under rotational-symmetry C 2x breaking.At a small flux ϕ 1 = 1 21000 Φ, the eigenstate shows the remaining localized corner state, and the corner state is characterized by the marker χ + (r) which is sufficiently localized along the entire edge despite the small permeated values towards the bulk (Supplementary Figure 3).Quantitatively, χ + (r) exhibits an exponential decay as exp[−α(x − x 0 )] with α = 0.10, which is smaller than α = 0.52 of the exact HOTIs due to the symmetry breaking (Fig. 4b).The exponential localization of χ + (r) from the edge for the corner states is in stark contrast to a linear delocalization ∝ −β(x − x 0 ) of χ + (r) along the whole geometry for the trivial state that occurs at, for example, ϕ = 13ϕ 1 (Fig. 4b).Such localization characteristics of the markers serve as a hallmark to identify nontrivial bulk topology, which fundamentally originates from the action of the projected symmetry operator, as in the generic topological crystalline insulating phases protected by spatial symmetries [42,43] (see also Methods for the detailed explanation for the real-space behavior of the HOTI marker).
Our HOTI marker captures the replica HOTI phases as well, at ϕ = p 14 Φ (p ∈ Z; p ̸ = 7Z) (Figs.4d-i).χ + (r) at ϕ = p 14 Φ show robust edge localization, consistent with the corner-localized eigenstates.The line profiles of χ + (r) (Fig. 4c) exhibit exponential decay (see also Supplementary Figure 6).The replica HOTIs can be viewed as the copies of exact HOTIs disordered by the fractional flux quantum acquired when electrons travel through the minimal Peierls path because the composite symmetry U 0 C 2x becomes exact when the interlayer coupling is turned off.We also verify that replica HOTIs generally appear at other large angles.Figure 5 shows the HO states and HOTI markers at the other twist angles θ = 13.2 • (m = 2, n = 3) and 9.4 • (m = 3, n = 4).We find that both the corner localized states in real space and the localization behavior of the calculated HOTI markers support the existence of the replica HOTI states at the flux where the flux periodicity is given by Φ m,n = N rep Φ G with the flux periodicity of graphene Φ G .Here, N rep = 19 and 37 for θ = 13.2 • and 9.4 • , respectively.We note that the localization strength of the HOTI markers (see the line profiles in Fig. 5) is weakened as we decrease the twist angle because the bulk gap is significantly reduced (see Supplementary Figure 9).We find that the out-of-plane rotational symmetry C 2x is essential to realize the re-entrant exact and replica HOTI phases in TBG under a magnetic field.In contrast to our model, there is no re-entrant corner state at half-flux periodicity in the magic-angle TBG model [30] with only C 2z T symmetry, where the flux pumps corner states into the bulk.The disappearance of the corner states at half periodicity confirms the inapplicability of the Stiefel-Whitney characterization for the HOTI states in the presence of a magnetic field.This indicates that additional crystalline symmetry, such as C 2x , is required to protect the corner states in TBG under a strong magnetic field.
In summary, we have demonstrated that HOTIs can occur without explicit protecting symmetries because of the self-similarity of Hofstadter butterflies as replicas of original HOTIs.We expect the distinct symmetry dependence of replica HOTIs can lead to distinct physical properties from the exact HOTIs (see Supplementary Note 6 for the detailed discussion).The HOTI marker is an invaluable tool for studying HOTI states in various situations beyond conventional methods using periodic boundary conditions.It offers the distinct advantage of being able to readily identify the HOTI phase, even at a small magnetic field in the open boundary condition.This is particularly advantageous compared to momentum-space methods relying on periodic boundary conditions, as they are computationally demanding at low magnetic fields, with their computational cost scaling inversely with the strength of the magnetic field.The exponents of our HOTI marker allows for quantitative analysis, which can be potentially useful for future study such as many-body disordered HOTIs.The observation of the replica HOTI at the fixed filling requires a huge magnetic field B ∼ 10 5 T, but replica topology may occur at different filling near low fields.Therefore, establishing an exact relationship between discrete scale invariance and band topology in this quantum fractal will be exciting future research with direct experimental implications.Additionally, a critical challenge that needs to be tackled in order to realize the observation is ensuring the stability of large TBG flakes under high magnetic fields.It would also be interesting to explore the Coulomb repulsion effect on replica phases at smaller angles, where the role of Coulomb repulsion is cru-cial [59][60][61][62][63][64].With much progress in synthesis of moiré materials [65][66][67][68][69][70] and measurement of Hofstadter energy spectrum [33,71], our results can pave the way for studying replica topology under magnetic field in generic moiré multilayer [67][68][69] and moiré quasiperodic [70] systems that host multiple interaction scales.

Tight-binding model
We employ the Moon-Koshino tight-binding model for twisted bilayer graphene in Ref. [44], which is written as where c † R i (c R j ) is a creation (annihilation) operator of an electron at the lattice site R i , and t ij (R i − R j ) is the hopping integral between the sites R i and R j .The hopping integral is given by Here, the hopping parameters are given as a decaying function of a hopping distance where a 0 ≈ 1.42 Å is the bond length of graphene, d 0 ≈ 3.35 Å is the interlayer distance, and δ = 0.319a 0 is the decay length.Here, we set V 0 ppπ = −2.7 eV and V 0 ppσ = 0.82 eV, which reproduce the band structure of 21.8 • twisted bilayer graphene with a bulk gap ∼ 9 meV in a HOTI state [47,48].Our tight-binding model under a periodic boundary condition (see atomic geometry used in Supplementary Figure 1) has 14 occupied p z orbital bands that consist of the same number of c 2x = +1 and c 2x = −1 bands, where c 2x is an eigenvalue of a twofold rotational symmetry operator C 2x about the x-axis.This implementation of the model successfully reproduces the nontrivial rotation-winding number (Fig. 3c) in line with the previous DFT results [47].For the calculations of Hofstadter butterflies and topological markers, we use the flake geometry with an open boundary condition (see Supplementary Figure 1).
To study the effect of the magnetic field, we incorporate a magnetic flux ϕ into the hoppings as an additional phase via Peierls substitution [45]: where the vector potential A 0 = B(0, x) = ϕ S min (0, x) for B = B ẑ and S min is the interior area of the minimal Peierls path.We prove that our twisted bilayer graphene lattice has exact flux periodicity (see Supplementary Note 2).The Hofstadter energy spectrum of our system is thus periodic under the translation by a magnetic flux quantum Φ = h e because the Hamiltonian can be gauge transformed according to [30] The unitary matrix r 0 A • dr) (r 0 : a fixed lattice site) is defined for the vector potential A(̸ = A 0 ) that leads to the flux quantum S min (∇ × A) • dS = Φ.

Kernel polynomial method
Hofstadter butterflies of twisted bilayer graphene can be efficiently calculated by using the kernel polynomial method [72].The essential idea of the kernel polynomial methods is to expand the density of states ρ(E) (E: energy) in terms of Chebyshev polynomials as, where U n (E) is the second kind n-th Chebyshev polynomials, Here, µ n is the moment for an operator Ô(E), which reads The targeting density of states operator ρ(E) is given by After putting ρ(E) into µ n , we obtain A stochastic approach is employed to obtain the trace by introducing the R-number of random vectors |r⟩, instead of (potentially unknown) exact eigenvectors: where R is set to a sufficiently large value to attain the converged density of states.Then, we take advantage of a recursive relation for the polynomial, to rewrite the trace as where As a result, the density of states is obtained as where g M m is the Jackson kernel, which is introduced to reduce the Gibbs oscillation [72].

HOTI topological marker
Topological marker is a local quantity in real space that characterizes the topological phases [38][39][40][41][42].The local Chern marker was first introduced as a topological marker whose spatial average in bulk in thermodynamic limit corresponds to the Chern number of the system [38].The topological marker was then generalized to the topological crystalline insulating (TCI) phases, in which the topological states are protected by the spatial symmetries [42].The generalized topological marker T G (r) related to the symmetry G is given by where G = P GP is a projected symmetry operator and a function F(P ) encodes the types of topological invariants.For example, F(P ) ∝ P [ X, P ] and P [[ X, P ], [ Ŷ , P ]] for 1D winding [39] and 2D Chern numbers [38], respectively, where X and Ŷ are position operators.
We extend the topological marker to a HOTI version in our twisted bilayer graphene system.
The extension is straightforward because the HOTI phase in twisted bilayer graphene is protected by the C 2x rotation symmetry resolved winding number, the rotation-winding number.Let us first see the topological marker χ (r) for the C 2x symmetry, which is given by where we used the relation Q = 1 − P in the last equality.By projecting the projection operators to the C 2x rotation ± subspaces as P = P + + P − and Q where we used the condition The C 2x rotation-resolved topological marker χ ± (r) serves as the real space local expression of the rotation-resolved Zak phase ν ± .

Real-space behavior of HOTI marker
To understand the real-space behavior of the HOTI marker, we first consider the localization property of the topological markers for TCI phases.The topological markers for TCI phases feature the exponential localization from the subspace restricted by the spatial symmetries [42,43].It is different from the case of the typical Chern insulators without symmetries where the localization sites of the topological marker are all the sites within the bulk [38,40,41].In a TCI phase, protected Here the length scale ζ is rough in the order of the inverse gap/localization strength.We note that the localization property is fundamentally arising from the action of the projected symmetry operator G: the projection matrix P is exponentially localized for the insulators [73][74][75] and the symmetry G restricts the localization site of the markers [42,43].The localization property is more general than the exponentially localized Wannier functions because the topological marker is localized even in the presence of a nonzero Chern number, which prohibits the construction of localized Wannier functions.
In the case of C 2x winding number, the presence of the projected symmetry operator C ± 2x gives rise to the exponential localization of the HOTI marker χ ± (r) from the edge.The absence of the winding number allows for the C 2x specification of the winding number, revealing the presence of the C 2x -protected metallic edge states.As in the case of the known HOTI phases characterized by C 2x winding number [47,48], the corner state is the Su-Schrieffer-Heeger type domain wall state, arising from the gap opening of the edge states.The edge is the bulk of the corner states in the On the contrary, HOTI trivial cases do not show such localization behavior from the edge.Instead, they are linearly delocalized over the entire geometry (∝ the position operator X) which follows from the form of the C 2x -marker formula proportional to X in Eq. 24.
given as,  We explicitly plug in Eq. ( 4).Since each x, y coordinates are the linear polynomial of (p 1 , p 2 , p 3 , q 1 , q 2 , q 3 ) ∈ Z 6 , the expression for the area of the triangles can be written in terms of the quadratic polynomials of (p 1 , p 2 , p 3 , q 1 , q 2 , q 3 ), S = √ 3 4N L |z 0 + i z p,i p i + i z q,i q i + i̸ =j z pp,ij p i p j + i̸ =j z pq,ij p i q j + i̸ =j z qq,ij q i q j | (6) where z are the non-zero integer coefficients of the polynomial.There are total 31 z coefficients and the explicit expression will be given later.To further simplify the above expression, we use the theorem of the linear Diophantine equation [1], which states as, By sequentially applying Theorem 1 to Eq. ( 6), we find that the set of the areas of the possible triangles is given as,

FIG. 1 .
FIG. 1. Replica HOTI states in TBG. a Atomic structure of TBG. S denotes the smallest area enclosed by a Peierls path for TBG at θ=21.8 • .For comparison, we illustrate the corresponding area for graphene S G min and the area of the graphene unit cell S G , where S = 1 7 S G min = 1 42 S G .b Schematic of re-entrant exact HOTI and replica HOTI phases as a flux ϕ function.Red and blue colors indicate the exact HOTI and replica HOTI phases, respectively.The k x and k y are pseudo-momenta, well defined in the presence of the magnetic translational symmetry under arbitrary rational flux.

FIG. 2 .
FIG. 2. Hofstadter butterflies of the atomistic tight-binding model of TBG.a-b Hofstadter spectrum of TBG a without and b with interlayer coupling calculated by using the kernel polynomial method (see Methods).Red dashed lines indicate the system's half and full flux periodicity.Blue dashed lines indicate quasi periodicities at ϕ = p 14 Φ (p = 1, 2, . . ., 6).In b, the red and green solid lines denote the highest occupied and lowest unoccupied states, respectively.c-j Magnified view of b at specific flux values c ϕ = 0,

FIG. 3 .
FIG. 3. Characterization of HOTI states.a Energy spectrum of TBG with θ = 21.8 • in OBC.Red Data points indicate corner states.b Topological corner state in OBC.The colored circle indicates the phases of eigenstate components.c Zak phase calculated using c 2x = +1 bands along the rotation-invariant k y = 0 line in periodic boundary condition.d HOTI marker χ + (r) for c 2x = +1 states in OBC.Here, χ + (r) is normalized to its maximum value.
by a spatial symmetry G, the eigenvalues of G classify the eigenstates of the Hamiltonian {|u i ⟩} and thus the projection matrix for occupied states P = n∈occ.|u n ⟩ ⟨u n | at the symmetry-invariant subspace S. As a result, the bulk topology of a TCI phase is encoded by the projection matrix P at the invariant subspace S, which allows the introduction of the real-space topological invariant, the topological marker T G (r) in Eq. 22.It is proven that the topological marker T G (r) exhibits the exponential localization from the fixed points r S (∈ S; S = {r|Gr = r}) of the spatial symmetry G as [42] |T G (r)| < O(e −|r S −r|/ζ ) when |r S − r| ≫ ζ.
HOTI state.Due to the action of the projected symmetry operator C ± 2x , our HOTI marker of HOTI states exhibits exponential localization from the edge r E as | χ ± (r)| < O(e −|r E −r|/ζ ) when |r E − r| ≫ ζ.

( 4 )
To find the minimal area of the triangular Peierls path, without loss of generality, we only need to consider the triangles, where the two of the vertices consist of the top layer sites [r t,1 = (x 1 , y 1 ), r t,2 = (x 2 , y 2 )] and the other consists of the bottom layer site[r b = (x b , y b ) = R(θ)(x 3 , y 3 )].The area of such triangles can be written as, S = 1 2 |x t,1 (y t,2 − y b ) + x t,2 (y b − y t,1 ) + x b (y t,1 − y t,2 )| (5)