Magic configurations in Moir\'e Superlattice of Bilayer Photonic crystal: Almost-Perfect Flatbands and Unconventional Localization

We investigate the physics of photonic band structures of the moir\'e patterns that emerged when overlapping two uni-dimensional (1D) photonic crystal slabs with mismatched periods. The band structure of our system is a result of the interplay between intra-layer and inter-layer coupling mechanisms, which can be fine-tuned via the distance separating the two layers. We derive an effective Hamiltonian that captures the essential physics of the system and reproduces all numerical simulations of electromagnetic solutions with high accuracy. Most interestingly, \textit{magic distances} corresponding to the emergence of photonic flatbands within the whole Brillouin zone of the moir\'e superlattice are observed. We demonstrate that these flatband modes are tightly localized within a moir\'e period. Moreover, we suggest a single-band tight-binding model that describes the moir\'e minibands, of which the tunnelling rate can be continuously tuned via the inter-layer strength. Our results show that the band structure of bilayer photonic moir\'e can be engineered in the same fashion as the electronic/excitonic counterparts. It would pave the way to study many-body physics at photonic moir\'e flatbands and novel optoelectronic devices.

Moiré structures have been of central interest in fundamental physics during the last few years. The most important milestone is the discovery of flatbands in the moiré patterns emerged when two graphene layers are overlapped at certain at magic twisted angles [1][2][3], leading to nonconventional superconductivity [4][5][6] and strongly correlating insulator states with nontrivial-topology [7,8]. Motivated by the electronic magic angles, photonic moiré has attracted tremendous research in light of shaping novel optical phenomena. Hu et al. have demonstrated [9,10] the topological transition of photonic dispersion in twisted 2D materials. However, the operating wavelength in these pioneering works are much larger than the moiré period, thus dispersion engineering is based on the anisotropy of an effective medium rather than the microscopic moiré pattern. On the other hand, Ye's group has recently reported on the realization of 2D photonic moiré superlattice [11]. Nevertheless, this work only focused on light scattering through the moiré pattern, but the lattice is on the same plane, and there is no bi-layer, neither twisting concepts. Most recently, numerical [12] and tight-binding [13] method have been proposed to investigate twisted bilayer photonic crystal slabs. In particular, Dong et al. has showed that local flatband would be achieved [13] in twisted bilayer photonic crystal at small twisted angle..
In this work, we report on a theoretical study of photonic band structures in moiré patterns that emerged when two mismatched 1D subwavelength photonic crystal slabs are overlapped. The essential physics of the system can be captured by an effective four-component Hamiltonian. Accompanying the analytical theory, numerical electromagnetic simulations are performed with a case study of silicon structures operating at telecom wavelength. The obtained band structure are resulted from an interplay between intra-layer and inter-layer coupling mechanisms which is tuned via the distance separating the two layers. Importantly, magic distances corresponding to the emergence of photonic flatbands within the whole Brillouin zone are demonstrated. The minibands of moiré superlattice can be described by a single-band tight-binding model with Wannier functions tightly confined within a moiré period. The tunnelling rate of light between nearest neighbor Wannier states is continuously modulated by the inter-layer distance and vanished at magic distance, leading to flatband formation and photonic localization. Despite its simplicity, this 1D setup captures much interesting physics of moiré systems of twisted two-dimensional materials. Our findings suggest that moiré photonic is a promising strategy to engineer photonic bandstructure for fundamental research and optoelectronic devices.
Our system consists of two 1D photonic crystal slabs which are two subwavelength high refractive index contrast gratings (Fig 1). These gratings have the same subwavelength thickness h and filling fraction κ and are separated by only a subwavelength distance L. Their periods a 1 and a 2 are slightly different but satisfying the commensurate condition arXiv:2104.12774v3 [physics.optics] 11 Jun 2021 a 1 /a 2 = N/(N + 1) for a natural number N . The period of the superlattice is given by Λ = (N + 1)a 1 = N a 2 , consisting of N + 1 periods of the upper grating and N period of the lower one. In the regime of N 1, a semi-continuous approach can be implemented: the two gratings are almost identical and the moiré pattern corresponds to a continuous shifting function δ(x) of the upper grating with respect to the lower grating, given by δ(0 ≤ x ≤ Λ) = x/N . The shifting δ sweeps an amount a 0 = (a 1 + a 2 )/2 when x varies across a moiré period. In other word, the moiré superlattice is obtained from the bilayer lattice by introducing a slight period mismatch: the period of the upper grating is shrunken from a 0 to a 1 and the period of the lower one is stretched from a 0 to a 2 . This configuration leads to a modulated relative displacement δ(x) with respect to the coordinate x. Two special configurations of δ/a 0 = 0 and 0.5 are referred to as AAand AB-stackings, resembling the terminology in Bilayer Graphene structure [14]. The moiré pattern is a period of a superlattice made of bilayer structures varying continuously from AA stacking to AB stacking. The period mismatch leads to a Brillouin zone mismatch and the size of the mini Brillouin zone K M is given by K M = K 1 − K 2 , where K 1 = 2π/a 1 and K 2 = 2π/a 2 .
In our perturbation approach, the dispersion characteristic of the moiré superlattice is derived from two coupling mechanisms among forward (φ 1+ , φ 2+ ) and backward (φ 1− , φ 2− ) fundamental guided waves of the two noncorrugated slabs with effective refractive index: i) Intra-layer coupling due to the diffractive processes [15] between counter-propagating waves from the same layer. ii) Inter-layer coupling via evanescence between co-propagating waves from separated layers. Using (φ 1+ , φ 1− , φ 2+ , φ 2− ) as basis, eigenmodes of the system are described by the following Hamiltonian (detailed derivation is given in the Supplemental Material): Here U 1,2 are the intra-layer coupling rates and V is the inter-layer one; v and ω 1,2 are the group velocity and offset energy of the guided waves at the Brillouin zone edge for each grating. A slight difference of values of the offset pulsation and the intra-layer coupling strength for each grating are due to the period mismatch, with ω 1 ≈ ω 2 ≈ ω 0 and U 1 ≈ U 2 ≈ U where ω 0 and U are the offset pulsation and the intra-layer coupling strength in the grating of period a 0 .
The energy-momentum dispersion is simulated numerically using Rigorous Coupled-Wave Analysis (RCWA) method [16][17][18]. The numerical results corresponding to N = 13 when increasing the separation distance L are presented in figures 2a-e. When L is comparable to a 0 , the band structure is simply the folding of single layer dispersions (Fig 2a). It suggests that the inter-layer coupling mechanism is negligible with respect to the intra-layer ones (i.e. V U ) for L a 0 .
In this configuration, a bandgap, purely due to the intra-layer coupling mechanism, is observed (Fig 2a). In analogy to semiconductor terminology, we refer to these upper/lower bands as conduction-like/valence-like. When L a 0 , the band hybridization due to the inter-layer coupling results in the formation of a pair of particle-hole minibands, referred to as electron-like/hole-like moiré band (Figs 2b-e). These two bands emerge within the bandgap of uncoupled layers and are well isolated from the conduction/valence-like continuum. In the following, we will pay particular attention to the behavior of these two bands when tuning the inter-layer interaction.
One may note that with the choice of a 0 = 300 nm, the spectral range of the these band is in the telecom (i.e.∼ 1.5µm). Intriguingly, there are some specific values of L at which the bandwidth of these bands becomes almost zero, and these moiré bands are nearly perfectly flat. Figures 2c and 2d depict the band structures with flat hole-like moiré band, and almost-flat electron-like band. Inspired by the analogy with the appearance of flatbands at magic angles in twisted bilayer graphene [2], we called these values magic distances.
The moiré band structure is calculated using the Hamiltonian model given by Eq. (1), taking v, U , ω 0 and V as input parameters. These parameters are retrieved from the simulation of single and bilayer lattice [19,20]. Figures 2f-j depict the band structures obtained by analytical calculations. These results reproduce quantitatively the numerical results presented in Figs 2a-e, showing the emergence of moiré states within the bandgap and their flattening at magic distances. Noticeably, there is a slight difference between simulation and analytical results: the RCWA suggest that the flattening of the electron-like band always takes place at a slightly smaller distance L than the one of the hole-like band, while the Hamiltonian model predicts that both bands become flat almost simultaneously.
The global spectral bandwidth, defined as ∆ω = max q (ω) − min q (ω), is used as the figure of merit to evaluate the flatness of moiré minibands. Figures 3a-c depict the global spectral bandwidth of the hole-like moiré band if different moiré superlattice (N = 9, 13 and 19) when scanning L. These results confirm the existence of magic distances, corresponding to the bandwidth vanishings. All of the analytical calculations are obtained with the same set of parameters that are previously presented. We highlight that the Hamiltonian model provides almost perfectly both the number of magic distances and its values.
For each moiré superlattice (i.e. a given N ), our design exhibits two adjustable parameters: i) The distance L for tuning the inter-layer coupling V (V = V 0 when L = 0 and decreasing exponentially when increasing L [20]) ; ii) The filling fraction κ, defined in Fig.1, for tuning the intra-layer coupling U (U = 0 when κ = 1 and increasing when decreasing κ [20]). Up to now, we have been investigating flatband emergence by scanning L while fixing κ = 0.8 (i.e. U = U 0 ). However, the direct parameters of the Hamiltonian (1) are U , V and N (from K M ). Thus a complete picture of magic configuration is captured when varying both V /U (i.e. competition Figure 2. (a-e) Simulated band structures corresponding to different L values. The design for the simulation uses silicon (n = 3.54) as the grating material, with h = 180 nm, κ = 0.8, a0 = (a1 + a2)/2 = 300 nm and N = 13. The photonic modes of uncoupled gratings are located below the light-line and the inter-layer coupling mechanisms, if not strong enough, would not be able to make these modes accessible for RCWA simulations. To solve this, a double period perturbation of 5% is implemented for the design of each grating. The unit-cell in RCWA simulation consists of two moiré patterns: One is shrunken to 0.95Λ, and the other one is dilated to 1.05Λ.  between inter versus intralayer coupling) and N (i.e. moiré pattern). Figure 3d presents the global bandwidth when scanning N and V /U within a reasonable range [21]. The observed "resonant dips" correspond to different magic config-urations. Dimensional analysis of Hamiltonian (1) suggests that our system is driven by two dimensionless ratios V /U , and U/K M ∼ N U [20]. Indeed, fitting the resonances of Fig.3d by a power law, we obtain a very simple empirical relation between there two dimensionless parameters: with the γ ≈ −1.42, η ≈ 12U 0 , and m is the counting order of the magic configuration. We note that N is the "moiré parameter" in our system and playing the same role as the twist angle in twisted bilayer graphene (each value of moiré parameter defines a moiré pattern) [22,23]. Therefore, the good metric for magic configurations is the magic number N m , and Eq.
(2) provides the design rule to achieve them. The analogy and similitude between this law and the one for magic angles in twisted bi-layer graphene [2] are striking and we expect an appealing interpretation for this simple relation.
Knowing that flatband states would give rise to an unconventional localization regime [11,27,28], we now investigate the localization of light at magic configurations. A closer look at the two moiré bands in Figs. 2 reveals that their dispersion characteristic are nearly single harmonic functions with the dominance of the first Fourier component with respect to higher-orders. Consequently, this suggests  Fig. 2. (e) The energy splitting between the two bound states as a function of L. Red circles are results from FDTD simulations, and the solid black line is the result from the effective Hamiltonian. For the analytical calculation, the boundary condition is chosen so that outside of the moiré molecule is bilayer structure of AB sites, and the bound states are calculated by the transfer matrix method [25,26]. For FDTD simulations, the structure only consists of two moiré cells.
that each moiré band may be described by a textbook singleband tight-binding model with only a few nearest neighbour couplings taken into account. It is of interest to compute the Wannier functions for the band under consideration since Wannier functions are the natural basis for the tight-binding model [29]. Figure 4(a) depicts the result of this calculation when scanning the ratio V /U , showing that more than 94% of the Wannier density is located within a single moiré cell. Such a concentration confirms the use of this Wannier function as a pseudo-orbital wave function for the tight-binding model with nearest neighbour couplings. However, it is important to stress that the high concentration of the Wannier function is not necessarily related to flatband formations. Yet, the physics of the moiré bands can be captured quite well by a simple tight-binding scheme in the Wannier basis. In this scenario, the moiré superlattice engenders a periodic potential landscape with minima at AA sites. Trapped photons in the Wannier states can tunnel to the nearest neibour ones with tunnelling rate J to form moiré bands of bandwidth 2 |J|. As consequence, when the couple (N U, V /U ) satisfies the Eq. (2) of magic configurations, the only way to obtain dispersionless bands is that the effective tunneling rate J becomes zero. This leads to the tightly localization of light within a single moiré cell at magic configurations. The compact localized states [30] of our localization is simply the Wannier function. We notice a resemblance of the flatband emergence in our system compared to the one in twisted bilayer graphene system [2,31,32]: both correspond to the good localization at the AA sites.
Keeping in mind the ability to localize light to a moiré period with high quality (albeit non-perfect), we investigate a much simpler problem: a "diatomic molecule" made of two moiré cells (Fig.4b). Figures 4c,d depicts the field distribution of the hole-like bound states with even ( Fig. 4c) and odd ( Fig. 4d) parity regarding the lateral mirror symmetry. The energy splitting when scanning the distance L is presented in Fig. 4e. Again, the results from the analytical model and numerical simulations show a very good agreement. Notably, these results demonstrate the crossing of these bound states exactly at the magic distances of the moiré superlattice from Fig. 3b. Consequently, it supports that the tunnelling rate J changes sign when scanning L across a magic distance value and vanishes when L takes a magic distance value.
In conclusion, we have investigated theoretically the 1D moiré superlattice of bilayer photonic crystal. All of analytical results derived from a simple effective Hamiltonian are in good agreement with numerical simulations, showing the emergence of flatband at magic configurations. The conditions for flatbands unify to a nontrivial relation between the counting order of the magic configuration and the magic number, given by N m ∼ m. The physics of the moiré minibands is captured by a simple tight-binding model, resulting in localization of photonic states within a single moiré period at flatband configurations when the tunnelling rate vanishes accidentally. As fundamental perspective, the implementation of nonlinearity via Kerr nonlinearity [33] or exciton-polariton platform [34], would pave the way to investigate the strongly correlated bosonic flatband physics [35-37] with intriguing phases of 1D matters [38]. For applications in optoelectronic devices, the design in this work uses silicon as dielectric material, operating in the telecom range with feasible fabrication [19,39,40], and is transferable to 1D integrated optics. The high sensitivity of dispersion band structure to the refractive index of surrounding medium (which determines the parameter U ) and spacing medium (which determines the parameter V ) can be harnessed for applications in sensing. Furthermore, the localization of light within the moiré period also suggests a unique way to engineer lattice of resonators of a high-quality factor for a phase-locked micro-laser array or high Purcell factor for tailoring spontaneous emission of quantum emitters. Another realization scheme is with dualcore fiber Bragg gratings [41, 42] to study soliton physics arising from photonic nonlinearity which will be greatly en-  In this section, we provide the detailed derivation of the effective Hamiltonian in the main text.
Hamiltonian of a single grating wave-guide Wave function a single grating wave-guide In perturbation theory, the eigenmodes in photonic grating are constituted by the coupling between forward ϕ + (k ≥ 0) and backward ϕ − (k ≤ 0) propagating waves of the non-corrugated waveguide of effective refractive index (see Fig S1a). Here the wave-function ϕ corresponds to the electric field E y for TE modes, and the magnetic field H y for TM modes. The dispersion characteristic ω + (k ≥ 0) and ω − (k ≤ 0), ω + (k) = ω − (−k), of these guided modes lies below the light-line (see Fig S1b) and are obtained by solving Maxwell equations of planar waveguide with effective refractive index. We can extend the definition of positive and negative wavefunctions for any k value by replacing ϕ ± (k) by Φ ± (k), given by: where Θ is the Heaviside function, With such definition, the spatial wave-function Φ ± (x) of positive and negative modes is obtained by the Fourier transform of Φ ± (k): With a spatial period a, the reciprocal lattice vector is given by K 0 = 2π a0 . High symmetry points in the momentum space are at wavevectors lK0 2 with l ∈ Z. A given odd(even) value of l corresponds to a X(Γ) point of the BZs. The effective wave-functions of positive (negative mode) near the high symmetry point lK0 2 (-lK0 2 ) are defined by: and We verify easily the relation between Φ ± (x) and Φ l,± (x), given by: Since band structures are mostly studied in the vicinity of a high symmetry point of the BZs, the most appropriate basis in real space and momentum space given by: Note that due to the fact that positive mode has positive wave-vectors and the negative mode has negative wave-vectors, only l ∈ N * appears in the definitions (S6). In the vicinity |q| K0 4 of high symmetry points (the blue points in Fig S1b) in momentum space, these relations can be approximated by We have the effective free Hamiltonian density in momentum space H free (l) (q) near the high symmetry points in the momentum space Diffractive coupling between counter-propagating waves Due to grating, the positive and the negative modes couple with each other via diffractive coupling where the diffractive coupling function U D (x) is periodic with period a: where U l = U −l because of the C 2 symmetry (reflection x → −x) of the grating . Due to the diffractive coupling, effectively, the positive mode couple with the negative mode that is shifted by lK 0 in the momentum space. Vice versa, one can think of the diffractive coupling is the negative mode couples with the positive mode that is shifted by −lK 0 in the momentum space.
Note that since positive mode has positive wave-vectors and the negative mode has negative wave-vectors, only l ∈ N * appears in the summation of Eq. (S11) and n runs from −l to l due to momentum conservation. However, the effective coupling becomes important when the energies of positive and negative bands are approximately identical, which corresponds to n = 0. Hence we rewrite the diffractive Hamiltonian as: Combining the (S8) and the diffractive coupling (S12), we can derive the effective Hamiltonian near the high symmetry point in the momentum basis From now on, we will concentrate on the high symmetry point corresponds to l = 1. We then replace ω 01 → ω 0 , U 1 → U and v 1 → v, thus: In the subsequent sections ,we will obmit the l indices and implicitly use Φ ± as Φ 1,± in (S4).
Effective Hamiltonian of bilayer To understand the inter-layer coupling mechanisms, an intuitive and informative example is the configuration of bilayer photonic lattice, referred to as the "fish-bone" structure in Ref [19]. Such a configuration consists of two identical gratings, one on top of the other with a relative displacement δ 0 (Fig S2). Two special configurations of δ 0 /a 0 = 0 and 0.5 are respectively the equivalent of AA and AB stackings in Bilayer Graphene structure [14]. We use the notation system with the implementation of index (1) and (2) to distinguish the upper and lower layer. We consider the basis made of effective wave-functions near the crossing point of the positive and the negative bands of each layer . (S15) Similar to the case of single layer, the Hamiltonian densities of uncoupled layers in these basis are: The evanescent coupling of the bilayer configuration is [44] In the regime in L a, we can assume that V f-b (x − y) = V.δ(x − y − δ 0 ). If we only consider the effective modes near the symmetry point corresonds to l = 1, we rewrite the inter-layer coupling Hamiltonian as: We then replace ∂ y → ∂ x and y → x − δ 0 in equation (S15) . The effective basis when working with both layers is given by: for real space and momentum space respectively. The matrix representation of inter-layer coupling Hamiltonian of Eq.(S18) in real space is written as: with the interlayer coupling matrix The bilayer Hamiltonian consists of the Hamiltonian of uncoupled layers and the inter-layer coupling Hamiltonian. Using effective Hamiltonians (S16) and the interlayer coupling (S20), we obtain the effective Hamiltonian for the bilayer system: Another form of the bilayer Hamiltonian in momentum space is reported in Ref. [19]. One can show that the two bilayer Hamiltonians are equivalent using a simple transformation of the basis.

