Topological Gaps by Twisting

It is shown that twisted $n$-layers have an intrinsic degree of freedom living on $2n$-tori, which is the phason supplied by the relative slidings of the layers and that the twist generates pseudo magnetic fields. As a result, twisted $n$-layers host intrinsic higher dimensional topological phases and those characterized by second Chern numbers can be found in a twisted bi-layer. Indeed, our investigation of phononic lattices with interactions modulated by a second twisted lattice reveals Hofstadter-like spectral butterflies in terms of the twist angle, whose gaps carry the predicted topological invariants. Our work demonstrates how multi-layered systems are virtual laboratories for studying the physics of higher dimensional quantum Hall effect and how to generate topological edge chiral modes by simply sliding the layers relative to each other. In the context of classical metamaterials, both photonic and phononic, these findings open a path to engineering topological pumping via simple twisting and sliding.

Twisted graphene bilayers were recently observed to have exceptional spectral characteristics [51][52][53] that host extremely rich single-and many-body physics [54][55][56][57]. Exciting phenomena occurring at magic angles have been observed, such as superconductivity in the flat bands [54] of twisted graphene bilayers and, very recently, hyperbolic and elliptic dispersion of polaritions in twisted photonic systems [58]. Definitely, the field of twistronics, consisting of spectral and dynamical engineering by twisting layered materials and metamaterials, is one of the most active research fields at the present time. Interestingly, in [59] it was found that the spectral gaps stabilized by the magic angles carry non-trivial (fragile) topological indices. Another exciting topological finding is that, under a modest magnetic field, the bilayered graphene opens topological gaps that host correlated Chern insulating states [60].
Moiré patterns, as any other aperiodic pattern, can be analyzed, at least formally, with the operator algebraic program [61][62][63] pioneered by Bellissard and others, long before the bilayers came to the full attention of the physics community. The power of this formalism was already demonstrated for incommensurate multilayered systems in [64], where a solution of the extremely difficult problem of computing the transport coefficients was given. However, the topological part of this comprehensive program, as applied to Moiré patterns, was missing. Zak [65] and Avron et al [66] thought us to look for winding and Chern numbers if an adiabatic parameter is given and lives on circles or tori, but one should contemplate that there is no obvious circle or torus associated with twisted bilayers. The work of Bellissard [61] indicates that such smooth manifolds could come in the form of the hull of the pattern, a concept that he introduced and referenced here as the phason space. For the bilayered systems described below, we compute this hull (more precise, the transversal of the pattern) and find it to be the 2-torus. Following further Bellissard's program, we prove the following statement: if one arranges and couples identical resonator in a periodic lattice and modulates the couplings by a second twisted lattice, the dynamical matrices always land in the non-commutative 4-torus, regardless of the details. This statement completely classifies the dynamics for this class of systems. As an automatic followup, we can announce that the bilayered systems support topological phases characterized by the 2 nd -Chern number, without any fine tunning or external magnetic fields. Furthermore, we show that topological edge chiral modes can be generated by simply sliding the layers relative to each other, hence supplying a robust and effective mechanism for topological pumping. Our analysis can be straightforwardly generalized to multi-layered systems such as tri-layers, where even higher virtual dimensional topological phases appear and the relative slidings of the layers can be achieved in more interesting ways.
To demonstrate the above principles, we consider three generic twisted phononic systems of different lat- Top view displaying lattice vectors a 1 and a 2 separated by an angle β. In this un-twisted configuration, the lattice sites are aligned with the peaks of the potential. (c) Also in a top view, the potential is twisted by an angle θ relative to the lattice. (d,e) Examples of truly aperiodic twisted bilayers for θ = π/50 and θ = π/25, respectively. The figures were generated for a 2 = 4 √ 2a 1 and β = π/ √ 7.
tice symmetries, for which we map the resonant spectra as a function of twist angle. Without any tuning, the computations reveal Hofstadter-like spectral butterflies and, as we shall see, the integrated density of states (IDS) evaluated inside the spectral gaps gives rise to well defined but intricate patterns of curves. The qualitative shape of these curves change from one example to another, yet we put forward one unifying prediction which fits every single pattern seen in our IDS plots. This prediction follows from the K-theory of the non-commutative 4-torus and the agreement with the numerical computations leave no doubt that the dynamical matrices for these systems fall into that algebra, which is the main statement of our work.
There are important practical consequences of our findings. The phason can be moved on its rightful space by sliding the bilayers relative to each other. The twist angle dictates how many chiral modes are generated by a cycle. This is by far the simplest and easiest way to manipulate the phason of an aperiodic pattern and produce topological pumping (but see also [40]). In the presence of an edge, looping the phason around the fundamental loops of this torus results in topological chiral edge modes. The reader will find in this paper a bulk-boundary principle, tested and confirmed by the numerical observations, which predicts the number of these chiral bands from the values of the bulk topological invariants. By that, we demonstrated that twistronics supplies new ways to generate and manipulate topological edge excitations with unprecedented control and precision. For bilayer graphene, for example, a simple vibration of the layers relative to each other should reveal the existence of the predicted topological edge modes.
On the computational side, let us recall that twisted bilayers away from the special angles are notoriously difficult to deal with because of lack of periodic approximants [64]. For the present context, things are made more difficult by the topological edge states which contaminate the bulk spectral gaps. This is a general prob-lem for aperiodic topological systems and is also encountered, for example, in topological quasicrystals [20,68]. We found that the periodic boundary conditions eliminate these topological edge states but impurity-like edge states still persist. The latter, however, have a low density and, as a consequence, the maps of the density of states, as opposed to the spectrum itself, supply remarkably clean pictures of the spectral butterflies [69]. The quality of the numerical simulations is reflected in sharpness of the integrated density of states inside the spectral gaps, which is the essential numerical outcome that is used to compare with the theoretical predictions and extract the topological invariants.

