Introduction

The engineering of band structures with non-trivial topological wave functions has achieved success in creating and controlling quantum phases in a variety of systems such as doped semiconductors1,2,3,4, ultracold atoms5,6, and metamaterials7,8. With the recent advance in twisted graphene heterostructores9,10,11,12 (i.e., “twistronics”) new, strongly interacting, solid state systems can now also be engineered with a rather weakly correlated two-dimensional semimetal (graphene)13,14,15. In these systems, as a consequence of the quenched kinetic energy, correlations dominate the physics and exotic many-body states may form. This interpretation relies on the reduction of the electronic velocity and large increase of the density of states (DOS) which was shown in twisted bilayer graphene (TBG) theoretically16,17,18,19 and experimentally20,21,22 prior to the recent groundbreaking discoveries in refs. 9,10,11. Understanding the essential single-particle ingredients necessary to build emulators of TBG can help shed light on the strong coupling regime where consensus about the form of an effective low-energy description remains elusive23,24,25,26,27,28,29.

In this manuscript, we develop the theory for twistronic emulators by first distilling the basic physical phenomena that create correlated flat bands out of two-dimensional Dirac cones. Generically, quasiperiodicity that respects the symmetry protecting the Dirac nodes creates flat bands in nodal, semimetallic band structures in a universal fashion near a previously unnoticed single-particle quantum phase transition (QPT)—what we call the “magic-angle” in analogy to TBG. At small angles in TBG, a single-scattering wavevector accounts for the majority of the band flattening17,30 but misses any QPT. With quasiperiodicity, an infinite sequence of higher wavevectors (i.e., Brillouin zone downfoldings) further flatten the bands and culminate into a QPT. This band flattening occurs irrespective of the topology present, and in fact, many of the models we study have topology distinct from TBG24. We demonstrate strong correlations by computing Wannier states within this series of bands; these lead to a Hubbard model with a quenched kinetic energy and relative to this, the interaction scale is increased dramatically. We therefore argue that the single-particle quantum critical state is unstable toward the inclusion of interactions, which form a correlated insulator at half filling.

Crucially, our findings are independent of many of the system’s details and, therefore, demonstrate the existence of a wide multitude of engineered, strongly coupled quantum systems that we call magic-angle semimetals. To demonstrate this, we classify the family of these models with symmetry protected nodes (including chiral TBG at moderate twist angles) as well as introduce and solve a series of models; most of which can be straightforwardly realized with existing ultracold atom, trapped ion, and metamaterial experimental setups. Thus, we propose a simple route to emulate the phenomena of magic-angle TBG in a wide variety of quantum many body systems15,31. Last, we show that the magic-angle effect can be observed at experimentally relevant time scales and temperatures in interacting ultracold Fermi gases through measurements of wave packet dynamics.

Results

“Magic-angle semimetals”

The whole class of magic-angle semimetal models are governed by Hamiltonians of the form:

$$\hat{H}=\hat{T}+\hat{V}+\hat{U},$$
(1)

containing single-particle hopping \(\hat{T}\), a quasiperiodic modulation \(\hat{V}\) (such as potential scattering or interlayer tunneling), and interparticle interactions \(\hat{U}\). The kinetic term \(\hat{T}\) has isolated nodal points in the Brillouin zone where the DOS vanishes in a power law fashion (i.e., semimetallic). The quasiperiodicity in \(\hat{V}\) is encoded in an angle originating either from twisted bilayers or the projective construction of quasicrystals32, and it is characterized by an amplitude W and an incommensurate modulation Q (or twist angle θ).

Generalizing the physics of the first magic angle of TBG to magic-angle semimetals results in the phenomena summarized by Fig. 1. First, increasing W quenches the kinetic energy, reducing the Dirac velocity v until it ultimately reaches zero at the single-particle quantum critical point (where the DOS becomes nonanalytic). The velocity vanishes in a universal manner characterized by critical exponents that are distinct in each dimension. Second, the DOS and wave functions display a transition from a ballistic semimetal to a metallic phase; this is a so-called “unfreezing” transition in momentum space, which represents a non-standard form of delocalization33. For a subset of magic-angle semimetals (including Eqs. (2) and (3)), the semimetal reenters at a second transition \(W^{\prime}_{c}\) with a reversed sign of the helicity at each Dirac node34; for general Q (or θ), multiple semimetal-metal-semimetal transitions can appear as W is tuned, see Fig. 1b, c. Third, the quenched kinetic energy implies a divergence of the dimensionless interaction coupling constant, Fig. 1d, leading to exotic many-body states. Importantly, these effects occur generically under the necessary condition that the quasiperiodic modulation respects the symmetries which protect the semimetallic touching points (see Supplementary Note 3).