Hamiltonian of uncoupled layers
We now consider a moiré bilayer of parameters as discussed in the main text (see Fig. S3). With such geometrical design, the Hamiltonian of the uncoupled layers from moiré configuration has the same form as the ones of uncoupled layers as in the bilayer configuration. The only difference to the bilayer configuration is the mismatch between BZ-sizes of the two layers (K 1 = 2π/a 1 for the upper layer, and K 2 = 2π/a 2 for the lower layer). The decomposition of wavefunctions corresponding to positive and negative modes is given by: l,± (q)e iqy .
Since the mismatch between BZ-sizes K M K 1 , K 2 , we expect the interlayer coupling to play an important role when the upper and lower modes are near the symmetry point with the same index l. We consider the effective theory near the symmetry point l = 1, and omit the l index by implicitly use Φ . (S25) The Hamiltonian densities of uncoupled layers in these basis are: The parameters of the Hamiltonians (S26) are determined from the simulation and experiment fitting for a single-layer unidimensional photonic crystal slab.

Hamiltonian of inter-layer coupling: moiré configuration
Co-propagating waves of the same momentum but from different layers are coupled via evanescent coupling. The evanescent mechanism is written in term of the coupling Hamiltonian as When L a, we can assume that V(x − y) = V δ(x − y − δ 0 ) where δ 0 is the offset shift between the two layers. Moreover, as discussed in the main text, the value of δ 0 is not relevant for the moiré structure, and we can assume it to be zero. Considering the effective model near the symmetry points corresponding to l = 1. We rewrite the inter-coupling Hamiltonian (S27) as the coupling of effective basis Φ Some remarks are in order. We now understand the origin of the spatial dependent phase shift φ(x) in Eq. (??) in the main text by looking at the expansions (S23) and (S24). Due to the mismatch between BZ-sizes, there is a different phase between the upper and the lower modes near the symmetry points corresponding to the same m. We then choose an effective basis when working with both layers We then replace ∂ y → ∂ x in Eq. (S26) and obtain the matrix representation of the effective Hamiltonian of Eq.(S28) in the effective basis (S29): .
(S30) the interlay coupling matrix Figure S4. (a)Phase-matching condition (conservation of momentum) for inter-layer coupling between co-propagating waves. (b) Inter-layer coupling mechanism in momentum space between different moiré B: Modes in the upper(lower) layer with Bloch momentum q couple to modes in the lower(upper) layer with Bloch momentum q − K M 2 and q + K M 2 . Each moiré BZ is indicated by its index p.
The difference of period would lead to a slight difference of values of the offset ω 0 and the intra-layer coupling strength U for each grating and a small modification of V with respect to the case of Bilayer lattice. However, since ω 0 U, V , in the first approximation, only ω 0 varies when switching from upper to lower layer.
The decomposition (S31) shows two types of inter-layer coupling in momentum space: • The positive mode with effective momentum q in the upper layer will couple to the positive mode with effective momentum q + K M 2 in the lower layer via T 1 .
• The negative mode with effective momentum q in the upper layer will couple to the negative mode with effective momentum q − K M 2 in the lower layer via T 2 .
We demonstrate this coupling mechanism in the momentum space explicitly in Fig S4b. This situation is similar to the inter-layer coupling model suggested by Bistrizer and Mac Donald in twisted bilayer graphene [1]; the only difference is that in twisted bilayer graphene, there are three couplings T 1 , T 2 and T 3 corresponding to three momentum shifts instead of just two.