A. The mechanical system defined
We consider a lattice L 1 of identical masses of fixed (x, y) coordinates r n = n 1 a 1 + n 2 a 2 , n = (n 1 , n 2 ) ∈ Z 2 , where a 1 and a 2 are arbitrary lattice vectors separated by an angle β (Fig. 1). We denote their magnitudes by a 1 and a 2 , respectively. The masses move along the z direction and interact via two-body potentials, while each mass experiences an external potential, represented by the colored surface in Fig. 1. Generically, such system is described by a Lagrangian The functional form of V θ is whereR θ is the rotation matrix by θ of the (x, y)-plane and V 0 is a periodic potential V 0 (r + r n , z) = V 0 (r, z) for all r n ∈ L 1 . We define L 2 =R −1 θ [L 1 ] to be the twisted lattice, such that V θ (r + r n ) = V θ (r) for all r n ∈ L 2 .
In Figs. 1(d,e), we illustrate two configurations of the system, corresponding to a generic lattice L 1 and two twist angles θ such that the periodicity of the system cannot be restored no matter what super-cell is used. The latter can only happen for a discrete set of θ-s, hence, the generic cases are those of purely aperiodic configurations. Our analysis will cover both the special and generic cases on equal footing. Let us recall that a fine sampling of the twists by commensurate angles requires notoriously large supercells [70] and we want to assure the reader, and especially the experimentalists, that our analysis and conclusions do not rely on any such periodic approximants. Let us mention that z can be replaced with any other local degree of freedom, such as a rotation angle. In that case, laboratory models of the system introduced above can be implemented, for example, with the systems of magnetically coupled spinners introduced in [43,71]. This task, however, is left to the future for now.
In the regime of small oscillations, the dispersion equation of the collective resonant modes takes the form − W r n − r n ,z n −z n ζ n .
Here, ζ n = z n −z n withz n being the equilibrium zcoordinates of the masses and Note that, in general, both the potential and the coupling constants are perturbed by the L 2 lattice. This will be considered in our theoretical analysis but left aside in our numerical experiments.

