Topological gaps by twisting

Abstract

Twisted bilayered systems such as bilayered graphene exhibit remarkable properties such as superconductivity at magic angles and topological insulating phases. For generic twist angles, the bilayers are truly quasiperiodic, a fact that is often overlooked and that has consequences which are largely unexplored. Herein, we uncover that twisted n-layers host intrinsic higher dimensional topological phases, and that those characterized by second Chern numbers can be found in twisted bi-layers. We employ phononic lattices with interactions modulated by a second twisted lattice and reveal Hofstadter-like spectral butterflies in terms of the twist angle, which acts as a pseudo magnetic field. The phason provided by the sliding of the layers lives on 2n-tori and can be used to access and manipulate the edge states. Our work demonstrates how multi-layered systems are virtual laboratories for studying the physics of higher dimensional quantum Hall effect, and can be employed to engineer topological pumps via simple twisting and sliding.

Introduction

Engineering topological states using aperiodic principles is an extremely active area of research, spread over different fields such as condensed matter^{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20}, photonics^{21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37}, acoustics^38,39,40, and mechanics^{41,42,43,44,45,46,47,48,49,50}. Common to all these is the existence of an intrinsic degree of freedom, the phason, which in many instances is experimentally accessible and fully controllable. The phason space augments the physical space and supplies additional virtual dimensions³, hence enabling physical phenomena beyond what can be ordinarily observed in our physical space. In particular, the phason can be used as an adiabatic knob to engineer topological pumping, which is one of the main applications of the formalism, as evidence by the literature cited above.

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⁵⁴ of twisted graphene bilayers, and very recently, hyperbolic and elliptic dispersion of polaritions in twisted photonic systems⁵⁸. 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⁵⁹ it was found that the spectral gaps stabilized by the magic angles carry nontrivial (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⁶⁰.

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⁶⁴, 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⁶⁵ and Avron et al.⁶⁶ 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⁶¹ 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 precisely, 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 noncommutative 4-torus, regardless of the details. This statement completely classifies the dynamics over this class of systems. As an automatic follow-up, we can announce that the bilayered systems support topological phases characterized by the 2nd-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 multilayered 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 lattice 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 noncommutative 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.

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⁶⁴. For the present context, things are made more difficult by the topological edge states which contaminate the bulk spectral gaps. This is a general problem for aperiodic topological systems and is also encountered, for example, in topological quasicrystals^20,67. We found that the PBC 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 (we learned of this effective solution from²⁰). The quality of the numerical simulations is reflected in sharpness of the IDS inside the spectral gaps, which is the essential numerical outcome that is used to compare with the theoretical predictions and extract the topological invariants.

There are several 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 ref. ⁴⁰). 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. Therefore, our findings open new paths for exploring higher dimensional Quantum Hall effects and manipulating edge states in multilayered systems via simple twisting and sliding.

Results

The mechanical system defined

We consider a lattice ${{\mathcal{L}}}_{1}$ of identical masses of fixed (x, y) coordinates r_n = n₁ a₁ + n₂ a₂, ${\boldsymbol{n}}\,=\,({n}_{1},{n}_{2})\,\in\, {{\mathbb{Z}}}^{2}$, where a₁ and a₂ are arbitrary lattice vectors separated by an angle β (Fig. 1a–c). We denote their magnitudes by a₁ and a₂, 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. 1a. Generically, such system is described by a Lagrangian

$${\mathcal{L}} = \, \mathop{\sum}\limits_{{\boldsymbol{n}}}\left(\frac{1}{2}m\ {\dot{z}}_{{\boldsymbol{n}}}^{2}\,-\,{V}_{\theta }({{\boldsymbol{r}}}_{{\boldsymbol{n}}},{z}_{{\boldsymbol{n}}})\right)\\ \, -\frac{1}{2}\mathop{\sum}\limits_{{\boldsymbol{n}},{\boldsymbol{n}}^{\prime} }W({{\boldsymbol{r}}}_{{\boldsymbol{n}}}\,-\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}^{\prime} },{z}_{{\boldsymbol{n}}}\,-\,{z}_{{\boldsymbol{n}}^{\prime} }).$$

(1)

The functional form of V_θ is

$${V}_{\theta }({\boldsymbol{r}},z)\,=\,{V}_{0}\left({\hat{R}}_{\theta }[{\boldsymbol{r}}],z\right),$$

(2)