Fig. 1: Magic-angle transition.
figure 1

A quasiperiodic potential or tunneling generically drives an eigenstate quantum phase transition from a semimetal (SM) to metal (M). a For many models, the velocity at the Dirac node v decreases with the strength of the potential W until it reaches v = 0 at the transition, Wc; this is an indication of the flattening of the bands. In some cases an intermediate metallic phase (see inset) separates a reentrant semimetal with a reversed helicity (depicted by the Dirac cones). b, c We construct a phase diagram in terms of potential strength W (interlayer tunneling for cTBG) and quasiperiodic modulation Q (twist angle θ for cTBG) by computing the density of states at zero energy ρ(0); analytical perturbative results (see Eq. (4), Supplementary Notes 3 and refs. 17,30) are represented by the green dashed lines. Cuts along the dashed white lines are presented in Fig. 2c, d. Color bars represent ρ(0) and with widths b 5, and c 1.25 and dark purple represents the value 0 on both. d An infinite number of semimetal minibands form as the transition is approached; each has higher effective interaction than the last as we approach the transition. For 2D SOC, we construct exponentially localized Wannier states on the first four minibands (see Fig. 4) leading to a model with an effective, strongly renormalized Hubbard interaction Ueff/teff in terms of the bare interaction U/t.

Full size image

Effective models

A variety of effective models (defined in Supplementary Note 1) illustrate our proposal. Here, we focus on two models: a 2D tight-binding Hamiltonian of “perfect” spin-orbit coupling (SOC) on a square lattice and a lattice model of TBG at moderate twist angles (θ ≈ 9∘) in the chiral limit (cTBG) that disallows interlayer tunneling between equivalent sub-lattices30 (we fix the bare lattice spacing to unity and ħ = 1). Nonetheless, our main conclusions also apply to TBG beyond the chiral limit for similar twist angles. (Here, we consider the chiral limit of TBG as it provides the clearest presentation of magic-angle criticality but such a transition can also be shown to persist in the full TBG model. This study will appear elsewhere). The SOC model is given by a hopping \({\hat{T}}_{{\rm{SOC}}}=t/2{\sum }_{{\bf{r}},\mu }(i{c}_{{\bf{r}}}^{\dagger }{\sigma }_{\mu }{c}_{{\bf{r}}+\hat{\mu }}+{\rm{h}}.{\rm{c}}.)\) and a quasiperiodic potential

$${\hat{V}}_{{\rm{SOC}}}=W\sum _{{\bf{r}},\mu ={\bf{x}},{\bf{y}}}\cos (Q{r}_{\mu }+{\phi }_{\mu }){c}_{{\bf{r}}}^{\dagger }{c}_{{\bf{r}}},$$
(2)

where the σμ are Pauli matrices, cr are two-component annihilation operators, t is the hopping strength, and ϕμ is the offset of the origin. The lattice model that captures the low-energy theory of cTBG at incommensurate twist angles contains \({\hat{T}}_{{\rm{cTBG}}}\) that describes nearest neighbor hopping (amplitude t = 2.8 eV) on the honeycomb lattice. The interlayer tunneling in the chiral limit is given by:

$$\begin{array}{*{20}{l}}{\hat{V}}_{{\rm{cTBG}}}=W{\sum \limits_{{\bf{r}},\mu }}\left[\cos ({{\bf{q}}}_{\mu }\cdot \displaystyle\tfrac{{\bf{r}} \,+ \, {{\bf{r}}}_{\mu }}{2}+{\phi }_{\mu }){c}_{1A{{\bf{r}}}_{\mu }}^{\dagger }{c}_{2B{\bf{r}}}-\right.\\ \left.\mathop{\sum }\limits_{n=1}^{6}\frac{{(-1)}^{n}}{3\sqrt{3}}\sin ({{\bf{q}}}_{\mu }\cdot \frac{{\bf{r}} \, + \, {{\bf{r}}}_{\mu n}}{2}+{\phi }_{\mu }){c}_{1A{{\bf{r}}}_{\mu n}}^{\dagger }{c}_{2B{\bf{r}}}+(A\leftrightarrow B)\right]+{\rm{h}}.{\rm{c}}.\end{array},$$
(3)

