Intrinsic superflat bands in general twisted bilayer systems

Twisted bilayer systems with discrete magic angles, such as twisted bilayer graphene featuring moiré superlattices, provide a versatile platform for exploring novel physical properties. Here, we discover a class of superflat bands in general twisted bilayer systems beyond the low-energy physics of magic-angle twisted counterparts. By considering continuous lattice dislocation, we obtain intrinsic localized states, which are spectrally isolated at lowest and highest energies and spatially centered around the AA stacked region, governed by the macroscopic effective energy potential well. Such localized states exhibit negligible inter-cell coupling and support the formation of superflat bands in a wide and continuous parameter space, which can be mimicked using a twisted bilayer nanophotonic system. Our finding suggests that general twisted bilayer systems can realize continuously tunable superflat bands and the corresponding localized states for various photonic, phononic, and mechanical waves.

In this work, we discover the robust presence of a class of superflat bands in general twisted bilayer systems proved by the tight-binding model (TBM) with negligible next-nearest-neighbor intralayer hoppings. Using the effective macroscopic potential well model (PWM) with spatially modulated couplings, we show that for small twists, localized states definitely appear centered on the AA stacked region (with deepest potential well) at isolated lowest and highest energies, manifesting C 6 and C 3 symmetries, respectively. Such localized states present negligible inter-cell coupling, forming superflat bands for general twisted bilayer systems, which is corroborated by exact TBM calculations. We further implement superflat bands and the corresponding localized states via twisted bilayer nanophotonic platforms. Importantly, these superflat bands arise for a continuous set of small angles and do not require fine tuning to the specific magic angles, being readily implementable for various wave systems and introducing an extremely large density of states (DOS) for lasing 34 , sensing 35 , and light-matter interactions 36 .

