{'fraction_non_alphanumeric': 0.04809242342826437, 'fraction_numerical': 0.021527404621171413, 'mean_word_length': 3.8886882285339026, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 40, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fractal energy gaps and topological invariants in hBN/Graphene/hBN double moir\'{e} systems

We calculate the electronic structure in quasiperiodic double-moir\'e systems of graphene sandwiched by hexagonal boron nitride, and identify the topological invariants of energy gaps. We find that the electronic spectrum contains a number of minigaps, and they exhibit a recursive fractal structure similar to the Hofstadter butterfly when plotted against the twist angle. Each of the energy gaps can be characterized by a set of integers, which are associated with an area in the momentum space. The corresponding area is geometrically interpreted as a quasi Brillouin zone, which is a polygon enclosed by multiple Bragg planes of the composite periods and can be uniquely specified by the plain wave projection in the weak potential limit.

Recently, attention is also paid to systems where multiple moiré superperiods compete. The hBN/graphene/hBN stack 23-30 is a typical example, where the moiré pattern caused by graphene and upper hBN layer and that by graphene and lower hBN layer form an incommensurate doubly-periodic potential to graphene as shown in Fig. 1(a). A similar situation is also found in twisted bilayer graphene on hBN [31][32][33] and in twisted trilayer graphene. [34][35][36][37] The hBN/graphene/hBN system is realized when a monolayer graphene is encapsulated by top and bottom hBN substrates. There the dual moiré effect is relevant only when the lattice orientations of upper and lower hBN layers are nearly aligned to graphene, since otherwise the moiré wavelength is too short and hardly affects the low-energy electronic states of graphene. Nearly-aligned hBN/graphene/hBN superlattices were experimentally fabricated using various techniques, [23][24][25][26]29,30 and it was shown that the coexistence of the different super periods gives rise to multiple minigaps in the spectrum, which can never be seen in a single moiré potential. 24,25 Theoretically, double moiré systems are generally hard to treat because the two superlattice periods are incommensurate in general and then the Hamilto- nian is essentially quasiperiodic. The band structure of the hBN/graphene/hBN system was calculated using large-scale numerical simulations 27,28 , where several major gaps and pseudo gaps were found as traces in the energy spectrum against the twist angle.
Here we ask: How can we topologically characterize energy gaps in quasiperiodic systems? In a usual periodic system, the electronic spectrum is separated arXiv:2105.04948v1 [cond-mat.mes-hall] 11 May 2021 into the Bloch subbands accommodating equal electron density, and the number of the subbands below a given gap is a topological invariant. In a doublyperiodic system, however, the absence of the rigorous unit cell prevents the definition of the Brillouin zone, so the topological characterization is not obvious. In one-dimension (1D), an energy gap in a double peroid with wavenumbers G α and G β is characterized by a pair of integers p and q, where the electron density below the gap is given by n e = (pG α + qG β )/(2π). This is regarded the Bragg gap of the (|p| + |q|)thorder harmonics. The integers p and q are directly related to the the topological properties such as the adiabatic pumping [38][39][40][41] and also the quantum Hall effect. 42 In the hBN/graphene/hBN system, similarly, some of the gaps can be associated with the Bragg gap of a composite reciprocal lattice vector, pG α 1 +qG α 2 +rG β 1 +sG β 2 where indeces α, β label the two different moiré patterns. 24,25,27,28,43 This scheme successfully explains a few gaps in the low-energy region, while does not generally work for all the gaps in the spectrum.
In this paper, we calculate the electronic structure of the hBN/graphene/hBN system in changing the twist angle, and identify the topological numbers of all the energy gaps by using a different scheme. First, we compute the band structures for a series of the commensurate approximants to simulate a continuous change of the twist angle. We find that the electronic spectrum actually contains a number of minigaps, exhibiting a recursive fractal structure when plotted against the twist angle. The topological characterization for the energy gaps is employed as follows. Now the system has the four distinct reciprocal lattice vectors G α 1 , G α 2 , G β 1 , G β 2 , and we can define a momentumspace area element (G λ i × G µ j ) z by combining two distinct vectors out of them. As a result, we have four linearly-independent areas A 1 , · · · , A 4 as shown in Fig. 2(a), which can be viewed as projected areas of the four-dimensional hypercube. We find that each energy gap is characterized by a set of integers m 1 , · · · , m 4 such that the electron density below the gap is given by n e = i m i A i /(2π) 2 . Moreover, we show that the area i m i A i is geometrically interpreted as a quasi Brillouin zone, which is a certain polygon composed of multiple Bragg-plane segments as shown in Fig. 2(b). The quasi Brillouin zone for a given gap can be identified by the plain wave projection in the weak potential limit. The bandgap characterization proposed in this work would be useful in other quasi-periodic 2D systems, such as twisted trilayer graphene, twisted bilayer graphene on hBN above mentioned, and also 30 • twisted bi- (a) Independent unit area elements A1, A2, A3, A4 obtained by cross product of the reciprocal lattice vectors G α 1 , G α 2 , G β 1 , G β 2 in hBN/graphene/hBN double-moiré system [Eq. (19)]. The A5 and A6 (dashed areas) can be expressed by others as A5 = −A3 − A4 and A6 = A3. (b) Example of quasi Brillouin zone (thick lines), which is composed of the Bragg planes for composite reciprocal lattice vectors (thin lines). See, Fig. 9 for more details.
layer graphene. [44][45][46][47] The paper is organized as follows. In Sec. II, we define the commensurate approximants and introduce the effective continuum Hamiltonian for the hBN/graphene/hBN system. We calculate the energy spectrum in Sec. III A, and specify the topological numbers of the band gaps in Sec. III B. In Sec. III C, we identify the quasi Brillouin zone associated with the topological numbers by using the plain wave projection. A brief conclusion is given in Sec. IV.

A. Atomic structure
We consider a hBN/graphene/hBN trilayer system as illustrated in Fig. 1, where the top (λ = α) and bottom (λ = β) hBN layers are rotated by θ α and θ β , respectively, relative to the middle graphene layer. Graphene and hBN share the same honeycomb structure with different lattice constants, a ≈ 0.246 nm and a hBN ≈ 0.2504 nm, respectively. 48 We define A and B as sublattices for graphene, and A λ and B λ as nitrogen and boron sites of λ-th hBN layer, respectively. The geometry θ λ = 0 is defined by the AB bond and the A λ B λ bond are parallel to each other.
The lattice vectors of graphene are given by a 1 = a(1, 0) and a 2 = a(1/2, √ 3/2), and those of hBN layers of λ = α, β by where R(θ λ ) is the two-dimensional rotation matrix by θ λ , and M = (1 + ε)1 represents the isotropic expansion by the factor 1 + ε = a hBN /a ≈ 1.018. In the following, we assume the twist angles θ α and θ β are small enough (a few degree or less) such that the moiré super period is much greater than the atomic lattice constant a. The primitive lattice vectors of the moiré pattern of the layer l are given by 18,49 The corresponding reciprocal lattice vectors are where a * i is the reciprocal lattice vectors for graphene which satisfies a i · a * j = 2πδ ij . The moiré superlattice period is given by The moiré rotation angle, or the relative angle of L λ i to a i is given by Figure 3 plots (a) the moiré superlattice period L and (b) the moiré rotation angle φ as a function of the twist angle θ. The super period L is 13.8 nm at θ = 0 • , and it decreases in increasing θ. The rotation angle φ is zero at θ = 0 and rapidly increases in the negative direction in increasing θ.

B. Commensurate moiré approximation
Generally, the two moiré superperiods of α and β are incommensurate and hence there is no unit cell in the trilayer systems as a whole. In any (θ α , θ β ), however, we always have a certain pair of lattice points of the two moiré patterns which happen to be very close to each other. The situation is expressed as where n λ j are integers and ∆L is the difference. When ∆L is much smaller than the moiré periods, the electronic structure of such the system can be  approximated by an exactly-commensurate system with ∆L neglected. Specifically, it is obtained by slightly rotating and expanding / shrinking the moire patterns so that ∆L vanishes. Figure 4(a) shows an actual example of commensurate approximant for (θ α , θ β ) = (0, 1.1908 • ), where (n α 1 , n α 2 ) = (1, 1) and (n β 1 , n β 2 ) = (−1, 3). When ∆L is neglected, Eq. (6) gives a primitive lattice vector of the commensurate super moiré structure, L SM 1 . The other primitive vector L SM 2 is obtained by rotating L SM 1 by 60 • . As a result, we have Correspondingly, the reciprocal superlattice vectors  In the following, we consider two series of hBN/graphene/hBN trilayer systems, In each case, we find a set of (θ α , θ β ) satisfying that ∆L is less than 1% of |n α 1 L α 1 + n α 2 L α 2 | and n α 1 , n α 2 ≤ n max , where n max = 12 and 17 for series I and II, respectively. The full list of (θ α , θ β ) in series I (II) is presented in Table I (II and III) in Appendix A. In series II, the list is dominated by exactly commensurate systems (i.e., ∆L = 0) which appear when the moiré periods of α and β are equal. For later reference, we label those commensuratel cases by [(n α 1 , n α 2 ), (n β 1 , n β 2 )] as p mn ≡ [(m, n), (n, m)], q mn ≡ [(m, n), (m + n, −n)], r mn ≡ [(m, n), (m, −m − n)]. (10)

C. Effective Hamiltonian
Since the hBN has a semiconducting gap, the lowenergy spectrum of the hBN/graphene/hBN system is dominated by the Dirac cones of graphene. We can derive the continuum Hamiltonian of the trilayer system in a similar manner to that for graphene-hBN bilayer. [14][15][16][17][18][19][20][21][22] It is written in 6 × 6 matrix form as which works on the basis of {A, The H G (2 × 2 matrix) is the Hamiltonian for graphene, which is given by where ξ = ±1 is the valley index graphene which correspond to the wave point K ξ = −ξ(2a * 1 + a * 2 )/3, k is the relative wave number measured from K ξ point, and σ ξ = (ξσ x , σ y ) with Pauli matrices σ x and σ y . The H hBN in the second and third diagonal blocks is the Hamiltonian for monolayer hBN, Here we adopt an approximation only considering the onsite potential as 14,18 The off diagonal matrix U λ is the interlayer Hamiltonians of the twist angle θ λ , which is given by 18 where t 0 ≈ 150 meV is the interlayer coupling energy, and r λ 0 is the origin of the moiré pattern of layer λ, which can be changed by sliding the hBN layer relative to graphene. 41 The low-energy effective Hamiltonian for graphene can be obtained by eliminating the hBN bases by the second order perturbation. It is explicitly written as, where with and V 0 ≈ 29 meV, V 1 ≈ 21 meV, and ψ ≈ −0.29(rad). 18 Using the effective Hamiltonian of Eq. (11), we calculate the band structure of the approximate commensurate systems introduced in the previous section. The set of wavenumbers hybridized by the commensurate double moiré pattern is given by q m1,m2 = k + m 1 G SM 1 + m 2 G SM 2 , where m 1 and m 2 are integers and k is a residual wavenumber defined inside the first super-moiré Brillouin zone spanned by G SM 1 and G SM 2 . We construct the Hamiltonian matrix in the bases for graphene, {|q m1,m2 , A , |q m1,m2 , B }, with k-space cut-off |q m1,m2 | < q c . Here we take q c = 2|G β 1 |, which is about 0.54 eV for θ β = 0 • and 1.2 eV for θ β = 2 • . Finally, the band diagram is obtained by plotting the eigenvalues of the Hamitonian matrix as a function of k.

A. Electronic spectrum
As a typical example, we show the band structure of the commensurate approximant for (θ α , θ β ) = (0, 1.1908 • ), which was considered in Fig. 4. Here we set the origins of the moire potentials, r α 0 , r β 0 , to zero. Figure 5(d) shows the energy band plotted along the symmetric line of the super moiré Brillouin zone. For comparison, we also present the band structures (a) with no moiré potential (intrinsic graphene), (b) with only the top moiré potential and (c) with only the bottom moiré potential plotted on the same path. In all the panels, we set the origin of energy (vertical axis) at the Dirac point of graphene. In the single moiré systems in Fig. 5(b) and (c), the biggest gap in the valence band (red / blue regions) is the first order moiré gap corresponding to the electron density of one electron (per valley and per spin) for a moiré unit cell. In the double moiré system, on the other hand, we see a series of the higher order gaps (green) due to the coexistence of the different moiré periods.
To study the twist-angle dependence of the electronic spectrum, we perform the band calculations for all the systems of the series I and II [Eq. (9)]. In commensurate systems, the band structure depends on the relative translation of the moire potentials, ∆r 0 = r α 0 − r β 0 . The dependence on ∆r 0 is generally greater in the systems of smaller L SM , and it quickly vanishes in increasing L SM . Here we average the DOS over 16 grid points of ∆r 0 for the systems with L SM < 50 nm, and otherwise we just take ∆r 0 = 0, since the dependence is minor. (c) shows the lower part of (a), where the firstorder gaps of the single moiré pattern λ = α and β are highlighted by red and blue curves, respectively, and typical higher-order gaps are marked by green curves. Figure 6(b) is the corresponding map of the energy gaps with vertical axis converted to the electron density, where the size of the black dots represent the gap width. In these plots, we see that the spectrum continuously change as a function of the twist angle, even though the adjacent approximants in the series have completely different super moiré periods and thus different numbers of minibands. Figure 7 shows similar plots for the series II plotted against θ α = −θ β . The vertical lines labeled by p mn , q mn , r mn represent the commensurate angles defined in Eq. (10), and the numbers on the top (0, 30, · · · , 120) indicate φ β − φ α , or the relative angle between the two moiré patterns. The r 10 (θ α ≈ 0.5972 • ) and q 01 (θ α ≈ 1.8377 • ) are special cases where the relative angle of the two moiré patterns is 60 • and 120 • , respectively, and hence the two moiré periods completely overlap. There we have a relatively small number of the subbands because of the coincidence of the double-period, but once moving away from these angles, we see that a number of tiny levels branching out just like Landau levels in a magnetic field. As a whole, we observe a recursive pattern ruled by the commensurate lines such as p n,n±1 , q n,n±1 , r n,n±1 . The red dashed curve in Figs. 7 (b) and (c) indicate the positions of the firstorder gaps of the two moiré patterns, which exactly match because of |θ α | = |θ β |. We observe that the first-order gap closes throughout the figure (dashed line) only leaving a small-DOS region around. The reason for the absence of the first-order gap will be explained in the next section.

B. Topological invariants for band gaps
The microgap structure observed in Figs. 6 and 7 resembles the Hofstadter butterfly 50 , which is the energy spectrum of the two-dimensional periodic lattice in magnetic field. The Hofstadter system is essentially equivalent to the one-dimensional Hamiltonian with double period 51,52 , where the fractal minigap structure emerges when the two periods are changed relatively to each other. There each minigap is characterized by a pair of integers p and q, such that the electron density below the gap is given by n e = (pG α + qG β )/(2π) where G α and G β are the wavenumbers for the two periods. The present hBN/graphene/hBN system is a two-dimensional version of this, where the double period is specified by . Actually, as shown in the following, all the gaps observed in Figs. 6 and 7 can be uniquely characterized by four topological integers associated with a specific k-space region.
Let us consider a general situation where the two moire patterns are incommensurate.
We can define four independent unit-areas by combining the four independent reciprocal lattice vectors which are illustrated in Fig. 2(a). Here (· · · ) z represents the z-component perpendicular to the plane, and it can be negative depending on the relative angles between the two vectors. The A 1 and A 2 are the Brillouin-zone areas of the individual moiré patterns of λ = α and β, respectively, while A 3 and A 4 are cross terms which combine the reciprocal vectors of the different moiré patterns. We can also define two more unit areas which are shown as dashed parallelgrams in Fig.  2. In hBN/grahpene/hBN system, however, they are not independent but can be expressed as A 5 = −A 3 − A 4 and A 6 = A 3 , considering that the angle between G λ 1 and G λ 2 is fixed to 120 • . Therefore, a complete set of independent unit areas is given by (A 1 , A 2 , A 3 , A 4 ). The areas A 1 , · · · A 6 can be regarded as the projection of faces of four-dimensional hypercube onto the physical 2D plane, which is analogous to the general argument of the quasicrystal. 53 In a conventional periodic 2D system with primitive reciprocal lattice vectors G 1 and G 2 , the electronic spectrum is separated into Bloch subbands, each of which accomodates the electron density |G 1 × G 2 |/(2π) 2 . In a doubly-periodic 2D system, in contrast, the areas A 1 , · · · , A 4 all serve as units of the spectrum separation. More specifically, we find that the electron density (per spin and valley) from the Dirac point to any gap in the hBN/graphene/hBN system can be uniquely expressed with four integers m 1 , m 2 , m 3 , m 4 as These integers are topological invariants i.e., they never changes as long as the gap survives in a continuous change of the moire pattern. Figure 6(c) shows m 1 , m 2 , m 3 , m 4 found for some major gaps in the case I. Figure 6(b) is the same plot but with the vertical axis being the electron density n e , and the black dots represent spectral gaps with size indicating gap width. Here the integers m 1 , · · · , m 4 are identified from the commensurate approximants as follows. In a commensurate case, A 1 , A 2 , A 3 and A 4 have the greatest common divisor A SM = (G SM 1 × G SM 2 ) z , so they can be written as A i = s i A SM with integers s i (i = 1, 2, 3, 4).
The n e is also quantized in units of A SM /(2π) 2 , and each band gap is characterized by an integer t = n e /[A SM /(2π) 2 ], which is the number of occupied subbands measured from the Dirac point.
is a continuous function of the twist angle, while A SM (and thus t, s i ) can only be defined for commensurate systems and it discontinuously changes in changing the twist angle. This result indicates that the same (m 1 , m 2 , m 3 , m 4 ) are shared by infinitely many commensurate approximants (with A SM ranging from 0 to infinity) which exist in a close vicinity of a specific (θ α , θ β ), and hence it is valid in the limit of A SM → ∞, i.e., incommensurate systems. Figures 7 (b) and (c) are similar plots for the case II. Here the condition |θ α | = |θ β | forces A 1 = A 2 , and then m 1 and m 2 becomes indeterminate. We can resolve the two integers by considering an infinitesimal rotation of either top or bottom hBN layer, and it turns out that m 1 = m 2 for any gaps of the case II. This is explicitly proved as follows. By starting from a case-II system (θ α , θ β ) = (θ, −θ), we can consider two distinct systems X : (θ α , θ β ) = (θ + δθ, −θ) and X : (θ α , θ β ) = (θ, −θ − δθ). The system X and X are actually identical by turning the whole system by 180 • with respect to an in-plane axis, and hence they have the exactly the same energy spectrum. The same energy gap is labeled by a different set of integers as m i and m i for X and X , respectively, which satisfy i m i A i = i m i A i . By considering the layer λ = α, β are interchanged in the 180 • -rotation process, the unit areas of X and X are related by (A 1 , A 2 , A 3   Band structure of a 1D Hamiltonian H = −∂ 2 /∂x 2 + 2V0 cos Gx with V0 = 0, 5, 10. The black solid lines represent the band dispersion ε nk in the extended zone scheme, and the size of blue points represents the spectral weight projected to the plain wave, A(q, ε).
The constraint among the six unit areas A 1 , · · · , A 6 can be broken by uniformly distorting either top or bottom hBN layer such that 120 • symmetry is broken. If we extend the parameter space to such the distorted systems, we should need six topological integers (m 1 , · · · m 6 ) to characterize minigaps, where the electron density is given by 6 i=1 m i A i . This is similar to situation in the series II, where m 1 and m 2 can be resolved by breaking the condition A 1 = A 2 .

C. Quasi Brillouin zones
Actually, the area m 1 A 1 + m 2 A 2 + m 3 A 3 + m 4 A 4 can be associated with a specific region in the momentum space, which is referred to as quasi Brillouin zone. In a conventional periodic 2D system defined by G 1 and G 2 , the Brillouin zones (n = 1, 2, 3 · · · ) are defined by a series of certain regions bounded by the Bragg planes, i.e., the perpendicular bisectors of the reciprocal vectors n 1 G 1 + n 2 G 2 . 54 There all the Brillouin zones have equal area of |G 1 × G 2 |, and therefore the carrier density below any gap is quantized to an integer multiple of the area. In a doubly-periodic 2D system, similarly, we can define a quasi Brillouin zone as an area bound by the Bragg planes for composite reciprocal vectors pG α 1 +qG α 2 + rG β 1 + sG β 2 . In conventional 3D quasicrystals such as Al-Mn alloys, the idea of the quasi Brillouin zones was used to explain the pseudogaps and the stability of the system. 55 In an incommensurate case, generally, the momentum space is filled by infinitely many Bragg planes, and there is no systematic way to define quasi Brillouin zones as in the periodic case. But here, we claim that each single gap in the spectrum can be associated with a specific figure, and the area is equal to m 1 A 1 + m 2 A 2 + m 3 A 3 + m 4 A 4 . Such figures include a simple hexagon defined by a single reciprocal vector as considered in the previous works 25,27,28 , but more generally, it can be a nonconvex polygon composed of multiple segments of different Bragg planes as shown in Fig. 9(a).
The shape of the quasi Brillouin zone for a given gap can be specified by the plain wave projection with the zero potential limit as follows. Let us explain the scheme using a simple one-dimensional Hamiltonian with a single periodic potential, H = −∂ 2 /∂x 2 + 2V 0 cos Gx, where G = 2π. The eigenenergy and the eigenfunctions are labeled as ε nk and |ψ nk , respectively, where n is the band index and k is the Bloch wavenumber in the first Brillouin zone (−π ≤ k ≤ π). Figure 8 shows the band structures calculated for different potential amplitudes, V 0 = 0, 5, 10. The black solid lines represent the band dispersion ε nk plotted in the extended zone scheme, and the size of overlapped blue points represents the spectral weight projected to the plain wave, or where |q = e iqx is the plain wave with −∞ < q < ∞, and the summation in k is taken over the first Brillouin zone. The pink regions indicate the first and the second energy gaps. In decreasing the potential amplitude V 0 , the gaps are narrowing, and the spectral weight approaches a simple parabola ε = q 2 .
In the limit of V 0 → 0, we can specify the points on the parabola, at which the energy gap opens in an infinitesimal V 0 (marked by red circles). These points actually determines the Brillouin zone boundary. The same strategy works for the double-period system as well. In our hBN/graphene/hBN system, we define the spectral weight as where ε α and |ψ α are the eigenenergy and the eigenstates of the system, and |q, X is the plain wave basis of the sublattice X = A, B of the monolayer graphene. For example, we take the commensurate approximant for (θ α , θ β ) = (0 • , 1.1908 • ) considered   Figure 9(c) shows the band structures from η = 0 to 1, calculated on a path from X (graphene's Dirac point) to a certain point Y shown in Fig. 9(a). The black solid lines represent the band dispersion plotted in the extended zone scheme, and the blue dots represent the spectral weight A(q, ε). At η = 0, we just have the graphene's Dirac cone. By tracing the gaps in the spectral weight in decreasing η from 1 to 0, we can specify the gap opening points just as in the one-dimensional case.
In Fig. 9(c), we consider four gaps with different indeces of (m 1 , ..., m 4 ). The (−1, 0, 0, 0) is the first-order gap of the moiré potential λ = α, and others are double-moire gaps caused by the coexistence of the two moiré patterns. In the limit η → 0, we find the gap-opening wave numbers P 1 , · · · , P 4 for these gaps. By following the same procedure for paths in different directions, we finally obtain the quasi Brillouin zone on the (k x , k y ) plane as the traces of P 1 , · · · , P 4 , which are illustrated as thick colored lines in Fig. 9(a). The figures are composed of segments of the Bragg planes, which are shown as thin lines. The first-order gap (−1, 0, 0, 0) gives a regular hexagon, which is the first Brillouin zone of the moiré potential of λ = α. The double-moire gap (−1, −1, −1, 1) also gives a hexagon but with a smaller size, which corresponds to the first Brillioin zone of a small reciprocal lattice vectors G α 2 +G β 1 . In contrast, the gaps (2, 1, 2, −2) and (−4, −2, −2, −2) are associated with flower-like complex figures composed of multiple Bragg line segments. In any cases, the area of the figure is shown to be exactly equal to m 1 A 1 + m 2 A 2 + m 3 A 3 + m 4 A 4 . Just as the conventional Brillouin zone in a periodic system, the quasi Brillouin zone is also a closed object, in that any sides of the boundary are precisely sticked to the other side and one can never go out of the region by crossing the boundary.
The quasi Brillouin zone continuously changes in changing the twist angle, regardless of the unit cell size of the commensurate approximants. Figure 10 shows the same plot calculated for a slightly different angle (θ α , θ β ) = (0, 1.2967 • ) [(iii) in Fig. 6]. The super moiré unit area of the system is about 10 times greater than that of Fig. 9(c), and accordingly we see much more band lines due to the band folding into the smaller Brillouin zone. If we see the spectral weight (blue dots), however, we find that it exhibits a similar structure to Fig. 9(c) (except that the gap (−1, −1, −1, 1) is not fully open), and the gaps close at the Bragg planes with the same indeces in the limit of η → 0. As a result, we end up with nearly the same shape of the quasi Brilluoin zone as shown in Fig. 11(iii). In Fig. 10, we see a number of extra band lines are just overlapping but hardly contribute to the spectral weight, and therefore they are neglected in the identification of the zone boundary. Because of this, the quasi Brillouin zone obtained here is generally different from one obtained by sort-ing all the eigenvalues in energy and tracking the same level index in the limit of the zero potential 56 , which is fully affected by all the overlapping band lines.
In Fig. 11, we show the continuous evolution of the quasi Brillouin zones as a function of the twist angle from (i) to (v) [corresponding to the labels in Fig.  6], where the figure continuously changes regardless of the discontinuous change of the rigorous period of the approximants. The areas of the figures are always equal to

IV. CONCLUSION
We theoretically studied the electronic structure of the hBN/graphene/hBN double-moiré system as a function of the top and bottom twist angles, and demonstrated that the spectrum consists of a number of fractal minigaps characterized by the non-trivial topological numbers. Specifically, each energy gap is characterized by a set of integers (m 1 , · · · , m 4 ), where the electron density below the gap is given by n e = i m i A i /(2π) 2 with characteristic momentum-space areas A 1 , · · · , A 4 . The area i m i A i corresponds to a quasi Brillouin zone bounded by multiple Bragg planes, which can be uniquely identified by the spectral distribution in the zero potential limit. In changing the twist angles, the quasi Brillouin zone also changes continuously regardless of the commensurability of the double moiré pattern.
We neglect the lattice relaxation effect throughout this work for simplicity, while the general theoretical scheme to characterize the band gap is valid as long as the system has a well-defined double period. The topological band-gap characterization proposed in this work should also be useful in other quasi-periodic 2D systems, such as twisted trilayer graphene [34][35][36][37] , twisted bilayer graphene on hBN [31][32][33] , and also 30 • twisted bilayer graphene. [44][45][46][47]

V. ACKNOWLEDGMENTS
This work was supported in part by JSPS KAK-ENHI Grant Number JP20H01840 and JP20H00127 and by JST CREST Grant Number JPMJCR20T3, Japan.  Quasi Brilluoin zones similar to Fig. 9(a) calculated for five different angles. The indeces (i) to (v) correspond to the labels in Fig. 6. TABLE I. List of approximate commensurate systems of of series (I)(θ α , θ β ) = (0, θ). The L SM [nm] is the super moiré period, and ∆L[nm] is the correction from the original incommensurate structure (see the text).   Table  III).   Table II).