where clA/Br annihilates an electron on layer l, sublattice A/B, and position r. The index μ = 1, 2, 3 labels nearest neighbors such that r1 − r = (0, 1) [q1 = kθ(0, −1)] with rμ − r [qμ] being 120∘ rotations of the previous vector. The positions rμn = rμ + an where \({{\bf{a}}}_{1}=(\sqrt{3}/2,3/2)\) and each subsequent an is a 60∘ rotation of the last (i.e., labeling nearest neighbors on the triangular Bravais lattice). Last, \(| {{\bf{q}}}_{\mu }| ={k}_{\theta }=\frac{8\pi }{3\sqrt{3}}\sin (\theta /2)\) encodes the twist angle, and ∑μϕμ = 0 to satisfy C6 symmetry. Typically the offsets ϕμ in either model are averaged over. The kinetic part \({\hat{T}}_{{\rm{SOC}}}\) (\({\hat{T}}_{{\rm{cTBG}}}\)) has a momentum space dispersion with four (two) Dirac nodes and a velocity v0 = t (v0 = 3t/2), see Fig. 2a, b inset. Returning to Eq. (2), we see that Q replaces the role of the twist angle in Eq. (3); unless otherwise stated, we highlight incommensurate effects by taking Q = 2π/φ2 (\(\theta =2\arcsin (\sqrt{3}/2{\varphi }^{5})\approx 8.9{6}^{\circ }\)) where φ is the golden ratio, and in numerical simulations we employ rational approximants Qn ≡ 2πFn−2/Fn (kθ is approximated using continued fractions, see Supplementary Note 1 for details) where the system size L = Fn is given by the nth Fibonacci number34. Other values, in particular smaller π—Q and θ, are discussed in the Supplementary Information and below. In the low-energy approximation this model is identical to the continuum model studied in ref. 30 where exact flat bands are uncovered and explained; this makes this model ideal to study incommensurate effects on the lattice.

Fig. 2: Eigenstate transition as manifested in the single-particle spectrum.
figure 2

a, b DOS ρ(E) in units of \({(t{L}^{2})}^{-1}\) averaged over 300 realizations of phases ϕμ and random twisted boundary conditions (explained in more detail in Supplementary Note 2). The gray shading represents the number of states in the first miniband and matches the area of the mini-Brillouin zones around each Dirac point produced by the leading scattering vectors depicted in the inset of a, b (we chose a rhombic representation of the Brillouin zone of TBG such that k = k1G1 + k2G2 for reciprocal lattice vectors G1,2 of graphene). c, d Cuts along the dashed white lines of the phase diagram in Fig. 1b, c, displaying ρ(0) and \({{\mathcal{I}}}_{M}(q=2,L)\) [Eq. (5)]. These illustrate sequences of semimetallic and metallic transitions concomitant with momentum space delocalization (see Fig. 3). The twist dispersions illustrate the difference between semimetallic phases (e, i, f, j) and the metallic phase (g, h) as well as the remarkably reduced bandwidths (note the reduced scale). The 2D SOC (cTBG) data were obtained for Q = 2πFn−2/Fn (\(\theta =2\arcsin (\sqrt{3}{F}_{n-5}/[2{F}_{n}])\)) at L = 144 (L = 377) and KPM expansion order NC = 212 (NC = 213) in the calculation of the DOS while L = 233 in e–j.

Full size image

In addition to Eqs. (2) and (3) we have studied a multitude of other d-dimensional models in an incommensurate potential: the π-flux model and the honeycomb model in 2D, a 3D variant of Eq. (2) (studied previously in ref. 34), and a 1D long range hopping model with a power law dispersion \(E=-t\ \,\text{sign}\,(\cos k)| \cos k{| }^{\sigma }\) with σ < 135—in this 1D case, v is not a velocity (details on 1D case can be found in the last part of Supplementary Note 2). Each of these models generates flat bands and magic-angle physics similar to TBG. Importantly, these semimetallic 2D Dirac points have been realized in cold atomic setups using either a honeycomb optical lattice36,37 or artificial gauge fields38,39,40, whereas the 1D model we consider can be realized using trapped ions41. The 3D variant of Eq. (2) is theoretically possible to implement42,43,44, but has not been experimentally realized yet. In each of these experimental setups, quasiperiodic potentials can then be realized, e.g., by additional lasers45, programmable potentials46, or a digital mirror device47. Alternative emulators of Dirac semimetals can also be realized in metamaterials, e.g., in topolectrical circuits7 or in arrays of electromagnetic microwave resonators48. Quasiperiodicity can then be encoded through the spatial modulation of the electrical circuit elements.