where ${\hat{R}}_{\theta }$ is the rotation matrix by θ of the (x, y)-plane and V₀ is a periodic potential V₀(r + r_n, z) = V₀(r, z) for all ${{\boldsymbol{r}}}_{n}\,\in\, {{\mathcal{L}}}_{1}$. We define ${{\mathcal{L}}}_{2}\,=\,{\hat{R}}_{\theta }^{-1}[{{\mathcal{L}}}_{1}]$ to be the twisted lattice, such that ${V}_{\theta }({\boldsymbol{r}}\,+\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}}^{\prime})\,=\,{V}_{\theta }({\boldsymbol{r}})$ for all ${{\boldsymbol{r}}}_{{\boldsymbol{n}}}^{\prime}\,\in\, {{\mathcal{L}}}_{2}$.

Fig. 1: Mechanical systems defined.

a Mechanical lattice (top layer) with an underlying potential surface (bottom layer). b Top view displaying lattice vectors a₁ and a₂ 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}\,=\,\root {[4]}\of {2}{a}_{1}$ and $\beta \,=\,\pi /\sqrt{7}$.

In Fig. 1d, e, we illustrate two configurations of the system, corresponding to a generic lattice ${{\mathcal{L}}}_{1}$ and two twist angles θ such that the periodicity of the system cannot be restored no matter what supercell 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⁶⁸ 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,69. 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

$$m{{{\Omega }}}^{2}{\zeta }_{{\boldsymbol{n}}} = \, \left({K}_{{\boldsymbol{n}}}\,+\,V^{\prime\prime}_{\theta} \left({{\boldsymbol{r}}}_{{\boldsymbol{n}}},{\bar{z}}_{{\boldsymbol{n}}}\right)\right){\xi }_{{\boldsymbol{n}}}\\ \, -W^{\prime\prime} \left({{\boldsymbol{r}}}_{{\boldsymbol{n}}}\,-\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}^{\prime} },{\bar{z}}_{{\boldsymbol{n}}}\,-\,{\bar{z}}_{{\boldsymbol{n}}^{\prime} }\right)\ {\zeta }_{{\boldsymbol{n}}^{\prime} }.$$

(3)

Here, ${\zeta }_{{\boldsymbol{n}}}\,=\,{z}_{{\boldsymbol{n}}}\,-\,{\bar{z}}_{{\boldsymbol{n}}}$ with ${\bar{z}}_{{\boldsymbol{n}}}$ being the equilibrium z-coordinates of the masses and

$${K}_{{\boldsymbol{n}}}\,=\,\mathop{\sum}\limits_{{\boldsymbol{n}}^{\prime} }W^{\prime\prime} \left({{\boldsymbol{r}}}_{{\boldsymbol{n}}}\,-\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}^{\prime} },{\bar{z}}_{{\boldsymbol{n}}}\,-\,{\bar{z}}_{{\boldsymbol{n}}^{\prime} }\right).$$

(4)

Note that, in general, both the potential and the coupling constants are perturbed by the ${{\mathcal{L}}}_{2}$ lattice. This will be considered in our theoretical analysis but left aside in our numerical experiments.

Numerical simulations

In our numerical applications, we considered a short-range two-body interaction such that

$$W^{\prime\prime} ({{\boldsymbol{r}}}_{{\boldsymbol{n}}}\,-\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}^{\prime} })\,=\,{e}^{-3| {{\boldsymbol{r}}}_{{\boldsymbol{n}}}\,-\,{{\boldsymbol{r}}}_{{\boldsymbol{n}}^{\prime} }{| }^{2}}.$$

(5)

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

$${V}_{\theta }^{^{\prime\prime} }({{\boldsymbol{r}}}_{{\boldsymbol{n}}})\,=\,0.1\ \left(\cos \left({{\boldsymbol{b}}}_{1}\cdot {\hat{R}}_{\theta }[{{\boldsymbol{r}}}_{{\boldsymbol{n}}}]\right)\,+\,\cos \left({{\boldsymbol{b}}}_{2}\cdot {\hat{R}}_{\theta }[{{\boldsymbol{r}}}_{{\boldsymbol{n}}}]\right)\right),$$

(6)