Superflat bands and localized states
General twisted bilayer systems display alternating patterns between AA and AB/BA stacked lattices (i.e., the A (B) site from the upper layer is perfectly aligned with the A/B (A) site from the lower layer), as illustrated in Fig. 1a. In momentum space, rotated unit cells in two layers cause a relative rotation (θ) of first Brillouin zones (BZs), generating an effective moiré BZ (see Fig. 1b). Periodic moiré superlattice has the lattice constant a M ¼ a 2sinðθ=2Þ , where a is the lattice constant of primitive unit cells (with a hexagonal p6m symmetry of wallpaper groups). We assume that the hopping rate between every two sites (i ≠ j) decays exponentially as a function of distance |r ij |, i.e., t ij $ A 0 e Àγjr ij j , because the classical electronic and photonic systems always allow the overlap of exponential-type wave functions [37][38][39][40] . Here γ represents the decay rate and A 0 is the normalized coefficient constraining the energy scale. In addition, negligible nextnearest-neighbor hoppings of intralayer sites restrict the range of t ij in the following form Without losing generality, we set the unit hopping t 0 = 1 in the following analysis. To ensure the dominance of nearest-neighbor hoppings accurately, we further choose γa~30 corresponding to t ij (a)~10 −5 ≪ t 0 . An exact hopping strength curve is displayed in Fig. 2a, where the spatial distance jr ij j ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi l 2 þ ρh 2 q , l and h represent the intralayer and interlayer distances, respectively. ρ = 0 (ρ = 1) stands for i and j located at the same (distinct) layers. We model general spinless twisted bilayer systems with the TB Hamiltonian where c ðyÞ i corresponds to the creation (annihilation) operator at the site i, and ϵ is the inherent potential which is considered as zero in general systems. This allows us to perform the exact analysis for moiré superlattices and provide numerical support for the following detailed models.
Furthermore, we calculate the AA and AB/BA stacked band structures under the above hopping relation using the analytical TBM 41,42 . The lowest/highest energies of first BZs (located at Γ point, i.e., the center of blue and red   . The insets show the field distributions of S 1 (highest energy) and S 2 (lowest energy) located at Γ point, implying C 3 and C 6 symmetries, respectively hexagons in Fig. 1b) can be reduced to Algebraic derivation reveals that for highest (or lowest) bands AA stacked lattices always have higher (or lower) energies than AB/BA stacked lattices, forming a natural potential difference, i.e., jE AA Such a relevant energy difference provides a spatial potential well where the deeper potential is located at the AA stacked region with effective masses m*~±2ℏ 2 /t 0 . These effective masses remain almost constant at arbitrary locations within the moiré superlattices, with the details given in Supplementary Note 1. Here we show a specific case with h ¼ a= ffiffi ffi 3 p (see Fig. 2b), where band structures of AA and AB/BA stacked lattices match well with our analysis. Two states (S 1 and S 2 ) with highest/lowest energies at Γ point of AA stacked lattices present C 3 and C 6 symmetries, respectively, preserved by irreducible representations in the orthogonal eigenspace, which are the crucial prerequisite for forming superflat bands as following discussions.
In the vicinity of lowest/highest energies of AA stacked lattices, the previous low-energy theory describing moiré bands is invalid 29 . A concise physical picture can be constructed to depict this system as illustrated in Fig. 3. The distorted lattices along the azimuth θ c = nπ/3, n = 1, 2, . . . , 6, centered around AA stacked lattices, reflect essential characteristics of the potential well. Specifically, for the distorted lattice with a distance from the center of AA stacked region r o , the coordinates of lattice center are (c x , c y ) = r o (cos(θ c ), sin(θ c )). The geometric center of A and B sites is shifted and projected on a specific circle with radius r c = 2r o sin(θ/4). The distance in x-y plane from one center to another center for two layers is d c = 2r o sin(θ/2). Here, r c and d c are independent of θ c . In the vicinity of AA stacked region, dislocated lattices for any r o and θ c allow for modeling on a scale of unit cells. A typical case for n = 0 is displayed in Fig. 3 (right panel). The Hamiltonian around Γ point characterizing lattice distortions of the system, Φ ¼ fϕ 1 A ; ϕ 1 B ; ϕ 2 A ; ϕ 2 B g, takes the form 43 ffiffi ffi 3 p a t 1 þt 2 2 k y Þ and the wavevector k = {k x , k y }. σ x,y are the Pauli matrices acting in sublattice space of single layers. t 1 and t 2 correspond to inter-cell hoppings between A and B sites for single layers along two distinct basis vectors, respectively, which are equal for zero θ or unequal (and exchanged in another layer) for nonzero θ. Besides, the exact derivation manifests that t 1 (t 2 ) only grows as θ decreases (increases) (Supplementary Note 2). The off-diagonal function F = {f 11 , f 12 ; f 21 , f 22 } represents the spatially modulated interlayer hoppings, which can be obtained analytically according to Fig. 3 and single depends on r o under a given θ, as described in Supplementary Note 2. By diagonalizing Eq. (4), that is, E(Γ) = PH(Γ)P −1 (P is an invertible matrix), the spatial potential V(r o ) is given by the function fmin E ð ðΓÞ; maxðEðΓÞg, which is related to the energies of S 1 and S 2 in distorted lattices. In Figs. 4a, b, we show a specific cross section P1P2 (with length ffiffi ffi 3 p a M ) for θ c = 0 or π, where θ = 6.01 ∘ and h ¼ a= ffiffi ffi 3 p . The negative m* matches with S 1 and has positive potential energies, while the positive m* matches with S 2 and has negative potential energies. One sees that potential exhibits local valley (peak) characteristic for positive (negative) m*. The potential difference between AA and AB/BA stacked lattices always holds making the central AA stacked lattice become the global extrema of potential, which supports a 2D potential well of finite depth.
Consider the isotropy distortion approximation in the vicinity of central AA stacked region. The system can be regarded as the evolution of a spinless particle with effective mass m* in a given V(r o ) potential well. We describe this process using the time-independent Schrödinger-like equation with eigenstates Ψ, given by The solutions of Eq. (5) are shown in Fig. 4c. Discrete energy levels correspond to different orders of Ψ manifesting the arrangement of s, p x,y , ..., p x,y , s states from lowest to highest energies. The first half of these states (E < 0) is composed of S 2 with C 6 symmetry, while the second half (E > 0) is composed of S 1 with C 3 symmetry. At lowest and highest energies, s states isolated from the continuous bulk energy spectrum exhibit ideal confinement, which can be understood from the confining V(r o ) induced by intrinsic spatial hopping modulations. To further demonstrate the properties of general periodic twisted bilayer systems, we calculate band structures of moiré superlattices using the TBM with the hopping function of Fig. 2a. A representative result for θ = 6.01 ∘ and h ¼ a= ffiffi ffi 3 p is plotted in Fig. 5a. Four subbands (red curves) near the zero energy for spinless particles are fully consistent with typical moiré bands, corresponding to the divergent DOS, see the right panel of Fig. 5a. Whereas for the lowest and highest energies, superflat bands (blue curves) emerge in isolation accompanied by extremely large DOS, labeled as Ξ − and Ξ + . Figure 5b shows typical eigenstates at Γ S point of Ξ − , Ξ + and their adjacent bands. Ξ − (A) and Ξ + (F) correspond to s states formed by S 2 and S 1 , respectively. The eigenstates for E < 0 (A-C) and E > 0 (D-F) cases are consistent with the solution of the above continuous PWM in Fig. 4c. We further study the energies of Ξ − and Ξ + with different h and θ both in TBM and PWM, as displayed in Fig. 5c. Since such superflat bands are constrained by the potential of AA stacked lattices, i.e., E AA Γ , the energies of Ξ − and Ξ + vary exponentially with h in a wide range of θ. As h → +∞, the energies of Ξ − and Ξ + tend to −3t 0 and 3t 0 , respectively, merging into the bulk energy spectrum progressively (Supplementary Note 3).