Single-particle spectrum and velocity renormalization

We first discuss the spectral characteristics of magic-angle semimetals probed through the DOS, defined as ρ(E) = 1/NH∑iδ(E − Ei) where Ei is the ith eigenenergy and NH is the size of the single-particle Hilbert space. At weak quasiperiodic modulation, the semimetal is stable, i.e., ρ(E) vanishes at zero energy with the same power law as in the limit of W = 0, while hard spectral gaps and Van Hove singularities develop at finite energy. For Weyl and Dirac Hamiltonians the low-∣E∣ DOS obeys ρ(E) ~ v−d∣E∣d−1, and as W increases, the (d − 1)st derivative of the DOS [ρ(d−1)(0) ∝ 1/vd] increases, see Fig. 2a, b for the model in Eqs. (2) and (3), respectively. These weak coupling features may be understood at the level of perturbation theory.

We find that gaps appear at finite energy due to the hybridization around Dirac nodes a distance Q (or \(\sqrt{3}{k}_{\theta }\)) away in momentum space, see the insets in Fig. 2a, b, inset. For the SOC (cTBG) model, this process “carves out” a square (hexagon) around each Dirac cone which contains 2[(π − Q)L/2π]2 (\(4{[3\sqrt{3}{k}_{\theta }L/4\pi ]}^{2}\)) states. For a given incommensurate Q or θ, there is an infinite sequence of relevant orders in perturbation theory that produce quasi-resonances and open up gaps near zero energy, forming minibands; this is in contrast to the commensurate case when this sequence is finite. For example, for 2D SOC and Q = 2π/φ2, the infinite sequence is given by half the even Fibonacci numbers F3n/2, which is the sequence 1, 4, 7, 72, 305, … (see Supplementary Note 3). In the incommensurate limit, near the magic-angle transition this sequence of gaps produces a corresponding sequence of minibands, shown in Fig. 1d for the second, third, and fourth. We explore the effect of this sequence of minibands using superlattices in “Commensurate superlattices and Hubbard models”.

Similar to TBG, the renormalization of the velocity in the 2D SOC model can be analytically determined using fourth-order perturbation theory (details in Supplementary Note 3)17. In terms of the dimensionless coupling constant \(\alpha =W/[2t\sin (Q)]\) for Eq. (2) this yields:

$$\frac{v(W)}{v(0)}=\frac{1-2{\alpha }^{2}[1-\cos (Q)]+{\alpha }^{4}\frac{4-5\cos (Q)+6\cos (2Q)}{\cos (Q)}}{1+4{\alpha }^{2}+{\alpha }^{4}\{16+{[2+1/\cos (Q)]}^{2}\}}.$$
(4)

The root of the numerator captures the first magic-angle transition line well when Q > π/2, see Fig. 1b, independently of whether Q is commensurate or incommensurate. To describe additional magic-angles, as observed in our numerical data in Fig. 1b, c, higher order perturbation theory is required. For reentrant semimetallic phases, Eq. (4) indicates the reversal of the Berry phase, consistent with the inversion of miniband states in 3D34. In each model we have considered for d > 1, we have found that the perturbative expression for the velocity (summarized in Supplementary Table 1) has a magic-angle condition where the velocity vanishes.

As the magic-angle is approached, higher perturbative corrections become relevant. To go beyond perturbation theory, we compute the DOS using the numerically exact kernel polynomial method (KPM), on sufficiently large system sizes across a range of models of various dimensions. At a critical α = αc ~ 1 the DOS becomes nonanalytic and a metallic spectrum with finite ρ(0) develops for α > αc, see Fig. 2c, d (for cTBG \(\alpha =\frac{W}{2t\sin (3{k}_{\theta }/4)}\)). In particular, for d > 1 and fixed Q or θ, ρ(E) ~ ∣W − Wc∣−β∣E∣d − 1 implying the velocity v(W) ~ ∣W − Wc∣β/d. Surprisingly, we find β ≈ 2 in each model investigated above 1D (see Supplementary Note 2)34, indicating that this exponent is universal. In 1D, this magic-angle effect also exists but is modified by the form of the dispersion such that ρ(E) ~ ∣W − Wc∣−β∣E∣1/σ−1, and for the case σ = 1/3 we find β = 4.0 ± 0.8.