A change of basis:
The effective basis (S29) was chosen in the same manner as in the twisted bilayer graphene literature [2]. Consequently, the Hamiltonian (S30) shares the same pattern as the Hamiltonian derived by Bistritzer and MacDonald in Ref [1,2] as expected. Notice that in the effective basis (S29), the origins of the effective momenta are different. The wave-function in the coordinate space is given by The wave-function of the electromagnetic wave near the vicinity of X point on the upper layer and lower layer can be read off from (S32) as One can use (S33) solved from effective Hamiltonian (S30) to compare directly with the electromagnetic wave in coordinate space of simulations and experiments. However, it is helpful to introduce another effective basis such that the wavefunction in the coordinate space is which impliesΦ (2) We are able to rewrite the moiré Hamiltonian (S30) in the new effective basis Since K 1 = (N + 1)K M , the momentum of the effective basis Φ where the characteristic wavevector is q 0 = K M /2 with K M is the moiré wavevector. Here ⊗ denotes the Kronecker product, σ x,y,z are Pauli matrices defined by σ x = 0 1 1 0 , σ y = 0 −i i 0 , σ z = 1 0 0 −1 , and σ ± = (σ x ± iσ y )/2. The last term in Hamiltonian (S38) only leads to minor quantitative corrections; for qualitative analysis, one can set ∆ U = 0. We see then that the equation (S38) is characterised by parameters (U, ∆, q 0 , V ). All of these quantities have the same dimension of energy (since v = 1). One can effectively set one of them, e.g., U , to be the unit. Moreover, when we specialise to the particular realisation of the effective Hamiltonian (S38), as in Appendix , we see that the parameters ∆ = ω (1) − ω (2) , ∆ U = U (1) − U (2) and q 0 = K M 2 are in fact physically dependent through the straining parameter in the system. In this case, we therefore only have three independent physical parameters (U, q 0 , V ). The model is specified by two dimensionless ratios between the independent parameters.