B. Numerical simulations
In our numerical applications, we considered a shortrange two-body interaction such that Let us be clear that the pair interactions were not truncated to first nearest neighbors and, instead, all pairs of masses interact even though the interaction might be exponentially small. We chose a potential such that where b i -s are L 1 's reciprocal vectors. The resulting dynamical matrix for the system of equations (3) was exactly diagonalized on a 100 × 100 resonator lattice with periodic boundary conditions (PBC), while sampling θ over 1000 equally spaced points in the interval [0, π]. The mass m was set to 1. We chose three representative lattices such that, at one end, we have a square lattice with a large point group symmetry and, on the other end, a lattice where both β/2π and a 1 /a 2 are irrational numbers such that the point group is trivial. The main reason for this is to convince the reader that our statements are independent of the point symmetry of the lattice, as it should for topological phases from class A. A deeper reason for these choices will be revealed in section II D. Figs. 2(a), 3(a) and 4(c) report the resonant spectra of these systems as functions of θ. Large bulk spectral gaps contaminated by edge spectrum are visible in all cases and, overall, the spectra project the same kind of fractality seen in the Hofstadter butterfly [72]. Let us specify that PBC prevents the topological edge modes but impurity edge states still persist because periodicity is broken by the twisted potential. Nevertheless, PBC are preferred because, while the impurity states still contaminate the bulk gaps, they do not display any spectral flow with θ, as it is evident in all our results. An important numerical finding is that the spectra can be almost entirely cleared of the edge states contamination by computing the corresponding density of states. This is exemplified in Figs. 2(b), 3(b) and 4(b). The symmetries of the spectral butterflies are now more apparent. For example, for the square lattice we have reflection symmetries about mid horizontal and vertical lines, as well as θ → π/2 − θ. This is why we only labeled gaps in the left half of the spectral butterfly. For the lattice with a 1 = a 2 and β = π/ √ 7, the reflection symmetry w.r.t. the mid vertical line is still present, and we still focus on spectral gaps from the left side of the spectral butterfly. For the most generic lattice a 1 a 2 and β = π/ √ 7, all the symmetries are lifted and we will investigate spectral gaps from all parts of the spectral butterfly.
The most interesting outcomes of our simulations are the integrated density states (IDS), defined as where Spec(D) is the spectrum of the dynamical matrix, i.e. the set of eigenvalues counted with their degeneracies. Throughout, | · | represents the cardinal of a set. The maps of the IDS as function of θ and Ω 2 are reported in Figs. 2(c), 3(c) and 4(c) for the three lattices considered in our study. Since the IDS(E) is constant when E takes values in the spectral gaps, the 3-dimensional IDS plots have an abrupt variation with respect to E whenever E traverses a gap. In the color maps shown in our figures, these variations appear as abrupt changes of the color and these features were further enhanced by using appropriate lightning. As a result, the values of the IDS inside the gaps can be easily identified by the dark lines visible in all our plots. As one can see, there are drastic changes in the pattern of these lines from one lattice to another. Yet, as we shall see, all IDS curves in these three figures are described by one unifying equation, where the topological invariants appear as integer coefficients.