This velocity renormalization is the manifestation of the aforementioned reconfiguration of the band structure and the appearance of a sequence of minibands. Of course, broken translational symmetry precludes a standard band structure of dispersive Bloch waves. In Fig. 2e–j, we therefore illustrate this reconfigured band structure, at a fixed rational approximant, in the form of the twist dispersion (obtained by exact diagonalization in the presence of twisted boundary conditions) along high symmetry lines for the models defined in Eqs. (2), (3). We performed the analogous analysis for a multitude of models and plotted the velocity v(W) near the semimetallic touching points in Fig. 1a. The velocity v(W) as determined by computing the twist dispersion agrees with the calculation of ρ(d − 1)(0), see Supplementary Note 2.

Critical single-particle wave functions

Magic-angle semimetals are intimately linked to the physics of Anderson transitions in momentum space; this is captured by the eigenfunctions near the Dirac node energy, E = 034.

We compute the low-energy wave functions using Lanczos for large L reaching up to L = 377 and 610 in the cTBG and SOC models, respectively. Qualitatively, we find that the structure of the wave functions in the semimetallic phase is stable and adiabatically connected to the ballistic W = 0 limit, with isolated ballistic spikes in momentum space, see Fig. 3a, b. In contrast, the form of the wave functions is completely different in the metallic state, see Fig. 3c, d, as it appears delocalized both in momentum and real space with non-trivial structure (see details in Supplementary Note 5). Finally, in the reentrant semimetal, the wave functions are again ballistic, see Fig. 3e, f. Crucially, in all models that we studied, the positions of the transitions in the spectral properties of the DOS coincide with the transitions of the wave functions characteristics within numerical resolution, see Fig. 2c, d.

Fig. 3: Eigenstate transition as manifested in momentum space wave functions at the Dirac node energy E = 0.
figure 3

a–f Wave function characteristics as described by the scaling exponent τM(q) averaged over 100 random phases and twisted boundary conditions. For W < Wc and \(W\, > \,{W}_{c}^{\prime}\) the wave functions are ballistic [with a frozen τM(q)] while for \({W}_{c}\, < \, W \, < \, {W}_{c}^{\prime}\) they are critical in momentum space [τM(q) is weakly non-linear in q]. Inset of a–f corresponding momentum space wave functions. The 2D SOC (cTBG) data were obtained for Q = 2πFn−2/Fn (\(\theta =2\arcsin (\sqrt{3}{F}_{n-5}/[2{F}_{n}])\)) at L = 144 (L = 377).

Full size image

In order to quantify the eigenstate QPTs of the wave functions, we generalize the multifractal wave function analysis33 to momentum space. We define the inverse participation ratio of the energy eigenstates in momentum space34 ψE(k) at a given energy E:

$${{\mathcal{I}}}_{M}(q,L)=\sum _{{\bf{k}}}| {\psi }_{E}({\bf{k}}){| }^{2q}\ \sim {L}^{-{\tau }_{M}(q)}.$$
(5)

We can now apply properties of the scaling exponent τM(q), typically used to analyze real space localization, to momentum space. It monotonically increases [obeying τM(0) = −d and τM(1) = 0] and distinguishes delocalized wave functions [τM(q) = d(q − 1)] from exponentially localized peaks [τM(q > 0) = 0] and critical states with non-linear “multifractal” τM(q). A variant of multifractal states, which are called “frozen,” display τM(q > qc) = 0 for a given qc ∈ (0, 1]; their peak height is system size independent, as in standard localized states, but show multifractal correlations in their tails33. We employ the standard binning technique (varying the binning size B) to numerically extract the scaling exponents τM(q) in systems of a given finite size, see Supplementary Note 5 for details.

Focusing on q = 2, as shown in Fig. 2c, d for the SOC and cTBG models, respectively, the momentum space wave function at the Dirac node energy delocalizes upon crossing the magic-angle in the incommensurate limit. The momentum space delocalization can not occur in the commensurate case; Bloch’s theorem ensures the existence of states with well defined (i.e., well localized) crystalline momenta. For example, consider Eq. (2) in the commensurate limit where Q/2π = a/b (a and b are coprime integers). In this case, \({{\mathcal{I}}}_{M}(q,L)\) is bounded from below by 1/bd(q − 1) and hence τM(q) = 0 in the thermodynamic limit L/b → ∞ preventing momentum space delocalization (see Supplementary Note 5). In contrast, we here numerically access the incommensurate limit using finite size scaling of rational approximants corresponding to L = b → ∞.