Periodicity and the Bloch Hamiltonian
It is perhaps surprising when one notices that the Hamiltonian (S38) seem to be periodic with the double supercell period 2π/q 0 = 2Λ, which we refer to as apparent period. Accordingly, naively solving these Hamiltonian one obtains a band structure with the apparent Brillouin zone of size K M /2. The Hamiltonian is in fact of higher translational symmetry. Indeed, let T Λ be the translation operator of one moiré period. Then one can easily verify that [45] the Hamiltonian is invariant under the generalized translational operator T Λ (σ z ⊗ I). The operator T Λ (σ z ⊗ I) generates the commutative group of generalised translational operators, under which the Hamiltonian is invariant. This shows that the actual period of the system is, not surprisingly, the moiré period Λ.
With the apparent period of 2Λ, the Bloch theorem states that we can assume the eigenstate of the Hamiltonian (S38) to be of the form where q is the moiré Bloch vector, −q 0 /2 ≤ q ≤ +q 0 /2 and the four-spinor u q (x) is periodic with the apparent period 2Λ. This leads to the Bloch Hamiltonian for u q (x), This Hamiltonian is to be solved for eigenvalues E q with periodic eigenstates u q (x), where the latter is also denoted by u E,q (x) when the explicit energy value is necessary for the clarity. The periodicity of the Bloch wavefunction u q (x) allows for the solution of the eigenvalue problem to be found through Fourier expansion.
It is important to emphasize again that when using the apparent period 2Λ of the Hamiltonian to calculate the band structure, the Bloch momentum q in equation (S40) is folded within [−q 0 /2, q 0 /2]. In order to unfold the band to the full moiré Brillouin zone [−q 0 , q 0 ], one simply solves the Bloch Hamiltonian for q in the full moiré Brillouin zone, but maintains only solutions that satisfy the generalised Bloch theorem T Λ (σ z ⊗ I)Ψ(x) = e iqx Ψ(x). In this way, the unfolded band structure such as in Fig. 2 can be obtained.