C. Theoretical Interpretation
In this section we explain the features seen in the numerical experiments using the K-theoretic tools developed in [61][62][63]. These works demonstrated the existence of a standard formalism, but the statements are not entirely constructive, hence, every new application poses certain computational challenges. As such, it is remarkable that Bellissard's program can be carried out entirely and explicitly for the present context.
Algebra of dynamical matrices. We encode the degrees of freedom in the vector |Z = n∈Z 2 ζ n |n and transform the dispersion equations into ω 2 |Z = D|Z with the dynamical matrix D = m,n w m,n (P) |m n| (8) written here in the most generic form. As we shall see, the details of the coupling coefficients are not important for this analysis. What is important is that they are fully determined by the pattern P, which consists of the union of the resonator lattice L 1 and potential lattice L 2 . In other words, if the potential and the type of resonators are fixed and no external intervention is allowed, there are pre-defined functions w m,n (P) of variable P that supply the couplings. Upgrading the coupling coefficients to coupling functions is a strategic point in the theory of dynamics over patterns. For our explicit model, these functions are already specified in (3). In typical metamaterial experiments, these functions can be mapped entirely by exploring the pattern space as it was done, for example, in [43]. Once these functions are cataloged, one can evaluate them on a particular pattern and generate the dynamical matrix. The next important observation is Galilean invariance, which says that if we rigidly translate P, hence both lattices, the coupling functions must display the following covariance relations [73]: w m−a,n−a (τ a P) = w m,n (P) or w m,n (P) = w m−n,0 (τ n P), (9)  where τ a P is the rigid shift of the pattern which brings the resonator labeled by a ∈ Z 2 to the origin. τ a P is also the pattern seen by an observer which jumped from the origin to the site a of the resonator lattice L 1 . After these observations are in place, something magic happens [73]. Indeed, we can drop one redundant index and, using q = m − n as well as the shift operator S q |n = |n + q , D takes a very particular form The extraordinary conclusion is that any Galilean invariant dynamical matrix over P is generated from a small algebra generated by the elementary shift operators S i , i = 1, 2, and by diagonal operators T f = n f (τ n P) |n n| with f a function on the space of patterns. Furthermore, one can check the commutation relation (CR) n f (τ n P)|n n|S q = S q n f (τ n+q P) |n n|, which can be written more compactly as In the following, we compute this algebra explicitly and show that it is isomorphic to the non-commutative 4torus. For this, we need first to parameterize the space of patterns. As in Fig. 5, it is useful to imagine an observer sitting on top of a resonator. Looking around, one sees a certain pattern P and if the observer jumps to another resonator, say a hundred lattice units away, one will perceive a completely different pattern. The question is then, what is the minimum information the observer needs to reproduce the entire pattern, if we place the observer on top of an arbitrary resonator. Of course, the observer knows that one is dealing with a bilayer and that L 2 is twisted by θ relative to L 1 . The answer is quite simple. The observer projects hers/his location onto the L 2 plane and this projection ξ necessarily falls in one primitive cell of L 2 . The observer sees the same pattern if this point lands on the opposite sides of the primitive cell, hence this primitive cell should be wrapped as a torus. In fact, the best strategy is to think that the entire second plane has been folded over one single primitive cell, e.g. the one shaded in Fig. 5. In other words, we work with the torus R 2 /L 2 , which is parametrized as Ξ = (R mod a 1 ) × (R mod a 2 ), hence its points are ξ = (ξ 1 , ξ 2 ), ξ i ∈ R mod a i . Now, the only information the observer needs in order to re-draw P is the position of its projection ξ in this primitive cell, hence ξ is the phason of the aperiodic pattern. Indeed, suppose that both layers have been erased. In this case, the observer will use a i to redraw L 1 . Then, using the coordinates (ξ 1 , ξ 2 ) together with the given twist angle, the observer re-traces the primitive cell of L 2 immediately above her/him and then tiles the plane by periodic translations of this primitive cell. An important question is if the observer explores the whole Ξ or only a part of it, as she/he jumps from one resonator to another. Of course, this question is equivalent to asking if the dynamical system τ a P described above is topologically ergodic.