The scaling analysis of \({{\mathcal{I}}}_{M}(q,L)\) at the energy of the Dirac node E = 0, presented in Fig. 3a–f for Eqs. (2) and (3), demonstrates three phases of distinct wave function structures in momentum space. A frozen spectrum τM(q) occurs in the two semimetal regimes. In sharp contrast, the function τM(q) unfreezes in the metallic phase with finite ρ(0). Surprisingly, throughout the metallic phase the spectra appear to be weakly multifractal in both momentum and real space (Supplementary Note 5), we find for the SOC model that τM(q) ≈ 2(q − 1) − 0.25q(q − 1) and for the cTBG model we obtain τM(q) ≈ 2(q − 1) − 0.15q(q − 1) (in the region ∣q∣ < 1 and within the limits of our numerical precision) in Fig. 3c, d, which are both non-linear in q. The observation of similar behavior in all models that we investigated (as listed in Supplementary Note 1) corroborates the interpretation of the magic-angle phenomenon in the incommensurate limit as one of eigenstate quantum criticality and generalizes the quasiperiodic 3D Weyl semimetal-to-diffusive metal QPT34 to arbitrary dimensions. In two dimensions, we do not find any signatures of diffusion (consistent with the marginality of two dimensions49,50) and in one dimension the semimetal transitions directly to an Anderson insulator (shown in Supplementary Note 2). Lastly, when d > 1 and W is substantially larger than the magic-angle transition, all investigated models undergo Anderson localization in real space (e.g., at W > 1.75t in the case of the 2D SOC model at Q = 2π/φ2).

Commensurate superlattices and Hubbard models

So far, our analysis regarded non-interacting magic-angle semimetals in the strict incommensurate limit. We now turn to the interparticle interaction term \(\hat{U}\) in the Hamiltonian in Eq. (1) in commensurate superlattices. In order to illustrate how the appearance of flat bands enhances correlations, we construct a series of emergent Hubbard models near the magic-angle transition for Eq. (2) at ϕμ = π/2 supplemented by:

$${\hat{U}}_{SOC}=U{\sum \limits_{{\bf{r}}}}{n}_{{\bf{r}},\uparrow }{n}_{{\bf{r}},\downarrow },$$
(6)

with \({n}_{{\bf{r}}\sigma }={c}_{{\bf{r}}\sigma }^{\dagger }{c}_{{\bf{r}}\sigma }\). In contrast to the previous discussion, we take commensurate approximations in order to build translationally invariant Hubbard models. In particular, we still use the rational approximants Qn = 2πFn−2/Fn, only now we take the size of the system L = mFn for some integer m, effectively taking the thermodynamic limit in L before the limit of quasiperiodicity Qn → Q. This is reminiscent of moiré lattices used to model TBG, and similarly, we can unambiguously define a supercell of size ℓ = Fn and isolate bands in k-space.

In particular, these bands are intimately related to the hierarchy of minibands derived with perturbation theory: when ℓ = F3a + b for integers a and b = 1, 2, the gap for the central band opens at order F3a/2 in perturbation theory (for ℓ = F3a, the Dirac nodes gap at order F3a/2. See Supplementary Note 3 for details). The series of superlattices indicated by ℓ correspond to the sequence of gap openings in “Single-particle spectrum and velocity renormalization”—making the notion precise—with (downfolded) Brillioun zones depicted in Fig. 4b. Near Wc, hard gaps open and the minibands form as illustrated in Fig. 4a for ℓ = 13, 55, 233 (respectively, the 2nd, 3rd, and 4th minibands). We conjecture that all of these minibands (as ℓ → ∞) achieve gaps near Wc as evidenced by Fig. 4a, c in the incommensurate limit, indicating something akin to the singular continuous spectrum of the Aubry–André model at criticality51. Furthermore, the central band becomes flatter, as indicated by the reduction in bandwidth seen in Fig. 4c which we track until the dispersion loses its semimetallic character.

Fig. 4: Supercell analysis and Wannier functions.
figure 4