Symmetry analysis
Since the moiré system has spatial refection and time-reversal symmetries, one expects that the Hamiltonian (S38) also carries these symmetries. This is indeed the case: • Spatial reflection: Consider the reflection along the x-axis. Let P denote the pure spatial coordinate reflection operator.
Since the reflection of the x-axis also changes the signs of the momenta within each chain, it also exchanges the two basis wavefunctions chosen is Section . Therefore one can expect that the full reflection operator to be (I ⊗ σ x )P . One can easily verify that the Hamiltonian (S38) is indeed invariant under this full reflection operator (I ⊗ σ x )P . As also expected, the spatial reflection (I ⊗ σ x )P brings the Bloch Hamiltonian H q in Eq. (S40), to H −q , implying the energy bands are symmetric under reflecting the Bloch wavevector, E q = E −q , and the Bloch wave functions obey u E,q (x) = u E,−q (−x). Being even and periodic with respect to the moiré wavevector K M = 2q 0 , an energy band E(q) is completely described by Fourier coefficients f p = 1 2q0 +q0 −q0 dq cos(pπ/q 0 )E(q). Fig. S5 presents the Fourier components f p of the flat band as indicated in Fig. 2 in the main text, but with slightly different dimensionless parameters as indicated in the caption. It shows that the first coefficients of the Fourier dominate over higher Fourier components, suggesting in an effective tight-binding model, the nearest neighbour coupling dominates. Importantly, higher Fourier coefficients, although small, do not vanish when the first coefficient vanishes (near the flat band). This suggests that the band, although becomes highly flat at the magic coupling, is not perfectly flat.