Nanophotonic implementation
To realize superflat bands and the corresponding localized states in nanophotonic systems, we propose a twisted bilayer photonic crystal (PC) composed of an air layer . Two superflat bands (Ξ − and Ξ + ) are labeled in blue located at lowest and highest energies with extremely sharp DOS. Besides, typical moiré bands labeled in red appear around the zero energy with divergent DOS. b Various eigenstates for Ξ − and Ξ + bands and adjacent bands, which are arranged from the lowest to highest energy, i.e., A-F, corresponding to s, p x,y , ..., p x,y and s states. c The energies of superflat bands that vary with h around a= ffiffi ffi 3 p . Different cases of θ (i.e., 3.48 ∘ , 3.89 ∘ , 4.41 ∘ , 5.09 ∘ , and 6.01 ∘ ) have also been represented in different colors. Faint circles and dark dotted lines correspond to the results in TBM and PWM, respectively, while faint gray dotted lines indicate the bulk energies of AA stacked lattices for infinite h and two twisted PC slabs, as shown in Fig. 6a. Single PC slab has a C 6v lattice with lattice constant a Si = 1.5 μm filled with air, where the sublattices are composed of silicon triangular prisms (refractive index n Si = 3.46) with sidelength l Si = 0.35a Si and height h Si = 0.5a Si (see the left inset of Fig. 6a). The air layer with a thickness of d Si = 0.2a Si is sandwiched between two twisted PC slabs (see the right inset of Fig. 6a). The entire structure is embedded in perfect metal in the stacking direction forming a conservative system (here the transverse magnetic (TM) polarization is considered). We also provide the design under open systems as support (Supplementary Note 4).
Owing to the long-wavelength limit of dielectric PCs, the lowest band with linear dispersion near Γ point exhibits a fixed lowest frequency 0 leading to the absence of superflat bands with C 6 symmetric states (Supplementary Note 5) 44 . So we only present the case possessing C 3 symmetric states. Figure 6b shows band structures near Γ point for AA and AB/BA stacked PCs with given parameters in Fig. 6a. Because the electromagnetic fields are concentrated in carefully designed triangular prisms, this system can match well with TBMs 44 . The C 3 symmetric eigenstates of these two bands preserve particular frequency difference ensuring that the states located in AA stacked lattices is isolated from bulk spectra of twisted bilayer PCs (see the insets of Fig. 6b). Then, we calculate the band structure of twisted bilayer PCs with twist angle 6.01 ∘ , as plotted in Fig. 6c. The superflat band (blue) is observed at the frequency 116.3 THz, describing well-confined s states with C 3 symmetry, as shown in the top panel of Fig. 6d. Adjacent bands exhibit multipole states of moiŕe superlattices accompanied by worse localization capabilities. For example, p x,y states form crossed and nonflat bands, see Fig. 6c and the middle and bottom panels of Fig. 6d.
Note that such a design process exactly focuses on a single mode of the triangular prism (e.g., the fundamental mode above, which is therefore located in several lower bands). Despite the robustness of localized states, the interaction of different order modes of the triangular prism may merge the superflat bands into upper adjacent bands, which should be avoided when setting essential parameters of the system (Supplementary Note 6).