The color coding matched across a–c (and Fig. 1d) indicates the second (orange), third (maroon), and fourth (purple) minibands. a The dispersion of Eq. (2) in the mini-Brillouin zone for superlattices (ℓ, W) = (13, 0.5), (ℓ, W) = (55, 0.5244), (ℓ, W) = (233, 0.5244) (from top to bottom); this illustrates successive emergence of minibands (from top-to-bottom) as a consequence of consecutive downfoldings of the Brillouin zone. b The corresponding mini-Brillouin zones (logarithmic scale). c The dramatic reduction in bandwidth near the critical point for each miniband. d For (ℓ, W) = (13, 0.5) and L = 104, computed Wannier function ψ(x, y) that is sitting upon the local density of states ρband(r) = ∑n∣〈r∣En〉∣2 (shown as a density plot) for eigenstates of the (orange) band \(\left|{E}_{n}\right\rangle\), on a 104 × 104 lattice. (Inset). The exponential localization of the Wannier state.

Full size image

We exploit this miniband formation and the existence of hard gaps to build symmetric Wannier functions in the semimetallic regime, see Fig. 4d. To build the Hubbard models, we perform approximate joint diagonalization on the position operators (\({\hat{x}}_{\mu }\)) projected (with projection operator P) onto a given band \({\hat{X}}_{\mu }^{{\rm{MB}}}\equiv P{\hat{x}}_{\mu }P\) in order to determine the Wannier states52 (for details and code, see Supplementary Note 6). We have checked that not only are the computed Wannier states exponentially localized to numerical precision (Fig. 4d, inset), but that they are also symmetric. Therefore, the minibands formed from the SOC model and pictured in Fig. 4 are not topological53, fragile54,55 or otherwise.

As a clear example, when W = 0.5t and (ℓ, m) = (13, 8), we see a clear band around E = 0 in Fig. 4a, and we find Wannier centers in a well-defined grid (Fig. 4d, main panel) corresponding to exponentially localized Wannier states (Fig. 4d, inset). The projected Hamiltonian has the approximate form of Eqs. (2) and (6) with a renormalized Ueff, teff and Weff = 0. With this approach, we can identify successive gaps leading up to the metallic transition from either side along with dramatic enhancements of interactions, which reach up to a massive Ueff/teff ~ 4100U/t for the fourth miniband with supercell ℓ = 377, as shown in Fig. 1d. This can also been shown analytically using a one-step renormalization group calculation, which yields the divergence Ueff/teff ~ U(1/ℓ)d−1Z2/v ~ 1/∣W − Wc∣, (\(\sqrt{Z}\) is the wave function renormalization), as shown in detail in Supplementary Note 3. Due to finite size, the apparent location of Wc can artificially shift, therefore in Fig. 1d we use \({W}_{c}={\tilde{W}}_{c}\frac{\sin Q}{\sin {Q}_{n}}\) where \({\tilde{W}}_{c}\) is the transition point when n → ∞. In Supplementary Fig. 16, we present the data for a large set of (ℓ, m) corroborating our findings.

Away from E = 0, nearly flat (semimetallic) bands can form well before the magic-angle transition with similarly large Ueff/teff, see Fig. 4a. In very close proximity to the transition, multi-orbital Hubbard models appear (see Supplementary Note 6).

Experimental cold atomic realization

All sufficient ingredients for emulating magic-angle phenomenon are available in ultracold atomic gas and metamaterial48,56 experiments. In particular for ultracold atomic gases, to probe fermionic strong correlations, the atomic species 40K and 6Li are prime candidates; we estimate that the underlying lattice can be relatively shallow (around eight lattice recoil energies), and need temperatures relative to the Fermi temperature (of the entire gas) T/TF ≈ 0.25 to ensure fermion population fills but does not exceed the first miniband. To see large correlations, trap sizes should accommodate at least roughly 30 × 30 optical lattice sites. In addition to any spectroscopic measurements that probe the DOS (e.g., radiofrequency spectroscopy57), we propose and demonstrate (in more detail in Supplementary Note 4) that the analysis of wave packet dynamics is an indicator of magic-angle physics. In the absence of interactions, we numerically predict a non-monotonic spreading of the wave function for increasing W (see Supplementary Fig. 10) in the regime with multiple magic angles. We have also studied the interacting model in the hydrodynamic regime by using a generalization of the Boltzmann kinetic equation58 (see details in Supplementary Note 4). Its solution confirms the drastic decrease of the expansion velocity and a substantial enhancement of diffusive dynamics near the magic angle, see Fig. 5. The observation of these effects is possible within experimentally realistic observation time of 50t−1 (~10–100 ms). Moreover, our work demonstrates an experimental protocol for realizing strong correlations by first cooling the gas to quantum degeneracy and then applying a quasiperiodic potential to create flat bands without the need to cool the system in a Mott insulator phase or load the atoms into a flat band.