Probability density distribution of the Bloch wave functions and Wannier functions
To study the Bloch wave function near the flat transition, we compute the density (S41) Figure S6 (left) demonstrates the probability densities of the Bloch wavefunctions with varying coupling V across a flat transition. It is important to notice that while the Bloch wave functions tend to concentrate within a moiré period, they do not vanish anywhere (also when the band is flat). In particular, there is no qualitative change in the density of the Bloch wave function as the band is crossing the flat transition.
To consider the possibility of concentrating light in the moiré lattice, we compute the Wannier functions for the Bloch Hamiltonian (S40). The computation of the Wannier function requires fixing the arbitrary phase in the numerical solution of the eigenvectors of the Bloch Hamiltonian (S40). This is a known difficulty in computing Wannier functions with maximal localization [24]. Fortunately, in one-dimensional systems, there is a known gauge fixing procedure, the twisted parallel transport gauge, that allows for the computation of Wannier functions of maximal localization [24].
Upon fixing the twisted parallel transport gauge, the Wannier function is then obtained directly as Notice that the integral runs over the full moiré BZ [−q 0 , q 0 ], that is, twice as much of the apparent BZ [−q 0 /2, q 0 /2]. Figure S6 (right) plots the probability density of the Wannier wavefunction (S42). While being highly concentrated, see Fig S7, one should notice that the Wannier function extends beyond a single Moiré period. There is also no qualitative change in the density of the Wannier wave function as the band is crossing the flat transition.