Potential Layer
Resonator Layer An observer can reproduce the entire bilayer from the coordinates (ξ 1 , ξ 2 ) of the origin of L 1 relative to the shaded primitive cell of L 2 . If the observer moves to different resonators on the L 1 lattice, this results in shifts over the folded R 2 /L 2 space generated by τ 1 (ξ) and τ 2 (ξ). From Fig. 5, we can see that these are just translations of the torus. In terms of coordinates ξ i , the generators of these translations are given by with: As it is well known, if at least two A i j 's are irrational numbers, then the orbit of one single point under repeated translations (13) fills the torus densely, hence the dynamical system (Ξ, τ) is topologically ergodic. One should not be surprised by the existence of this dynamical system because (Ξ, τ) is just the hull of P, predicted to exist for any point pattern in [61][62][63].
We now have all the information to compute the algebra which generates the dynamical matrices. Since any continuous function over torus accepts a discrete Fourier decomposition, written slightly differently below, the algebra of T f operators has two generators T j corresponding to the elementary functions: f → u j (ξ) = e ı2πξ j /a j , j = 1, 2.
Furthermore, since (u j • τ i )(ξ) = e ı2π A i j u j (ξ), the CR's (12) become S i T j = e −ı2π A i j T j S i . As such, the algebra which contains all Galilean invariant dynamical matrices is generated by four unitary elements: Let us point out that, by considering two degrees of freedom per resonator, additional entries of the Φ-matrix can be populated by ± 1 2 values, as explained in [23] for a spinful model.
Predictions via K-Theory. In K-Theory [75], the set of projections that can be (stably) deformed into a given projection is called the K 0 -class of that projection. It is the complete topological invariant associated to that projection, in the sense any other topological invariant is already determined by its K 0 -class. These classes of projections can be added and subtracted, hence they form an abelian group, the K 0 -group of the algebra. Two homotopic projection P and P are also similar: P = U * PU for some unitary element from the same algebra. If T is a trace on the algebra, then automatically T (P) = T (P ) because we are allowed to make cyclic permutations inside a trace. This means that any trace is constant over the homotopy classes of projections. As such, traces are bona-fide topological invariants.
In the case of non-commutative 4-torus, the K 0 -group has 2 d−1 -generators {e J }, conveniently labeled by a sub- set J ⊆ {1, 2, 3, 4} of directions with |J| = even [74]. Furthermore, for generic Φ-matrices, the non-commutative 4-torus accepts a unique trace T , which coincides with the trace per volume [61]. Any gap projection P G of a dynamical matrix defines a K 0 -class and accepts a decomposition in terms of the generators [P G ] 0 = J n J [e J ] 0 . The integer numbers n J are called gap labels [62] and, in general, they represent the complete set of independent topological invariants that can be associated to a gap projection. They are related but not necessarily equal to the Chern numbers (see below). Since traces are linear maps, On the other hand, hence T (P G ) = IDS(G). As such, if we can resolve the values of the trace on the generators of the K 0 -group, we can make a prediction about the allowed values of IDS. For the non-commutative torus, this extremely useful piece of information was supplied in [76], and we have where Φ J is the matrix Φ restricted to indices J and Pf is the pfaffian of the resulting anti-symmetric matrix. In our case, this gives the prediction [77] When there are no linear relations with integer coefficients between A i j -s, we can compute all topological invariants supplied by n {i, j} -s by fitting (21) to the numerically obtained IDS curves in Figs. 2(c), 3(c) and 4(c). Unfortunately, Det(A) = 1, hence we can only determine the sum n ∅ + n {1,2,3,4} via this procedure.
The values of the Chern numbers on the K 0 -generators were computed in [67][p. 141]: for •-gap from Fig. 4(b).
However, due to the symmetry under π/2 rotations, n {1,3} = n {2,4} and n {1,4} = −n {2,3} , which follows directly from the expressions of the first Chern numbers. As such, (24) can be used to determine all topological invariants supplied by n {i,j} . We found that Eq. (24) fits perfectly all IDS curves seen in Fig 2(c). Fittings of the IDS curves inside the eight large gaps identified in Fig. 2(b) are reported in Fig. 6(a), together with the topological invariants extracted from the fittings. Let us remark that, by using the symmetries of the spectral butterfly, we can automatically fit many more IDS curves, 56 to be more precise.
For the bilayer analyzed in Fig. 3, one additional linearly independent term is present in the IDS expression: Form symmetry considerations, we found again that n {1,4} = −n {2,3} . Again, we have verified that Eq. (25) perfectly fits all IDS curves seen in Fig 3(c). Fittings of the IDS curves inside the eight large gaps identified in Fig. 3(b) are reported in Fig. 6(b), together with the topological invariants extracted from the fitting. The symmetry of the spectral butterfly can be used to automatically fit eight additional IDS curves in the right side of the IDS plot. Finally, for the bilayer analyzed in Fig. 4, we have five linearly independent terms present in the IDS expression: In this case there is no point symmetry left and the four topological numbers n {i, j} are all independent. We found again that Eq. (25) fits perfectly all IDS curves seen in Fig 4(c). Fittings of the IDS curves inside the eight large gaps identified in Fig. 4(b) are reported in Fig. 6(c), together with the topological invariants extracted from the fitting.
Let us point out that every single gap among the 24 gaps analyzed in Fig. 6 displays a non-zero n {i, j} but we have not yet able to demonstrate the existence on nontrivial top invariants n {1,2,3,4} . For this, we turn to the bulk-boundary correspondence for the twisted bilayers.