Fig. 5: Boltzmann wave packet spreading.
figure 5

Spreading of the mean square radius 〈r2〉 = ∑rr2ρ(r) of the particle density ρ(r) as a function of time in units of the inverse hopping rate 1/t (a α < αc, b α > αc). Here, we consider the interacting 2D SOC model, Eqs. (2) and (6), and we employ Eq. (4) to incorporate the magic-angle effect (occurring at αc ≈ 0.53 in this approximation) into a semianalytical hydrodynamic treatment. The initial steady state at finite temperature is defined by a particle [energy] density \(\rho ({\bf{r}})={e}^{-{{\bf{r}}}^{2}/[2{\xi }^{2}]}/{\xi }^{2}\) [\({\rho }_{E}({\bf{r}})={v}_{0}\left(1+3{e}^{-{{\bf{r}}}^{2}/{\xi }^{2}}\right)/{\xi }^{3}\)], with v0 ≡ v(α = 0) is the bare velocity and we chose ξ = 4 for the initial spread of the density profile. The hydrodynamic equations were numerically solved in the presence of an onsite repulsion U(α = 0) = 0.025t and Umklapp scattering rate 1/τ(α = 0) = 0.0075t (see Supplementary Note 4).

Full size image

Discussion

In summary, we introduced a class of magic-angle semimetals and demonstrated the general appearance of a single-particle QPT in the incommensurate limit at which, simultaneously, (1) the kinetic energy vanishes universally, (2) a non-zero DOS appears at zero energy, and (3) the wave functions display delocalization and multifractality in momentum space. In the presence of interactions we demonstrated that this eigenstate criticality leads to a strongly correlated Hubbard model by computing Wannier states on a superlattice. Lastly, we presented a detailed discussion of an experimental realization in cold atomic quantum emulators.

Regarding experimentally realized twisted graphene heterostructures at much smaller twist angles than we have considered here (θ ≈ 1.1∘), it has not been obvious whether incommensuration is an important ingredient59. Quasiperiodic effects rely upon weakly detuned processes at which the total transferred momentum wraps the Brillouin zone. In contrast, the momentum transfer induced by scattering off a small angle superstructure is minute. Therefore—it is often concluded—both effects of incommensurability and intervalley scattering are negligible as processes in higher order perturbation theory. As our numerics demonstrate, this results in the suppression of the width of the metallic sliver in Fig. 1b, c that scales like \(W^{\prime}_{c} -{W}_{c} \sim {\theta }^{3}\), making observing such a metallic phase exceedingly difficult at small twist angles. Nonetheless, we expect Anderson delocalization in momentum space even at small twist angles. This is because this physics is dominated by rare resonances (as manifested in the locator expansion60) and controlled by α, while perturbative processes are parameterized by W/t and are therefore small. Furthermore, the contiguous phase boundary in Fig. 1b, c may imply that the physics of small angles directly connects to large, incommensurate twists61,62,63. However, within present day numerics, we cannot exclude that this boundary of eigenstate phase transitions terminates at a finite, small angle, which would imply the existence of a critical Anderson delocalization end point in Fig. 1b, c. The coexistence of finite DOS with other features of this phase at larger angles suggests that the phase extends to θ → 0 (Q → π), but an end point is appealing as it would establish a theoretical paradigm of quasiperiodic Anderson tricriticality. Any rational approximant or commensurate angle truncates the infinite sequence of resonances and minibands which leads to a rounding of the QPT (akin to finite size effects in usual transitions) and the absence of momentum space delocalization. The amplified interactions due to flat bands and an enhanced DOS occur for both incommensurate and commensurate cases as Fig. 1d demonstrates. This enhancement coupled with eigenstate quantum criticality in the incommensurate limit characterizes magic-angle semimetals, including TBG at moderate twist angles.

Note Added: While this manuscript was under consideration following its announcement in arXiv:1809.04604, independent proposals to simulate twistronics in cold atoms appeared and were published in refs. 64,65.

Methods

Numerical methods

The numerical methods used are the KPM66 for the DOS, exact diagonalization and Lanczos for eigenstates, approximate joint diagonalization for Wannier functions67, and numerical partial differential equation solvers for the Boltzmann kinetics. These methods are explained in context in the Results section with additional details in the Supplementary information, particularly Supplementary Notes 2, 4, 5, and 6.