Dynamical signature of flatbands
From the above analysis, it is clear that the concentration of the probability density of the Bloch wave function or the Wannier function is not the signature of the flat band. In fact, the very physical meaning of localization in this context is a dynamic one.
Suppose the system has a flatband u q (x), that is for some energy level E q = E 0 , independent of q. Then this is nothing but saying that Ψ q (x) = e iqx u q (x) are having the same energy E q = E 0 for all q. This means that given any wave function in momentum space v q , the wave packet is also an eigen-wavefunction with energy E 0 . As a consequence, the probability density Ψ(x) † Ψ(x) is unchanged overtime. This is true for any wavepacket v q in the Bloch momentum space, in particular the Wannier function (S42).

Finite systems: tunnelling and resonances, bound states
To understand better the nature of the flat bands, we compute the tunnelling and bound states of light in a finite number of moiré periods. This calculation can be carried out employing the (generalised) transfer matrix method [25], particularly adapted to the case of Dirac-like equations in Ref. [26].
To do so, we rewrite the eigvenvalue equation where H(x) is a 4 × 4 matrix given by All possible x-evolutions of the x-dynamical equation (S45) is described by the 4 × 4 x-evolution operator G E (x 2 , x 1 ), which is the solution of subject to the initial condition G E (x 2 , x 1 ) = I. The function G E (x 2 , x 1 ) summarises all information about the eigenwave function Ψ(x) corresponding to the eigenvalue E of the Hamiltonian H. Therefore it is a convenient way to relate different properties of H, such as the existence of extended states, transmission amplitudes, probability distribution, the density of states, etc. On the other hand, with well-developed methods for the ordinary differential equations (ODEs) [46], the computation of G E (x 2 , x 1 ) is relatively easy. One should, however, notice that the x-dynamics is non-hermitian and sometimes numerical instabilities have to be addressed.

Boundary condition and the computation of tunnelling rate
To investigate the tunnelling phenomena through the finite moiré structure between x 1 and x 2 , one has to consider the realisation of the asymptic area outside the moiré structure. For convenience, we choose this to be of the type of fishbone structure [19]; that is, fixing the phase in the coupling between the two chains in the Hamiltonian (S38) to be e ±iq0x1 (constant) for x ≤ x 1 , and e ±iq0x2 (constant) for x ≥ x 2 .
For the fixed phases e ±iq0x1 or e ±iq0x2 , the eigenstate of the Halmitonian (S38) can be easily solved, resulted in the fishbone band structure [19]. Plugging a plane-wave solution Ce ikx into the resulted Hamiltonian, one finds the fishbone eigenvalue equation, where φ = q 0 x 1 or φ = q 0 x 2 , which are here simply constants. Fixing the energy E, we are interested in solving this equation for k. The resulted equation is a generalised eigenvalue problem. In general, the obtained generalised eigenvalues k are complex.
To fix an ordering, we order the four (generalised) eigenvalues k according to their increasing phases, that is, the angles with respect to the real axis, computed counterclockwise. Let us consider the possible solutions of equation (S48). One sees that if k is a solution, k * is also a solution (time-reversal symmetry). Also, if k is a solution, −k is also a solution (spatial reflection symmetry). In general, one has 4 different wavevectors satisfying (S48). If one of the solution k is generically complex (i.e., not pure real or pure imaginary), then by acting with the time-reversal symmetry and reflection symmetry, one obtains all the other three solutions k * , −k, −k * , which are also generically complex. On the other hand, if one of the solution k is real, then the time-reversal symmetry and the reflection symmetry only give −k as another solution. There are then two possibilities: the other two solutions can also be real, or they must be purely imaginary.
To consider the tunelling phenomena, we are interested in the energy range of +∆ + U ≤ E ≤ −∆ + U (for ∆ < 0). Here for a fixed energy E, there are two real wavevectors ±k (with the convention k ≥ 0), corresponding to the phases of 0 and π. Two other modes are of pure imaginary wavevectors iκ (with the convention κ ≥ 0) corresponding to exponential decaying or exponential amplifying modes and phases of π/2 and 3π/2. By W we denote the matrix of which the columns are the corresponding eigenvectors (ordered such that phases of the eigenvalues increase, here must be 0,π/2,π and 3π/2). The general wavefunction depends on 4 amplitudes of these different solutions, a ± and b ± , explicitly given by Derivation of the band structure of the infinite system As an interesting side remark, we mention that the band structure of the system can also be computed from the generalised transfer matrix G(x 2 , x 1 ). To this end, we choose x 2 − x 1 to be an apparent period of the potential (twice as much of the moiré period), x 2 − x 1 = 2Λ. Then from the fact that Ψ(x 2 ) = G(x 2 , x 1 )Ψ(x 1 ) and the Bloch theorem Ψ(x 2 ) = e iq2Λ Ψ(x 1 ) we obtain det[G(x 2 , x 1 ) − e iq2Λ ] = 0. This allows one to compute the Bloch wavevector corresponding to the energy under consideration E. By selecting the real wave vector q, the band structure of the system can then be derived.