Discussion
The intrinsic superflat bands in our work have the property of isolated energy spectra without mode hybridization between different bands, so that the corresponding eigenstates have a clear and highly symmetrical phase distribution, as shown in Figs. 4c and 5c  localized eigenstates are almost insensitive to periodic moiré superlattice boundaries, which is understood as the origin of superflat bands and can be described by the PWM. The carried C 3 and C 6 symmetries distinguished from moiré flat bands formed by the four-band reconstruction (moiré bands) near the zero energy have not been fully discussed before [9][10][11]26 . Recently, we notice that a displacement electric field is applied in specific twisted bilayer systems (e.g., graphene and boron nitride heterostructure) to study the valley topology of moiré bands 45,46 .
In our system, this is equivalent to yielding a nonzero |ϵ| with distinct signs for two layers. The energies of superflat bands will be corrected corresponding to a shift g(|ϵ|), where g(|ϵ|) ≥ 0 and grows as |ϵ| increases, see the details given in Supplementary Note 7. Apart from that, nonzero |ϵ| cannot affect the presence of superflat bands and localized states.
In conclusion, combining theoretical PWM analysis and TBM calculation, we have demonstrated a class of superflat bands with C 6 and C 3 symmetric states for small twists in general twisted bilayer systems. The dislocated lattices formed by the systematic hopping modulation create macroscopic effective potential wells centered around the AA stacked region, leading to the wellconfined states described by the PWM. We also mimic these two effects in nanophotonic systems displaying the unique electromagnetic wave confinement. Notably, superflat bands and the corresponding localized states can be realized for continuous twist angles (distinct from the discrete set of twist angles in magic-angle physics), showing a class of generalized effects of twisted bilayer systems distinguished from the fragile topology. The concept of generalized localized states may inspire a shortcut technology for generating zero-dimensional localization, avoiding complex boundary splicing of (higher-order) topological insulators, which will greatly benefit the wave trapping and manipulation. The frequencies of superflat bands and the configurations of localized states can be adjusted by twist angles, and this offers an advanced platform for reconfigurable devices. Our results can be extended to photonics [47][48][49] , phononics, and mechanical waves, where ideal transport can be realized for integrated chips in information technologies.

Nanophotonic simulation
Numerical simulations for nanophotonic systems in this work are all performed using the 3D electromagnetic module of commercial finite-element simulation software (COMSOL MULTIPHYSICS). In solving the eigenvalues and eigenstates of AA, AB/BA, and moiré lattices in our silicon-air platform, the calculation regions are selected as hexagonal unit cells with side lengths a ffiffi , respectively, with a being 1.5um and θ being 6.01 ∘ . Bottom and top boundaries along the stacking direction are set as perfect electric conductors. So only the transverse magnetic (TM) polarization is considered for the data in Fig. 6b and c, i.e., E z . Whereas 2D periodic directions satisfy Bloch's theorem E z (r + R) = e ik⋅R E z (r), where R is a real space lattice vector.