E. Bulk-boundary correspondence and existence of 2 nd -Chern states
Physical boundaries are created by restricting either one of n k coefficients of r n to non-negative values. If n k ≥ 0, then the boundary cuts the k-th direction and we will call it a k-boundary. We denote the resulting dynamical matrix by D k (ξ). The bulk-boundary for class A in higher dimensions states [67] that the surface states admit topological invariants in the form of odd Chern numbers and that there is a precise relation between all bulk and surface invariants. In particular, for our lower Chern numbers [67][p. 175], (27) where N ki is the net number [78] of eigenvalues of D k (ξ 1 , ξ 2 ) that cross an arbitrary reference line inside the bulk gap when ξ i is varied from 0 to a i while holding the other ξ fixed. Using (22), we can write the bulkboundary principle explicitly, Numerically, we generate a k-boundary by using open boundary conditions in that physical direction and periodic boundary condition in the remaining direction. As always, we will create a pair of boundaries, hence the numerically computed edge modes will always come in pairs. In Fig. 7(a-d), which was simulated on a 21 × 21 lattice, we analyze the bulk-boundary correspondence for the •-gap from Fig. 4(b). As one can see, D 1 (ξ) displays 21 positively sloped chiral bands when ξ 2 is varied [79], and no chiral bands are present for the other three cases. This is consistent with n {1,4} = 1 and trivial values for the other invariants. Similarly, in Fig. 7(e-f), which was simulated on a 26 × 26 lattice, we analyze the bulkboundary correspondence for the •-gap from Fig. 4(b).
In this case, D 1 (ξ) displays 27 positively sloped chiral bands when ξ 1 is varied and no chiral bands are present for the other three cases. This is consistent with n {1,3} = 1 and trivial values for the other invariants. These numerical findings confirm the predicted bulk-boundary correspondences based on (28) and the data from Fig. 6 and, furthermore, they enable us to actually conclude that n {1,2,3,4} = 0 for these two particular gaps.
Additional boundary spectra are reported in Fig. 8, which were simulated on a 23 × 23 lattice and θ = 0.2.
They correspond to the •-gap in Fig. 4(b). From the data in Figs. 8 and 6, and by assuming n {1,2,3,4} = −1, we have: As one can see, the predicted bulk-boundary correspondence (28) holds with a remarkable precision given the relatively small size of the lattice [80]. Similar agreements hold true for other spectral gaps from Fig. 4(b) and, for example, for the •-gap we found n {1,2,3,4} = 1.
As such, twistronics is capable of generating gaps with 2 nd -Chern numbers.

III. DISCUSSION
We have demonstrated that twistronics can be a simple yet extremely effective way to produce topological gaps and topological boundary modes. Indeed, twisted bilayers have a "hidden" degree of freedom, the phason ξ, which leaves on a torus and can be controlled by simple relative shifts of the layers. For generic twist angles, these shifts do not affect the bulk spectrum, hence the bulk gaps, but they generate dispersive chiral boundary modes in the presence of a boundary. The count of these modes agrees with a precise topological bulk-boundary principle.
To our knowledge, this is the first time when the algebra of dynamical matrices for a Moiré pattern has been explicitly computed. With that result at hand, the Ktheoretic machinary invented by Bellissard [61] enabled us to classify all topological phases from class A supported by these twisted bilayers and to produce a highthroughput of topological gap labels. Let us mention that the topological invariants can be computed directly using the algorithms developed in [81]. However, those algorithms require a substantial computational effort and, as such, the technique based on the IDS fitting is an important outcome of our work. Using the newer results from [67], we were able to also make precise predictions about the bulk-boundary correspondence for these Moiré patterns. To our knowledge, it is the first time when a bulk-boundary correspondence is observed for non-integer invariants.
Our analysis generalizes to the cases where there are more degrees of freedom per primitive cell, such as the honeycomb lattice, or when the coupling constants are modulated not by one but by multiple twisted lattices. Indeed, assume that a lattice L 3 is added on top of L 2 in Fig. 5. By following similar arguments, it is easy to see that reproducing the whole patterns requires the knowledge of the projection of the resonator where the observer sits on both R 2 /L 2 and R 2 /L 3 tori. As such, the phason space is a 4-torus and the dynamical system τ can be computed by similar methods. It follows that the algebra which generates the dynamical matrices is the non-commutative 6-torus, which hosts topological phases with 3 rd -Chern numbers [67].
In conclusion, we proved that the twisted layered systems, both classical and quantum, can be resourceful virtual laboratories for exploring completely new physics and topological states. For metamaterials research, our findings open new venues for engineering topological gaps and robust boundary modes without any need for fine-tuning or active components. Data Availability. The data and the numerical codes are available from the authors upon request. [2] J. He, Y. Liang, S.-P. Kou, Topological hierarchy insula-