PARAMETER RETRIEVAL FOR THE EFFECTIVE HAMILTONIANS
The effective Hamiltonian of the moiré structure is determined by the energies ω (1,2) 0 , U (1,2) , V and the group velocity v. These values are retrieved from the dispersion characteristics of the single layer structure (for ω 0 ,U and v), and of the bilayer structure (for V ) which are obtained by RCWA simulations. In the following, we will discuss in details these parameter retrieval methods.
Parameter retrieval of single grating structure The dispersion characteristic of a single grating structure is easily calculated from Eq. (S14) in the main text. It consists of two bands of opposite curvature ± v 2 2U , with corresponding band edge energies given by ω 0 ± U . As a consequence, ω 0 and U are directly extracted from the energy of resonances at q = 0 of the RCWA simulations. Then knowing U , the group velocity v is extracted from the curvature of these resonance. As shown in Fig. S9b, the band structure which is calculated by the effective Hamiltonian using the retrieved parameters reproduce perfectly the simulated one.
With the retrieval method presented above, we can explore the dependence of U , ω 0 and v on geometrical parameters of the system. In particular, two dependencies are studied in details: • Dependence on the period a when a is slightly different than a 0 : this dependence is responsible to the slight difference between U (1) , ω (1) 0 and U (2) , ω (2) 0 corresponding to upper and lower gratings of period a 1 and a 2 . The results of this study are shown in Fig. S9c. We notice that the linear dependence ω 0 (a) leading to a simple proportional relation between (2) 2 and 1 N ≈ a2−a1 a0 . As a consequence, the three parameters ∆,∆ U and q 0 of the Hamiltonian (S38) are connected and can be reduced to a single one, for example q 0 .
• Dependence on the filling fraction κ: the strong and almost linear dependence of U (κ) is shown in Fig. S9d. It suggests that the filling fraction is the parameter for tuning the intralayer coupling strength. (d) Dependence of the retrieved parameters when the filling fraction κ is scanned from 0.5 to 1. It shows that while the offset energy ω0 and the group velocity are slightly modified, the intra-layer coupling strength U is greatly modified from 3U0 to 0. .

Parameter retrieval of bilayer structure
The dispersion characteristic of bilayer structure can be analytically calculated from Eq. (S22) from the main text. The detailed of these eigenmodes has been reported in [19]. Here we only discuss how to retrieve the inter-layer coupling strength from these band structure and the validation of the method.
Since ω 0 and U are already retrieved from the simulation of single grating, only V left to be retrieved. One may show that, for AA stacking (i.e. δ/a = 0), the band structure consist of four bands with bandedge energies given by ω 0 ± U + V and ω 0 ± −V . As a consequence, V is directly extracted from the energy of resonance at q = 0 of the RCWA simulations for anyone from the four bands. Using this method, we can easily obtain the dependence of V as the function of the distance L separating the two grating. The results shown in Fig. S10b evidences the dependence law V = V 0 e −L/L0 used in the main text.
Finally, we confirm the validity of the retrieved parameters by using them to calculate the band structure of the bilayer for different relative shift δ/a, and for diffrent value pof L. The results presented in Fig. S10c show perfect aggreement between the calculated dispersion and the ones obtained by RCWA simulations, thus validate the retrieved parameters.

BAND EDGES OF MOIRÉ BANDS
To investigate the interplay between intra and inter-layer coupling in the formation of moiré bands, the band-edge energies (at Γ and X points) of electron-like and hole-like moiré bands are extracted from effective Hamiltonian calculations when scanning Figure S10. (a) Sketch of a bilayer grating structure. (b) Band-edge energies of the band structure of a bilayer grating as a function of the distance L between the two layers. The two grating are identical and aligned, with a = a0, κ = 0.8 and h = 0.6a0. The blues circles correspond to extracted data from RCWA simulation. The solid red lines are fittings, given by ω0 ± U − V and ω0 ± U + V . Here ω0 = Ω0 and U = U0, obtained from parameter retrieval of the single grating. And V (L) = V0e −L/L 0 with V0a0/2πc = 0.032 and L0/a0 = 0.34. (c) Band structure of bilayer grating structures of different relative displacement δ/a, obtained by RCWA simulation and by the effective Hamiltonian using retrieved parameter U, ω0, v and V . the ratio V /U for different moiré configurations with fixed value of U = U 0 . The results depicted in Fig S11 evidence two important features: • The magic configuration takes place at the crossings of band edge energies from the same miniband when tuning V /U (indicated by green arrows in Figs S11).
• The two moiré bands get closer when increasing V /U , as previously discussed when scanning L in the maintext. Interestingly, the gap between them is closed, and they merge together when V /U 2 for all value of N . Indeed, the bandgap when the two gratings are slightly different (i.e. N 1) and uncoupled (i.e. V U ) is given by the gap of a single grating, thus amounts to 2U . When the interlayer layer coupling V is implemented, the two moirés bands emerge and are separated to the corresponding continuum by a quantity ∼ V . Thus they would merge at the zero energy when V ∼ U . This feature is not revealed from the numerical simulation since the maximum value of V /U from our design is 1.76 (i.e. L = 0 and κ = 0.8).