where b_i-s are ${{\mathcal{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₁/a₂ 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 2.4. Figs. 2a–c, 3a–c, 4a–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⁷⁰. 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.

Fig. 2: Spectral properties of square lattice with β = π/2, a₁ = a₂ = 1.

a Eigenfrequencies (Ω²) of 100 × 100 lattice calculated with periodic boundary conditions as a function of θ. The same result is repeated in (b), but color-coded by the associated density of states (colorbar). Particular gaps are marked by colored dots and labeled (I–VIII) for reference. c Integrated density of states (IDS) as a function of θ, where the color gradient represents Ω².

Fig. 3: Spectral properties of lattice with $\beta \,=\,\pi /\sqrt{{\bf{7}}},{{\bf{a}}}_{{\bf{1}}}\,=\,{{\bf{a}}}_{{\bf{2}}}\,=\,{\bf{1}}$.

a Eigenfrequencies (Ω²) of 100 × 100 lattice calculated with periodic boundary conditions as a function of θ. The same result is repeated in (b), but color-coded by the associated density of states (colorbar). Particular gaps are marked by colored dots and labeled (I–VIII) for reference. c Integrated density of states (IDS) as a function of θ, where the color gradient represents Ω².

Fig. 4: Spectral properties of lattice with $\beta \,=\,\pi /\sqrt{{\bf{7}}},{{\bf{a}}}_{{\bf{1}}}\,=\,{\bf{1}}/\sqrt[4]{{\bf{2}}},{{\bf{a}}}_{{\bf{2}}}\,=\,{\bf{1}}$.

a Eigenfrequencies (Ω²) of 100 × 100 lattice calculated with periodic boundary conditions as a function of θ. The same result is repeated in (b), but color-coded by the associated density of states (colorbar). Particular gaps are marked by colored dots and labeled (I–VIII) for reference. c Integrated density of states (IDS) as a function of θ, where the color gradient represents Ω².

An interesting 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. 2b, 3b, 4b. 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 (with colored dots and numbers (I–VIII)) in the left half of the spectral butterfly. For the lattice with a₁ = a₂ and $\beta \,=\,\pi /\sqrt{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₁ ≠ a₂ and $\beta \,=\,\pi /\sqrt{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 IDS, defined as

$${\rm{IDS}}({{\Omega }})\,=\,{\left.\frac{\left|{\rm{Spec}}(D)\cap (-\infty ,{{{\Omega }}}^{2}]\right|}{| {{\mathcal{L}}}_{1}| }\right|}_{{{\mathcal{L}}}_{1}\,\to\, {{\mathbb{Z}}}^{2}},$$

(7)

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 Ω² are reported in Figs. 2c, 3c, 4c 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.

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 $\left|Z\right\rangle \,=\,{\sum }_{{\boldsymbol{n}}\,\in\, {{\mathbb{Z}}}^{2}}{\zeta }_{{\boldsymbol{n}}}\ \left|{\boldsymbol{n}}\right\rangle$ and transform the dispersion equations into ${\Omega }^{2}\left|Z\right\rangle \,=\,D\left|Z\right\rangle$ with the dynamical matrix

$$D\,=\,\mathop{\sum}\limits_{{\boldsymbol{m}},{\boldsymbol{n}}}{w}_{{\boldsymbol{m}},{\boldsymbol{n}}}({\mathcal{P}})\ \left|{\boldsymbol{m}}\right\rangle \left\langle {\boldsymbol{n}}\right|$$

(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 ${\mathcal{P}}$, which consists of the union of the resonator lattice ${{\mathcal{L}}}_{1}$ and potential lattice ${{\mathcal{L}}}_{2}$. In other words, if the potential and the type of resonators are fixed and no external intervention is allowed, there are predefined functions ${w}_{{\boldsymbol{m}},{\boldsymbol{n}}}({\mathcal{P}})$ of variable ${\mathcal{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 meta-material experiments, these functions can be mapped entirely by exploring the pattern space as it was done, for example, in⁴³. Once these functions are cataloged, one can evaluate them on a particular pattern and generate the dynamical matrix.

The next relevant observation is Galilean invariance, which says that if we rigidly translate ${\mathcal{P}}$, hence both lattices, the coupling functions must display the following covariance relations⁷¹:

$${w}_{{\boldsymbol{m}}\,-\,{\boldsymbol{a}},{\boldsymbol{n}}\,-\,{\boldsymbol{a}}}({\tau }_{{\boldsymbol{a}}}{\mathcal{P}})\,=\,{w}_{{\boldsymbol{m}},{\boldsymbol{n}}}({\mathcal{P}})\,\,{\rm{or}}\,\,{w}_{{\boldsymbol{m}},{\boldsymbol{n}}}({\mathcal{P}})\,=\,{w}_{{\boldsymbol{m}}\,-\,{\boldsymbol{n}},{\boldsymbol{0}}}({\tau }_{{\boldsymbol{n}}}{\mathcal{P}}),$$

(9)

where ${\tau }_{{\boldsymbol{a}}}{\mathcal{P}}$ is the rigid shift of the pattern which brings the resonator labeled by ${\boldsymbol{a}}\,\in\, {{\mathbb{Z}}}^{2}$ to the origin. ${\tau }_{{\boldsymbol{a}}}{\mathcal{P}}$ is also the pattern seen by an observer which jumped from the origin to the site a of the resonator lattice ${{\mathcal{L}}}_{1}$.

After these observations are in place, something magic happens⁷¹. Indeed, we can drop one redundant index and, using q = m − n as well as the shift operator ${S}_{{\boldsymbol{q}}}\left|{\boldsymbol{n}}\right\rangle \,=\,\left|{\boldsymbol{n}}\,+\,{\boldsymbol{q}}\right\rangle$, D takes a very particular form

$$D\,=\,\mathop{\sum}\limits_{{\boldsymbol{q}}}{S}_{{\boldsymbol{q}}}\mathop{\sum}\limits_{{\boldsymbol{n}}}{w}_{{\boldsymbol{q}}}({\tau }_{{\boldsymbol{n}}}{\mathcal{P}})\ \left|{\boldsymbol{n}}\right\rangle \left\langle {\boldsymbol{n}}\right|.$$

(10)

The extraordinary conclusion is that any Galilean invariant dynamical matrix over ${\mathcal{P}}$ is generated from a small algebra generated by the elementary shift operators S_i, i = 1, 2, and by diagonal operators ${T}_{f}\,=\,{\sum }_{{\boldsymbol{n}}}f({\tau }_{{\boldsymbol{n}}}{\mathcal{P}})\ \left|{\boldsymbol{n}}\right\rangle \left\langle {\boldsymbol{n}}\right|$ with f a function on the space of patterns. Furthermore, one can check the commutation relation (CR)

$$\mathop{\sum}\limits_{{\boldsymbol{n}}}f({\tau }_{{\boldsymbol{n}}}{\mathcal{P}})\left|{\boldsymbol{n}}\right\rangle \left\langle {\boldsymbol{n}}\right|{S}_{{\boldsymbol{q}}}\,=\,{S}_{{\boldsymbol{q}}}\mathop{\sum}\limits_{{\boldsymbol{n}}}f({\tau }_{{\boldsymbol{n}}\,+\,{\boldsymbol{q}}}{\mathcal{P}})\ \left|{\boldsymbol{n}}\right\rangle \left\langle {\boldsymbol{n}}\right|,$$

(11)

which can be written more compactly as

$${T}_{f}\ {S}_{{\boldsymbol{q}}}\,=\,{S}_{{\boldsymbol{q}}}\ {T}_{f\circ {\tau }_{{\boldsymbol{q}}}}\,\,\forall \,{\boldsymbol{q}}\,\in\, {{\mathbb{Z}}}^{2}.$$

(12)

In the following, we compute this algebra explicitly and show that it is isomorphic to the noncommutative 4-torus. 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 ${\mathcal{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 ${{\mathcal{L}}}_{2}$ is twisted by θ relative to ${{\mathcal{L}}}_{1}$. The answer is quite simple. The observer projects hers/his location onto the ${{\mathcal{L}}}_{2}$ plane and this projection ξ necessarily falls in one primitive cell of ${{\mathcal{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 ${{\mathbb{R}}}^{2}/{{\mathcal{L}}}_{2}$, which is parametrized as ${{\Xi }}\,=\,({\mathbb{R}}\ {\rm{mod}}\,{a}_{1})\,\times\, ({\mathbb{R}}\ {\rm{mod}}\,{a}_{2})$, hence its points are ξ = (ξ₁, ξ₂), ${\xi }_{i}\,\in\, {\mathbb{R}}\ {\rm{mod}}\,{a}_{i}$.

Fig. 5: Derivation of the phason space.

An observer can reproduce the entire bilayer from the coordinates (ξ₁, ξ₂) of the origin of ${{\mathcal{L}}}_{1}$ relative to the shaded primitive cell of ${{\mathcal{L}}}_{2}$. If the observer moves to different resonators on the ${{\mathcal{L}}}_{1}$ lattice, this results in shifts over the folded ${{\mathbb{R}}}^{2}/{{\mathcal{L}}}_{2}$ space generated by τ₁(ξ) and τ₂(ξ).

Now, the only information the observer needs in order to redraw ${\mathcal{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 ${{\mathcal{L}}}_{1}$. Then, using the coordinates (ξ₁, ξ₂) together with the given twist angle, the observer re-traces the primitive cell of ${{\mathcal{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 ${\tau }_{a}{\mathcal{P}}$ described above is topologically ergodic.

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

$${\tau }_{i}({\boldsymbol{\xi }})\,=\,\left(\begin{array}{l}({\xi }_{1}\,+\,{a}_{1}{A}_{i1})\ {\rm{mod}}\ {a}_{1}\\ ({\xi }_{2}\,+\,{a}_{2}{A}_{i2})\ {\rm{mod}}\ {a}_{2}\end{array}\right),\quad i\,=\,1,2,$$

(13)

with:

$$A\,=\,\left(\begin{array}{ll}\frac{\sin (\beta \,-\,\theta )}{\sin (\beta )}&\frac{{a}_{1}}{{a}_{2}}\frac{\sin (\theta )}{\sin (\beta )}\\ -\frac{{a}_{2}}{{a}_{1}}\frac{\sin (\theta )}{\sin (\beta )}&\frac{\sin (\beta \,+\,\theta )}{\sin (\beta )}\end{array}\right)\quad ({\rm{Det}}\ A\,=\,1).$$

(14)

As it is well known, if at least two A_ij’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 ${\mathcal{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,

$$f({\boldsymbol{\xi }})\,=\,\mathop{\sum}\limits_{{\boldsymbol{q}}\,\in\, {{\mathbb{Z}}}^{{\boldsymbol{2}}}}{f}_{{\boldsymbol{q}}}\,{\left({e}^{\imath 2\pi {\xi }_{1}/{a}_{1}}\right)}^{{q}_{1}}{\left({e}^{\imath 2\pi {\xi }_{2}/{a}_{2}}\right)}^{{q}_{2}},$$

(15)

the algebra of T_f operators has two generators T_j corresponding to the elementary functions:

$$f\,\mapsto\, {u}_{j}({\boldsymbol{\xi }})\,=\,{e}^{\imath 2\pi {\xi }_{j}/{a}_{j}},\quad j\,=\,1,2.$$

(16)

Furthermore, since $({u}_{j}\circ {\tau }_{i})({\boldsymbol{\xi }})\,=\,{e}^{\imath 2\pi {A}_{ij}}{u}_{j}(\xi )$, the CR’s (12) become ${S}_{i}{T}_{j}\,=\,{e}^{-\imath 2\pi {A}_{ij}}{T}_{j}{S}_{i}$. As such, the algebra which contains all Galilean invariant dynamical matrices is generated by four unitary elements: U₁ = S₁, U₂ = S₂, U₃ = T₁, U₄ = T₂, with CR’s ${U}_{i}{U}_{j}\,=\,{e}^{\imath 2\pi {\phi }_{ij}}{U}_{j}{U}_{i}$, where

$${{\Phi }}\,=\,\left[{\phi }_{ij}\right]\,=\,\left(\begin{array}{llll}0&0&-{A}_{11}&-{A}_{12}\\ 0&0&-{A}_{21}&-{A}_{22}\\ {A}_{11}&{A}_{21}&0&0\\ {A}_{12}&{A}_{22}&0&0\end{array}\right)$$

(17)

Let us point out that, by considering two degrees of freedom per resonator, additional entries of the Φ-matrix can be populated by $\pm \frac{1}{2}$ values, as explained in²² for a spinful model.

Predictions via K-Theory

In K-Theory⁷², the set of projections that can be (stably) deformed into a given projection is called the K₀-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₀-class. These classes of projections can be added and subtracted, hence they form an abelian group, the K₀-group of the algebra. Two homotopic projection P and $P^{\prime}$ are also similar: $P^{\prime} \,=\,{U}^{* }PU$ for some unitary element from the same algebra. If ${\mathcal{T}}$ is a trace on the algebra, then automatically ${\mathcal{T}}(P)\,=\,{\mathcal{T}}(P^{\prime} )$ 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 noncommutative 4-torus, the K₀-group has 2^d−1-generators {e_J}, conveniently labeled by a subset J ⊆ {1, 2, 3, 4} of directions with ∣J∣ = even (the empty set ${{\emptyset}}$ is also a valid choice for J). Furthermore, for generic Φ-matrices, the noncommutative 4-torus accepts a unique trace ${\mathcal{T}}$, which coincides with the trace per volume⁶¹. Any gap projection P_G of a dynamical matrix defines a K₀-class and accepts a decomposition in terms of the generators ${[{P}_{G}]}_{0}\,=\,{\sum }_{J}{n}_{J}\ {[{e}_{J}]}_{0}$. The integer numbers n_J are called gap labels⁶² 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,

$${\mathcal{T}}{[{P}_{G}]}_{0}\,=\,\mathop{\sum}\limits_{J}{n}_{J}\ {\mathcal{T}}{[{e}_{J}]}_{0}.$$

(18)

On the other hand,

$${\mathcal{T}}({P}_{G})\,=\,\mathop{{\rm{lim}}}\limits_{| {{\mathcal{L}}}_{1}| \,\to\, \infty }\frac{{\rm{Tr}}({P}_{G})}{| {{\mathcal{L}}}_{1}| }\,=\,\mathop{{\rm{lim}}}\limits_{| {{\mathcal{L}}}_{1}| \,\to\, \infty }\frac{\{\#\,{\rm{states}}\,{\rm{below}}\,G\}}{| {{\mathcal{L}}}_{1}| },$$

(19)

hence ${\mathcal{T}}({P}_{G})\,=\,{\rm{IDS}}(G)$. As such, if we can resolve the values of the trace on the generators of the K₀-group, we can make a prediction about the allowed values of IDS. For the noncommutative torus⁷³, this extremely useful piece of information was supplied in⁷⁴, and we have

$${\rm{IDS}}(G)\,=\,{\mathcal{T}}({P}_{G})\,=\,\mathop{\sum }\limits_{J\,\subset\, \{1,2,3,4\}}^{| J| \,=\,{\rm{even}}}{n}_{J}\ {\rm{Pf}}({{{\Phi }}}_{J}),$$

(20)

where Φ_J is the matrix Φ restricted to indices J and Pf is the pfaffian of the resulting antisymmetric matrix. In our case, this gives the prediction (technically, the Pfaffians should be positive but relaxing the signs makes the fitting easier. For bulk-boundary correspondence, however, one need to reassess the signs of the terms)

$${\rm{IDS(G)}}= \, {n}_{{{\emptyset}}}\,+\,{n}_{\{1,3\}}{A}_{11}\,+\,{n}_{\{1,4\}}{A}_{12}\,+\,{n}_{\{2,3\}}{A}_{21}\\ \,+{n}_{\{2,4\}}{A}_{22}\,+\,{n}_{\{1,2,3,4\}}\ {\rm{Det}}(A).$$

(21)

When there are no linear relations with integer coefficients between A_ij-s, we can compute all topological invariants supplied by n_{i, j}-s by fitting (21) to the numerically obtained IDS curves in Figs. 2c, 3c, 4c. Unfortunately, ${\rm{Det}}(A)\,=\,1$, hence we can only determine the sum ${n}_{{{\emptyset}}}\,+\,{n}_{\{1,2,3,4\}}$ via this procedure.

The values of the Chern numbers on the K₀-generators were computed in⁷⁵[p. 141]:

$${{\rm{Ch}}}_{J^{\prime} }{[{e}_{J}]}_{0}\,=\,\left\{\begin{array}{l}0\,{\rm{if}}\,J^{\prime} \,\nsubseteq\, J,\hfill\\ 1\,{\rm{if}}\,J^{\prime} \,=\,J,\hfill\\ {\rm{Pf}}({{{\Phi }}}_{J\,\setminus\, J^{\prime} })\,{\rm{if}}\,J^{\prime} \,\subset\, J,\end{array}\right.\quad J,J^{\prime} \,\subset\, \{1,2,3,4\}.$$

(22)

Since the Chern numbers are also linear maps, their values on the gap projection ${[{P}_{G}]}_{0}\,=\,{\sum }_{J}{n}_{J}\ {[{e}_{J}]}_{0}$ can be straightforwardly computed from (22):

$${{\rm{Ch}}}_{J^{\prime} }{[{P}_{G}]}_{0}\,=\,{n}_{J^{\prime} }\,+\,\mathop{\sum}\limits_{J^{\prime} \subsetneq J}{n}_{J}\ {\rm{Pf}}({{{\Phi }}}_{J\,\setminus\, J^{\prime} }).$$

(23)

As one can see, the top Chern number corresponding to $J^{\prime} \,=\,\{1,2,3,4\}$, also known as the second Chern number and denoted by Ch₂, is always an integer, but the lower Chern numbers may not be. We will use the above relations in our discussion of the bulk-boundary correspondence.

Theory meets numerics

The predictions based on Eq. (21) are reported in Fig. 6a–c. Specifically, for the bilayer analyzed in Fig. 2, there are linear dependencies and the prediction from Eq. (21) reduces to

$${\rm{IDS}}(G)\,=\,({n}_{{{\emptyset}}}\,+\,{{\rm{Ch}}}_{2}({P}_{G}))\, + \, ({n}_{\{1,3\}}\,+\,{n}_{\{2,4\}})\cos \theta \\ + \, ({n}_{\{1,4\}}\,-\,{n}_{\{2,3\}})\sin \theta .$$

(24)

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. 2c. Fittings of the IDS curves inside the eight large gaps identified in Fig. 2b are reported in Fig. 6a, while the topological invariants extracted from the fittings are reported in Table 1. 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.

Fig. 6: Integrated density of states fitting for gap labeling.

Integrated density of states (color maps) as a function of twist angle θ for lattices with (a) a₁ = a₂ = 1 and β = π/2, (b) a₁ = a₂ = 1 and $\beta \,=\,\pi /\sqrt{7}$, and (c) ${a}_{1}\,=\,1/\root {[4]}\of {2}$, a₂ = 1 and $\beta \,=\,\pi /\sqrt{7}$. The fittings according to the theory are overlaid as dashed lines with colors corresponding to different spectral gaps, illustrating excellent agreement.

Table 1 Gap labeling for β = π/2, a₁ = a₂ = 1 lattice.

For the bilayer analyzed in Fig. 3, one additional linearly independent term is present in the IDS expression:

$${\rm{IDS}}(G)= \, ({n}_{{{\emptyset}}}\,+\,{{\rm{Ch}}}_{2}({P}_{G}))\,+\,{n}_{\{1,3\}}\frac{\sin (\beta \,-\,\theta )}{\sin \beta }\\ \, +({n}_{\{1,4\}}\,-\,{n}_{\{2,3\}})\frac{\sin \theta }{\sin \beta }\,+\,{n}_{\{2,4\}}\frac{\sin (\beta \,+\,\theta )}{\sin \beta }.$$

(25)

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. 3c. Fittings of the IDS curves inside the eight large gaps identified in Fig. 3b are reported in Fig. 6b, while the topological invariants extracted from the fitting are reported in Table 2. The symmetry of the spectral butterfly can be used to automatically fit eight additional IDS curves in the right side of the IDS plot.

Table 2 Gap labeling for $\beta \,=\,\pi /\sqrt{7},{a}_{1}\,=\,{a}_{2}\,=\,1$ lattice.

Finally, for the bilayer analyzed in Fig. 4, we have five linearly independent terms present in the IDS expression:

$${\rm{IDS}}(G)=\, ({n}_{{{\emptyset}}}\,+\,{{\rm{Ch}}}_{2}({P}_{G}))\,+\,{n}_{\{1,3\}}\frac{\sin (\beta \,-\,\theta )}{\sin \beta }\hfill\\ +{n}_{\{2,4\}}\frac{\sin (\beta \,+\,\theta )}{\sin \beta }\,+\,{n}_{\{1,4\}}\frac{{a}_{1}}{{a}_{2}}\frac{\sin \theta }{\sin \beta }\,-\,{n}_{\{2,3\}}\frac{{a}_{2}}{{a}_{1}}\frac{\sin \theta }{\sin \beta }.$$

(26)

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. 4c. Fittings of the IDS curves inside the eight large gaps identified in Fig. 4b are reported in Fig. 6c, while the topological invariants extracted from the fitting are reported in Table 3.

Table 3 Gap labeling for $\beta \,=\,\pi /\sqrt{7},{a}_{1}\,=\,\sqrt[4]{2},{a}_{2}\,=\,1$ lattice.

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.

Bulk-boundary correspondence and existence of 2nd-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 ${\widehat{D}}_{k}({\boldsymbol{\xi }})$. The bulk-boundary for class A in higher dimensions states⁷⁵ 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⁷⁵[p. 175],

$${{\rm{Ch}}}_{\{k,i\}}({P}_{G})\,=\,{\left.\frac{{N}_{ki}}{\sqrt{| {{\mathcal{L}}}_{1}| }}\right|}_{| {{\mathcal{L}}}_{1}| \,\to\, \infty },\quad k\,\in\, \{1,2\},\,i\,\in\, \{3,4\},$$

(27)

where N_ki is the net number of eigenvalues of ${\widehat{D}}_{k}({\xi }_{1},{\xi }_{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 (the eigenvalues which cross the reference point from below/above are counted with ± signs, respectively). Using (22), we can write the bulk-boundary principle explicitly,

$${\left.\frac{{N}_{ki}}{\sqrt{| {{\mathcal{L}}}_{1}| }}\right|}_{| {{\mathcal{L}}}_{1}| \,\to\, \infty }\,=\,{n}_{\{1,2,3,4\}}\ {\rm{Pf}}\left({{{\Phi }}}_{\{1,2,3,4\}\setminus \{k,i\}}\right)\,+\,{n}_{\{k,i\}}.$$

(28)

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. 7a–d, which was simulated on a 21 × 21 lattice, we analyze the bulk-boundary correspondence for the orange -gap (V) from Fig. 4b. As one can see, ${\widehat{D}}_{1}({\boldsymbol{\xi }})$ displays 21 positively sloped chiral bands when ξ₂ is varied, and no chiral bands are present for the other three cases (the small gaps are due to the hybridization of the boundary modes from opposite sides and should be neglected). This is consistent with n_{1, 4} = 1 and trivial values for the other invariants. Similarly, in Fig. 7e, f, which was simulated on a 26 × 26 lattice, we analyze the bulk-boundary correspondence for the purple -gap (III) from Fig. 4b. In this case, ${\widehat{D}}_{1}({\boldsymbol{\xi }})$ displays 27 positively sloped chiral bands when ξ₁ 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.

Fig. 7: Observation of edge states and bulk-boundary correspondence for lattice with $\beta \,=\,\pi /\sqrt{{\bf{7}}},{{\bf{a}}}_{{\bf{1}}}\,=\,{\bf{1}}/\root {[4]}\of {{\bf{2}}},{{\bf{a}}}_{{\bf{2}}}\,=\,{\bf{1}}$.

Eigenvalues (Ω²) as a function of normalized phason coordinates ξ₁/a₁, ξ₂/a₂ for 100 × 100 lattice with periodic boundary conditions (black curves), or smaller lattices (red curves) with periodic-open (${\widehat{D}}_{2}({\xi }_{1},{\xi }_{2})$) and open-periodic (${\widehat{D}}_{1}({\xi }_{1},{\xi }_{2})$) boundaries along the a₁ and a₂ directions. In particular, (a–d) focus on the orange gap (V) using a 21 × 21 lattice with θ = 1.55, while (e–h) focus on the purple gap (III) using a 26 × 26 lattice with θ = 2.75.

Additional boundary spectra are reported in Fig. 8, which were simulated on a 23 × 23 lattice and θ = 0.2. They correspond to the red -gap (I) in Fig. 4b. From the data in Figs. 8, 6, and by assuming n_{{1, 2, 3, 4}} = − 1, we have:

$$\begin{array}{lll}&&\frac{{N}_{23}\,=\,4}{23}\,=\,0.17,\,| {{\rm{Ch}}}_{\{2,3\}}| \,=\,| \,-\,{\rm{Pf}}({{{\Phi }}}_{\{1,4\}})\,+\,{n}_{\{2,3\}}| \,=\,0.18;\\ &&\frac{{N}_{13}\,=\,2}{23}\,=\,0.09,\,| {{\rm{Ch}}}_{\{1,3\}}| \,=\,| \,-\,{\rm{Pf}}({{{\Phi }}}_{\{24\}})\,+\,{n}_{\{1,3\}}| \,=\,0.06;\\ &&\frac{{N}_{24}\,=\,2}{23}\,=\,0.09,\,| {{\rm{Ch}}}_{\{2,4\}}| \,=\,| \,-\,{\rm{Pf}}({{{\Phi }}}_{\{1,3\}})\,+\,{n}_{\{1,3\}}| \,=\,0.1;\\ &&\frac{{N}_{14}\,=\,6}{23}\,=\,0.26,\,| {{\rm{Ch}}}_{\{1,4\}}|\,=\,| \,-\,{\rm{Pf}}\left({{{\Phi }}}_{\{2,3\}}\,+\,{n}_{\{1,4\}}| \,=\,0.25\right..\end{array}$$

(29)

As one can see, the predicted bulk-boundary correspondence (28) holds with a remarkable precision given the relatively small size of the lattice (we purposely kept the lattice size small to be able to count the edge bands). Similar agreements hold true for other spectral gaps from Fig. 4b and, for example, for the cyan -gap (VII) we found n_{{1, 2, 3, 4}} = 1. As such, twistronics is capable of generating gaps with 2nd-Chern numbers.

Fig. 8: Observation of additional edge states and bulk-boundary correspondence for lattice with $\beta \,=\,\pi /\sqrt{{\bf{7}}},{{\bf{a}}}_{{\bf{1}}}\,=\,{\bf{1}}/\root{[4]}\of {{\bf{2}}},{{\bf{a}}}_{{\bf{2}}}\,=\,{\bf{1}}$.

Eigenvalues (Ω²) as a function of normalized phason coordinates ξ₁/a₁, ξ₂/a₂ for 100 × 100 lattice with periodic boundary conditions (black curves), overlaid to the eigenvalues of 23 × 23 lattice (red curves) with periodic-open (${\widehat{D}}_{2}({\xi }_{1},{\xi }_{2})$) and open-periodic (${\widehat{D}}_{1}({\xi }_{1},{\xi }_{2})$) boundaries along the a₁ and a₂ directions. The figures focuses on the red gap (I) for θ = 0.2.

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 K-theoretic machinery invented by Bellissard⁶¹ enabled us to classify all topological phases from class A supported by these twisted bilayers and to produce a high-throughput of topological gap labels. Let us mention that the topological invariants can be computed directly using the algorithms developed in⁷². However, those algorithms require a substantial computational effort and, as such, the technique based on the IDS fitting is a major outcome of our work. Using the newer results from⁷⁵, 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 ${{\mathcal{L}}}_{3}$ is added on top of ${{\mathcal{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 ${{\mathbb{R}}}^{2}/{{\mathcal{L}}}_{2}$ and ${{\mathbb{R}}}^{2}/{{\mathcal{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 noncommutative 6-torus, which hosts topological phases with 3rd-Chern numbers⁷⁵.

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.