Magic-angle semimetals

We introduce and characterize a class of `magic-angle' semimetal models and, as an application, propose multiple cold atomic quantum emulators of twisted bilayer graphene. The models are defined in one to three dimensions and all contain momentum space nodes as well as an incommensurate quasiperiodic potential. With both numerical and analytical tools, we demonstrate an undiscovered link between magic-angle phenomena and a single-particle eigenstate quantum phase transition. At criticality, we report a nonanalytic behavior of the density of states, the appearance of flatbands, and wave functions which Anderson delocalize in momentum space displaying a multifractal scaling spectrum. We outline the necessary conditions for magic-angle semimetals and construct effective Hubbard models on superlattices by computing Wannier states, which demonstrates that the effective interaction scale is dramatically enhanced across the single-particle transition. As a result, we argue that the eigenstate quantum criticality is unstable towards the inclusion of interactions. All sufficient ingredients of our proposal are available in present day cold atomic laboratories for which the magic-angle effect can be exploited to induce strong correlations in quantum degenerate gases.

Introduction. Since the discovery of quantum phases of matter, particularly strongly correlated phases, there has been a need to fundamentally understand and control them. In particular, engineering band structures with non-trivial topological wave functions has achieved success in creating and controlling novel quantum phases in a variety of systems such as doped strong spin-orbit insulators [1], ultracold atoms [2] and photonics [3] with artificial gauge fields [4], electronics with topological circuit Lagrangians [5], semiconducting multilayer heterostructures [6,7], and dynamically driven Floquet systems [8]. However, the ability to design band structures that, in the presence of interactions, reliably create unconventional superconductors through strong correlations has remained out of reach.
The discovery of a correlated insulator [9,10] and accompanying superconducting phases [11,12] in graphene heterostructures revolutionized this field of quantum engineering. New, strongly interacting, solid state systems can now be engineered with a rather weakly correlated two-dimensional semimetal (graphene). Qualitatively, the Moiré pattern at certain magic angles of twisted bilayer graphene (TBG) [13], graphene-boron nitride [14], and other van der Waals heterostructures [15] are believed to quench the kinetic energy of the electronic degrees of freedom below the scale of electronic interactions. As a consequence, 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 theoretically [16][17][18] and experimentally [19][20][21] prior to this year's groundbreaking discover- * YXF, EJK, and JHW contributed equally.
ies [9][10][11]. However, quantitative theoretical predictions in the strong coupling regime remain elusive, and there is currently no general consensus about the form of an effective low-energy description [22][23][24][25][26][27][28]. This paper presents a completely different perspective on this hotly debated problem by first distilling the basic physical phenomena that generically create correlated flat bands out of two-dimensional Dirac cones. We show that in arbitrary dimension quasiperiodicity generically, and in a universal fashion, creates flat bands in nodal, semimetallic, band structures due to a previously unnoticed single particle quantum phase transition (QPT)what we call the "magic angle" in analogy to TBG. Near this transition we demonstrate the existence of strong correlations by computing Wannier states within a series of band gaps; these lead to a Hubbard model with a dramatically increased interaction scale. We therefore argue that the single particle quantum critical state is unstable towards the inclusion of interactions, which form a Mott 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 which we call magic-angle semimetals. To demonstrate this, we classify the family of these models with symmetry protected nodes as well as introduce and solve a series of models; most of which can be straightforwardly realized with existing ultracold atom and trapped ion experimental setups. Thus, we propose a simple route to emulate the phenomena of magicangle TBG in a wide variety of quantum many body systems [15,29]. c. An infinite number of semimetal minibands form as the transition is approached. We construct exponentially localized Wannier states on the first four minibands (see Fig. 3) leading to a model with an effective, strongly renormalized Hubbard interaction U eff /t eff in terms of the bare interaction U/t. the formĤ =T +V +Û (1) containing single particle hoppingT , a quasiperiodic modulationV (such as potential scattering), and interparticle interactionsÛ . The kinetic termT has isolated nodal points in the Brillouin zone where the DOS vanishes in a power-law fashion (i.e. semimetallic). The quasiperiodicity inV is encoded in an angle originating either from twisted bilayers or the projective construction of quasicrystals [30], and it is characterized by an amplitude W and an irrational modulation Q.
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 unique "unfreezing" transition in momentum space, which represents a novel form of delocalization [31]. For a subset of magic-angle semimetals [including Eq. (2) below], the semimetal reenters at a second transition W c with a reversed sign of the helicity at each Dirac node [32]; for general Q, multiple semimetalmetal-semimetal transitions can appear as W is tuned, see Figs. 1b. Third, the quenched kinetic energy implies a divergence of the dimensionless interaction coupling constant, Fig. 1c, leading to exotic many-body states. Importantly, these effects occur generically under the necessary condition that the quasiperiodic potential respects the symmetries which protect the semimetallic touching points (see supplement [33]). Effective models. A variety of effective models illustrate our proposal. For concreteness, we formally introduce the following, a 2D tight-binding Hamiltonian of "perfect" spin-orbit coupling (SOC) on a square lattice (we fix the bare lattice spacing to unity and = 1) where σ µ are the Pauli matrices, r is the position on the lattice, c r are two-component spinors of fermionic annihilation operators (at the single particle level, they can also be bosonic), and φ µ is a phase offset typically averaged over. The kinetic partT has a momentumspace dispersion E 0 (k) = ±t sin 2 k x + sin 2 k y with four Dirac nodes and a velocity v 0 = t, see Fig. 2a inset. The role of the twist angle and interlayer hopping in TBG are replaced by the incommensurate wavevector Q and potential strength W , respectively. Unless otherwise stated, Q = 2π/ϕ 2 where ϕ is the golden ratio, and in numerical simulations we employ rational approximants, Q n ≡ 2πF n−2 /F n , where the system size L = F n is given by the nth Fibonacci number [32]. In addition to Eq. (2) 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. [32]), and a 1D long range hopping model with a power-law dispersion E = −t sign(cos k)| cos k| σ with σ < 1 (this corresponds to a hopping t ij ∼ |i − j| −(1+σ) for |i − j| 1 [34]-in this 1D case, v is not a velocity) [33]. 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 lattice [35,36] or artificial gauge fields [37][38][39], whereas the 1D model we consider can be realized using trapped ions [40]. The 3D variant of Eq. (2) is theoretically possible to implement [41][42][43], but has not been experimentally realized yet. In each of these experimental setups, quasiperiodic potentials can then be realized, e.g. by additional lasers [44], programmable potentials [45], or a digital mirror device [46]. Single-particle spectrum and velocity renormalization. We first discuss the spectral characteristics of magicangle semimetals probed through the DOS, defined as where E i is the ith eigenenergy. 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 Fig. 2a for the model Eq. (2). 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 of degenerate spin states a distance Q away in momentum space [32], see the gray lines in Fig. 2a, inset. This process "carves out" a square around each Dirac cone which contains 2[(π − Q)L/2π] 2 states. For the chosen Q, a second gap is numerically observable at larger W , see Fig. 2b, which corresponds to a fourth order Umklapp scattering and generates a miniband of 2[(4Q− 3π)L/2π] 2 states per cone. For a given incommensurate Q (along with the distance between nodes), there is an infinite sequence of relevant orders in perturbation theory that open up gaps near zero energy, forming minibands; this is in contrast to the commensurate case when this sequence is finite. For Q = 2π/ϕ 2 , it is given by half the even Fibonacci numbers F 3n /2, which is the sequence 1, 4, 7, 72, 305, . . . (see supplement [33] for proofs).
The renormalization of the velocity can be analytically determined using fourth-order perturbation theory by integrating out states at Manhattan distance Q and 2Q from each Dirac point [17,33]. In terms of the dimensionless coupling constant α = W/[2t sin(Q)] for Eq. (3) The root of the numerator captures the first magic-angle transition line well when Q > π/2 and is only inaccurate for the transition line not adiabatically connected to W = 0 at small Q, see Fig. 1b. For reentrent semimetallic phases, Eq. (3) indicates the reversal of the Berry phase, consistent with the inversion of miniband states in 3D [32].
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 non-analytic and a metallic spectrum with finite ρ(0) develops for α > α c , see Fig. 2b. In particular, for d > 1 and fixed Q, ρ(E) ∼ |W − W c | −β |E| d−1 implying the velocity v(W ) ∼ |W − W c | β/d . Surprising, we find β ≈ 2 across a broad range of models, Q values, and dimensions [32,33], in agreement with perturbation theory in 2D. This exponent also agrees with similar results in TBG [17,18] in-dicating that this phenomenon is within the same universality class. In 1D this magic-angle effect also exists but is modified by the form of the dispersion such that ρ(E) ∼ |W − W c | −β |E| 1/σ−1 , and for the case σ = 1/3 we find β = 4.0 ± 0.8. This velocity renormalization for a multitude of models is plotted in Fig. 1a, here v(W ) is determined by computing the twist dispersion of the energy using exact diagonalization and agrees with the calculation of ρ (d−1) (0), see [33].
Critical single-particle wave functions: Critical properties across the magic-angle transition appear not only in the spectrum but also in the eigenfunctions near E = 0 [32], see Figs. 2c,d. We show that magic-angle semimetals are intimately linked to the physics of Anderson transitions in momentum space, and we conjecture this to be true for TBG due to universality.
We compute the low energy wavefunctions using exact diagonalization on small system sizes and Lanczos for large L reaching up to L = 610. 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. 2h. In contrast, the form of the wave functions is completely different in the metallic state, see Fig. 2i, as it appears delocalized both in momentum and real space with non-trivial structure [33]. Finally, in the reentrant semimetal, the wave functions are again ballistic, see Fig. 2j. 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 Figs. 2c,d.
In order to quantify the eigenstate QPTs of the wave functions, we generalize the multifractal wave function analysis [31] to momentum space. We define the inverse participation ratio of the energy eigenstates in momentum space [32] ψ E (k) at a given energy E We can now apply properties of the scaling exponent τ M (q), typically used to analyze real space localization, to eigenstates in momentum space. It monotonically increases [obeying τ M (0) = −d and τ M (1) = 0] and distinguishes delocalized wave functions [τ M (q) = q(d − 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 > q c ) = 0 for a given q c ∈ (0, 1]; their peak height is system size independent, as in standard localized states, but show multifractal correlations in their tails [31]. 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 supplement [33] for details. The scaling analysis of I M (q, L), presented in Figs. 2df for Eq. (2), discloses three phases of distinct wavefunction structures in momentum space. A frozen spectrum τ M (q) occurs in the two semimetal regimes of smallest W and sufficiently large W . 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 [33], we find τ M (q) ≈ 2(q − 1) − 0.25q(q − 1) (in the region |q| < 1 and within the limits of our numerical precision) in Fig. 2f. The observation of similar behavior in all models that we investigated [33] corroborates the interpretation of the magic-angle phenomenon as one of eigenstate quantum criticality and generalizes the quasiperiodic 3D Weyl semimetal to diffusive metal QPT [32] to arbitrary dimensions. In two dimensions we do not find any signatures of diffusion (consistent with the marginality of two dimensions [33, 47,48]) and in one dimension the semimetal transitions directly to an Anderson insulator [33]. Lastly, when d > 1 and W is much larger than the magic-angle transition, all investigated models undergo Anderson localization in real space.
Wannier functions, Hubbard model.-We now turn to the interparticle interaction termÛ in the Hamiltonian in Eq. (1). In order to illustrate how the appearance of flatbands enhances correlations, we construct a series of emergent Hubbard models near the magic-angle transition for Eq. (2) at φ µ = π/2 supplemented bŷ with n rσ = c † rσ c rσ . In contrast to the previous discussion, we take systems that are not perfectly incommensurate in order to build translationally invariant Hubbard models. In particular, we still use the rational approximants Q n = 2πF n−2 /F n , only now we take the size of the system L = mF n for some integer m, effectively taking the thermodynamic limit in L before the limit of pure quasiperiodicity Q n → Q. This is reminiscent of Moiré lattices used to model TBG, and similarly, we can unambiguously define a supercell of size = F n and isolate bands in k-space. In particular, these bands are intimately related to the hierarchy of minibands derived with perturbation theory: when = F 3a+b for integers a and b = 1, 2, the gap for the lowest |E| band opens at order F 3a /2 in perturbation theory [33] (for = F 3a , the Dirac nodes gap at order F 3a /2). This allows us to systematically study each miniband opening.
To build the Hubbard models, we perform approximate joint diagonalization on the position operators (x µ ) projected (with projection operator P ) onto a given bandX MB µ ≡ Px µ P in order to determine the Wannier states [49] (for details and code, see along with dramatic enhancements of interactions, which reach up to a massive U eff /t eff ∼ 4100U/t for the fourth miniband with supercell = 377, as shown in Fig. 1c. This can also been shown analytically using a one step renormalization group calculation, which yields the di- . Due to finite size, the apparent location of W c can artificially shift, therefore in Fig. 1c we use W c =W c sin Q sin Qn whereW c is the transition point when n → ∞. In the supplement we present the data for a large set of ( , m) corroborating our findings [33].
Away from E = 0, nearly flat (semimetallic) bands can form well before the magic-angle transition with similarly large U eff /t eff , see Fig. 3a, and near the transition, multiorbital Hubbard models [50] also appear [33].
Summary and experimental implications. We introduced a class of magic-angle semimetals of which twisted bilayer graphene is the most prominent representative. Employing a combination of numerical and analytical methods we have demonstrated the general appearance of a single-particle quantum critical point at which, simultaneously, (i) the kinetic energy vanishes universally, (ii) a non-zero density of states appears at zero energy, and (iii) the wave functions display delocalization and multifractality in momentum space. Lastly, 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.
We conclude with the experimental relevance of our theoretical findings. As mentioned, all sufficient ingredients for cold atomic quantum emulation of the magic-angle phenomenon, i.e. semimetals in a quasiperiodic potential, are available in present day cold atomic laboratories. In addition to any spectroscopic measurements that probe the density of states (e.g. radiofrequency spectroscopy [51]), we propose the analysis of wavepacket dynamics as an indicator of magic-angle physics and numerically predict a non-monotonic behaviour of the wave function spreading as a function of W [33] in the parameter regime with multiple magic angles. Moreover, our work opens an entirely novel 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).
Regarding experimentally realized twisted graphene heterostructures, it has not been obvious whether incommensuration is an important ingredient. 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. However, since by definition α ∼ 1 at the magic angle, the argument of parametric suppression of high order processes is questionable and incommensurate effects are therefore in principle possible in TBG. Furthermore, analogous to the continuous phase boundary in Fig. 1b, the physics of small angles directly connects to large, incommensurate twists [52][53][54]. Of course, the amplification of interaction effects as a consequence of non-analytic DOS persists independently of strict or approximate incommensuration, but perfect eigenstate criticality involves Anderson delocalization in momentum space which only exists for perfect incommensuration. While any rational approximant leads to a rounding of the QPT (akin to finite size effects in usual transitions), our paper promotes eigenstate quantum criticality as the ultimate origin of flatbands in magic-angle semimetals, including twisted bilayer graphene.
under Grant Number W911NF-17-1-0482. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation herein. The authors acknowledge the In this section we define the models which we analyzed. The tight-binding Hamiltonians of models dubbed "perfect" spin orbit coupling (SOC) are given byT where t is the hopping matrix element, and W is the amplitude of the quasiperiodic potential. In the two-dimensional (2D) case µ = x, y and in the three-dimensional (3D) case µ = x, y, z and r takes values in the set of all lattice points on a square (cubic) lattice. The Hamiltonian for the π-flux model has the same potential term in 2D. The hopping term is modified as followsT where we choose the gauge with A x (r) = π/2 for all sites r on the square lattice, and A y (r) = −(−1) rx π/2. For the chosen gauge, periodic boundary conditions require the lattice size in x direction to be even. The spinless honeycomb lattice model is given by a Hamiltonian of the form The sum over r A is over one of the two sub-lattices, while r is over all points. The index i labels the three nearest neighbors of r A , and d i is the vector from r A to its nearest neighbor i. The vectors δ µ are a choice of each particular model and for numerics we choose The kinetic part of the Hamiltonian for the one dimensional model with power law disperion [S1] is given in momentum spacê We assume σ < 1, this expression can be readily Fourier transformed to a tight binding model with long range hopping (LRH). This yields a hopping amplitude for |i − j| 1 and Γ(x) is the Gamma function. The potential iŝ

II. SCALING OF THE DENSITY OF STATES
In this section we discuss the numerical method used in the analysis of the spectrum and the finite size effects of the method. We use the kernel polynomial method (KPM) to calculate the density of states ρ(E), which expands ρ(E) in terms of Chebyshev polynomials up to an order N c , and we use the Jackson kernel to filter out Gibbs oscillations due to the finite expansion order. To determine the velocity v, in two-dimensions for example, we fit the low-|E| asymptote ρ(E) ≈ ρ (E = 0)|E| to extract ρ (E = 0) ∼ 1/v 2 . Note that in 2D formally ρ (E = 0) is not just a single derivate due to the |E| scaling, but we use this notation to unify 2D and 3D; the latter it is simply a second derivative. For details on the KPM technique see Ref. S2. We use twisted boundary conditions and we average over random twists to reduce finite size effects. Now we discuss the effect of finite lattice size L and finite cutoff N c on ρ(E = 0) and ρ (E = 0).
As an exemplary case we present results here for the "perfect" SOC model defined in the main text. Results on the other models are similar and we also present results on the 1D model below. Figure S1 illustrates the dependence on L and N c . For smaller N c such as N c = 2 12 = 4096, ρ (0) for all choice of L ≥ 55 almost overlap for W ≤ 0.515. For N c = 2 14 = 16384, the ρ (0) data converges as a function of L only for L ≥ 144. Still, the L convergence is only valid for W ≤ 0.515. This demonstrates that the observed convergence in L is strongly dependent on N C and therefore requires studying the scaling in N C for fixed L.
When fixing L and varying N c , the semi-metal to metal transition becomes sharper as ρ(E = 0) rises more abruptly approaching a sharp step as shown in Fig. S1. This sharpening allows us to pinpoint the location of the transitions accurately, in this case we find W c = 0.525 ± 0.005 and W c = 0.551 ± 0.005. Importantly, the peak of ρ(0) does not shrink as we vary L or N c , providing strong evidence of the presence of the intermediate metallic phase. In addition, we find that ρ (0) does not saturate as we increase the expansion order, indicating within our numerical accuracy that at the transition ρ (0) diverges, similar to what was found in 3D [S3]. From the above data of ρ (0) we determine the scaling exponent β defined by ρ (0) ∼ |W c − W | −β . We use ρ (0) data obtained for N c = 2 14 and L = 144 and we extract β from a log-log fit of 1/ρ (0) versus |W − W c |, see Fig. S2.

A. Dispersion and velocity
In this section we demonstrate the identification of the kinetic velocity as obtained from the twist dispersion with the parameter entering the low-energy asymptote of the DOS. We also compare these numerical results with perturbation theory. We implement twisted boundary conditions by including a factor e iθ·r/L for each real space field located at r and twist vector θ. Each component of θ takes value in (0, 2π) and we compute the energy eigenstates E(θ) using exact diagonalization for various values of the twist θ. Such a change of boundary condition has no effect on the bulk physics, but effectively moves the origin of the finite size induced momentum grid, so that plotting the spectrum as a function of the twist shows a projection of the dispersion onto 1/Lth of the Brillouin zone. Figure S3 shows the twist dispersion for various models in one, two, and three dimensions, which clearly demonstrates the dramatic flattening of the bands at the transition. These results where obtained for system size L = 233 in 1D, L = 144 in 2D, and L = 21 in 3D. Using the twist dispersion we can estimate the velocity by fitting the lowest energy band near 0 twist to a straight line. We compare the velocity as calculated from the twist dispersion with the KPM result of the DOS and fourth order perturbation theory in Fig. S4, which all agree well. The parameter σ defined in equation (S6) determines the behavior of the dispersion relation near k = 0. This can be seen directly from the twist dispersion in Fig. S3. We present detailed results for the 1D LRH model in Fig. S5. The DOS depends on σ by ρ(E) ∼ |E| 1/σ−1 (Ref. S4), which is demonstrated in Fig. S5. In the following we present detailed results for σ = 1/3 and leave the full exploration of this 1D model for future work. Focusing on σ = 1/3 is numerically advantageous since as we approach the transition the scaling ρ(E) ∼ |W − W c | −β |E| 2 allows us to use the second derivative of the DOS ρ (0) to estimate β and we can compute ρ (0) accurately using the KPM. Notice that the power-law remains constant when varying W in the semimetal phase, showing the 1D model is also stable to a weak quasiperiodic potential. Upon approaching the transition we find ρ (0) displays a clear divergence with no sign of saturation as we increase the expansion order (see Fig. S5), similar to the 2D model we have discussed above. We find W c = 2.05 ± 0.03 and from the power-law scaling ρ (0) ∼ |W − W c | −β we extract β = 4.0 ± 0.8 for σ = 1/3, see Fig. S2. Distinct from our results in 2D and 3D the transition in the 1D model is accompanied by real space localization. To demonstrate this we calculate the IPR in real space and momentum space at zero energy. The real space IPR becomes finite and momentum space IPR vanishes near the critical W c . In addition, when the momentum space IPR goes to zero the DOS becomes non-zero demonstrating that the generation of DOS is tied with momentum space delocalization, similar to the higher dimensional models.

III. WAVEPACKET DYNAMICS AS A PROBE OF THE SINGLE PARTICLE QUANTUM PHASE TRANSITION
In this section we demonstrate that wavepacket dynamics can be used as a tool to observe the single particle quantum phase transition in the models that posses a reentrant semimetal in two dimensions. Note that in three dimensions due to the diffusive states at finite energy that dominate the dynamics this probe is not useful [S3]. We initialize a wavepacket to a single site |ψ 0 and use the KPM to time evolve the state to obtain |ψ(t) = e −iHt |ψ 0 and from this we compute the spread of the wavepacket as a function of time from δr(t) 2 ≡ ψ(t)|r 2 |ψ(t) − ψ(t)|r|ψ(t) 2 . From the long time dynamics we extract the dynamical exponent from the scaling δr(t) 2 ∼ t 2/z . We focus on the perfect SOC Hamiltonian from the main text. The KPM expands the time evolution operator in terms of Chebyshev polynomials up to an order N C , which dictates the final time that can be reached in the numerical calculations. Here we focus on a large linear system size of L = 987 and use a KPM expansion order N C = 2 12 , which allows us to time evolve the wavepacket until it spreads out across the sample. As shown in Fig. S6 we find rather unusual wavepacket dynamics which is a signature of a sequence of semimetal-metal-semimetal transitions. As a function of increasing W we find that the speed at which the wavepacket spreads out monotonically decreases until we reach the metallic phase where the dependence on W is rather weak. Then, upon reentering the semimetal phase at larger W the wavepacket spreading speeds back up. This is clearly demonstrated in the dynamical exponent z showing non-monotonic behavior in the inset. Interestingly, our estimate of z is not diffusive, consistent with the marginal nature of 2D. We expect that this wavepacket signature can be used to detect the transition in cold atom experiments.

IV. ANALYTIC RESULTS
This part of the supplement is devoted to the summary of details on analytical arguments presented in the main text.

A. Perturbative calculation of velocity renormalization
In this section we present the perturbative calculation of velocity renormalization using the language of retarded Green's functions,Ĝ and are interested in diagonal components G k,k with k = k , only (e.g. for the DOS we only need ρ(E) = −(1/π)Im k TrG k,k (E)). We define the self energy at momentum k by all diagrams which are G 0 (k, E) irreducible and write We expand about a given node K i of the dispersion T (K i + p) T (K i ) + h (Ki) (p) to leading order in p 1/a. For models which satisfy the symmetry constraints exposed in the main text (see also Sec. IV D) Σ(K i + p, E) = EΣ E + h(p)Σ p to leading order in E, p. Henceforth, we choose the energy offset such that T (K) = 0. Then, In this section we evaluate the self energy to leading and, for some models, next to leading order in powers of W and summarize them in Tab. I. A discussion of infinite order perturbation theory can be found in Sec. IV D.
To illustrate the procedure we analyze the model of 2D perfect SOC for which the states at small k with Hamiltonian H(k) = t(sin(k x )σ x + sin(k y )σ y ) tk · σ are connected to the states at k ± Qê x,y and therefore to leading order perturbation theory For the next to leading order, all states at Manhatten distance 2Q from the origin are integrated out and we obtain W t 4 4 cos(Q) + 10 cos(2Q) + 11) csc 4 (Q) sec 2 (Q) (4 − 5 cos(Q) + 6 cos(2Q)) csc(Q) 4 sec(Q) (S12) Graphic demonstration that the model of perfect SOC in 2D is a direct sum of two decoupled π flux models. The model of perfect SOC, on the left of the equality sign, is characterized by direction dependent hopping matrices. Using blue squares and red circles to depict the bipartition, hopping only connects | , ↑ with |•, ↓ , and separately | , ↓ with |•, ↑ . The hopping in y -direction is imaginary and directed (this results from the asymmetry of σy) and, in conclusion, leads to the inclusion of a flux π per plaquette. The onsite potential does note violate the described block-diagonalization. This is the origin of Eq. (4) of the main text. It turns out that the results obtained for the 2D model of perfect SOC directly apply to the π flux model. This is best graphically shown, see Fig. S7: the model of 2D perfect SOC is a direct sum of two π-flux models which in the absence of interactions completely decouple. By consequence, all single particle results obtained for model of 2D perfect SOC also hold for the π-flux model.

B. Renormalization of interactions
In this section we present an analytical estimate of the renormalization of the interaction upon projection onto certain minibands and approaching the transition, we concentrate on d > 1. Let the bare (W = 0) model in the continuum be written as (K i are various Dirac/Weyl nodes, with linear k · p Hamiltonian h (Ki) (p)) (S13) The spectrum of h (Ki) (p) has the form v 0 |p| with bare value v 0 ∼ ta and, for contact interaction (σ = 0), g {Ki} ∼ U a d , while for Coulomb interaction (σ = d − 1) g {Ki} ∝ δ K1,K2 . Perturbation theory indicates a dimensionless parameter for σ = 0 (contact interaction), (S14) S8. Divergence of contact interaction according to Eq. (S18) for the model of 2D SOC. Here, the fourth order perturbative self energy was employed and we used γ = 1/5.
We now consider the effect of integrating out high energy states and projecting onto a miniband with effective Brillouin zone size 1/a . This leads to (S16) The renormalizations Z and v/v 0 originate from scalar and matrix components of the self-energy and were calculated perturbatively above. We now first rescale p < = a a p with p ∈ (0, 1/a) and then define c . Under this rescaling, we restore the form of Eq. (S13), including its UV cut-off 1/a, but obtain the rescaling v 0 → va/a , g → g(a/a ) d−σ Z 2 . From this we obtain the final formula for renormalization of the dimensionless coupling constant Here, γ is an unknown constant of order unity which depends on details of the cut-off of the linearized theory. Importantly, the integration reduces the bare contact interaction by a factor (a/a ) d−1 , except in the closest vicinity where the vanishing velocity overtakes the reduction, see Fig. S8.

C. Relationship to number theory
In this section we show the relationship of the sequence of relevant perturbative processes with certain well known sequences from number theory. Starting from the scattering process of order l 1 = 1 we want to determine the sequence {l n } ∞ n=1 for which the l n th order momentum transfer carves out smaller minibands than the l n−1 th order. In formulae, this implies for the 2D model of perfect SOC of the main text and arbitrary incommensuration wavevector Q the condition sin 2 (l n Q) < sin 2 (l n−1 Q). We now concentrate on the specific case Q = 2π/φ 2 = π(3 − √ 5). For this situation, the defining condition on the sequence of l n is sin 2 (πl n √ 5) < sin 2 (πl n−1 √ 5). The sequence {l n } ∞ n=1 for which l n √ 5 successively approaches integers is the sequence of denominators of the leading rational approximants, i.e. the sequence of denominators of continued fraction convergents of √ 5 (OEIS ID A001076). This sequence is also half the value of the even fibonaccis l n = F 3n /2.
This sequence also connects to the formation of minibands as found with the finite size Q n = 2πF n−2 /F n . Intuitively, when F n is even (n is a multiple of 3), then at order F n /2, the Dirac nodes gap out, but then for F n+1 and F n+2 this perturbative gap must have moved to small but finite energy, forming the miniband. This motivates using Q 3n+1 and Q 3n+2 to study the effective model of successive minibands, as we do in Sec. VI.
Therefore, we offer a proof that connects this sequence to the perturbative minibands which requires the following facts about Fibonaccis • F 3m is even while F 3m+1 and F 3m+2 is odd.
To determine the order of perturbation theory where a gap is opened, we need to find an integer g n such that g n Q n = 2π × q/F n such that q is an integer closest to F n /2 modulo F n . This is accomplished by the following theorem Theorem. Let n = 3m + r for r = 0, 1, 2, then the integer g n = F 3m /2 is the smallest integer such that g n F n−2 ≡ Fn+δn 2 mod F n with integer |δ 3m | ≤ 1. In particular, δ 3m = 0, δ 3m+1 = (−1) m , and δ 3m+2 = (−1) m+1 .
Proof. We break this up into cases. Since F n−2 does not divide F n , we merely need to find a g n that satisifies the relevant cases in order to find the unique g n , as long as g n < F n .
Case 1: r = 0. For this case we can prove the above by noting that if g 3m = F 3m /2, and F n−2 is necessarily odd, then (We use N to represent an arbitrary, unimportant, integer.) with δ 3m = 0. The value g n is smallest since δ n must vary by 2 in order for the equation Fn+δn   2 to remain integer valued, and the condition |δ n | < 1 prevents that.

D. Generality of the magic-angle phenomenon -symmetry protection
In this part of the supplement we discuss the generality of our findings by highlighting the general condition for the appearance of the magic angle phenomenon, namely the stability of the semimetal at weak coupling.
We concentrate on nodes in the kinetic termT which are protected by a symmetry group G T . For example, this analysis applies to each model we have considered in 2D as well as Dirac semimetals in 3D. Note that in general G T is a subgroup of all symmetry operations of the kinetic term. Let U S T be the representation of S T ∈ G T in the (e.g. spinorial) Hilbert space, then the symmetry of the Hamiltonian implies T (k) = U † S T T (S T k)U S T . We concentrate on high symmetry points where S T K = K, ∀S T ∈ G T . Then, a non-trivial representation implies degeneracy in view of [T (K), U S T ] = 0 ∀S T ∈ G T (formally, two non-commuting U S T are needed). We further assume a group G V of spatial (point group) symmetries of the quasiperiodic backgroundV , such that Here now,Ũ S V is the representation of S V ∈ G V andê 0 is an arbitrary vector in R d .
In this section, we consider Eq. (S9) formally to all orders in perturbation theory. The semimetallic behavior persists if a) Σ(k) is hermitian and b) T (k) + Σ(k) has the same symmetry protected touching point as T (k), i.e. if Σ(k) respects the symmetries ensuring the semimetal. In view of the incommensuration, perfect resonance is formally absent to any order in perturbation theory and therefore, the decay rate 1/τ ∼ k |T k,k | 2 δ(E k − E k ) (more generally: the anti-hermitian part of the self-energy) vanishes (T k,k denotes the T-matrix). Thus a) is fulfilled and 1/τ = 0 signals the breakdown of perturbation theory (spontaneous unitarity breaking). We can then show to all orders in perturbation theory that the semimetal is stable provided G T is a subgroup of G V .
We proceed to the proof of Σ(k) = U † S T Σ(S T k)U S T under the outlined assumptions. To get a feeling, we first consider second order perturbation theory.
We compare to This expression is invariant provided the action of S T onto G V is a bijection of G V onto itself ∀S T ∈ G T , i.e. S T S V ∈ G V ∀S V ∈ G V and S T G V = G V as this allows to uniquely relable the summation index. Taking S V = 1 implies that S T ∈ G V and hence G T is a subgroup of G V . By consequence, the representation in the Hilbert space fulfills is invariant under the symmetries protecting the semimetal. We now continue with the next order Σ (4) , from there the generality of the statement becomes apparent, The exclusion S Vê 0 +S Vê0 = 0 ensures the irreducibility with respect to G k . Again we can apply an S T transformation and exploit the two conditions exposed above to relabel both S V and S V . This implies the invariance of Σ (4) . This procedure can be used to arbitrary order in perturbation theory.

V. MULTIFRACTAL ANALYSIS
In this section, we explain the procedure of extracting the multifractal exponents and the underlying physical interpretations.

A. Scaling exponent τM
We now discuss how we extract τ M (q) from inverse participation ratio (IPR). For a given wavefunction in momentum space, the finest grid is 2π/L × 2π/L. We introduce an integer binning factor B which controls the resolution of the momentum space wavefunction. The IPR constructed in this way has been dubbed the momentum space IPR [S3] and given by where k denotes the wavevector of the binned wavefunction grid andφ E (k ; B) is the binned wavefunction. The scaling behavior only holds when 1 B L/a. To extract τ M (q), we choose two consecutive values of binning FIG. S9. τR(q) with different values of W with L = 144 and 100 realizations. For W = 0.35 and W = 0.7, the multifractal spectra approach to plane wave, τR(q) = 2(q − 1), as increasing b. For W = 0.53, the spectrum is multifractal described by a nonlinear function. For W = 1.75, τR approaches to a localized spectrum, τR(q) = 0 for q > 0, as increasing b. The existing data still suffer from a strong finite size effect and do not clearly show a localized spectrum. We therefore turn to the finite size dependence of the real space wavefunction with W = 1.75, which indeed shows localization as the real space IPR is L independent in this regime. At last, we average τ M (q; B 1 , B 2 ) over phases in the quasiperiodic potential. For L = 144, we construct τ M with 100 realization. We only use 10 realization for L = 610. For real space scaling exponents (τ R ), we adopt similar procedures to extract their values.

B. Real space multifractality
As a complementary analysis, we study the multifractal spectrum of the real space wavefunctions at zero energy. We construct the corresponding real space scaling exponent spectrum τ R (q) which encodes plane wave, multifractal, and localized wavefunctions.
For W/t < 0.525, the real space wavefunction is a plane wave, characterized by the Fourier modes at Dirac points (k x , k y ) = (0, 0), (0, π), (π, 0), and (π, π). The multifractal spectrum τ R (q) is characterized by a straight line τ R (q) = 2(q − 1). Plane wave states can also be found for 0.55 < W < 1.5, the inverted semimetallic regime. For 0.525 < W/t < 0.55, the real space wavefunctions show multifractal behavior characterized by a nonlinear τ R (q) spectrum. For W/t 1.5, the real space wavefunctions become localized, τ R (q) = 0 for q > 0. We numerically compute the τ R (q) spectrum for various values of W in different phases. In Fig. S9, we compute W/t = 0.35, 0.53, 0.7, 1.75 for L = 144 and average over 100 realizations. W/t = 0.35 and W/t = 0.7 show plane wave spectrum as increasing the real space binning factor b. W/t = 0.53 demonstrate a weakly nonlinear multifractal spectrum with a fitting function τ R (q) ≈ 2(q − 1) − 0.16q(q − 1) for |q| < 1. A much larger system size is needed for quantitatively determining the spectrum which is beyond the scope of this work. For W/t = 1.75, the spectrum gradually approaches to localized like behavior as increasing b. In our L = 144 data, we do not find a clear localized spectrum, due to finite size effects. We therefore turn to the real space IPR data as a function of system size in Fig. S10. This indeed shows strong evidence for localization as the IPR is L independent for W/t ≥ 1.75.

C. Interpretation of unfreezing transition
As we have demonstrated in the maintext, the momentum space wavefunctions serve as a proxy for the semimetalmetal transition. The multifractal analysis in momentum space here is distinct from the conventional notion of wavefunctions in real space [S5]. We discuss the interpretation of the multifractal spectrum in depth here.
For W = 0, the zero energy wavefunction is composed of the Fourier modes at the Dirac points . The zero energy states are linear combinations of these four plane waves. The probability distributions (integrating over the spin degrees of freedom) of the momentum space wavefunction generically contains four peaks, which we call "ballistic peaks". In the multifractal τ M (q) spectrum, these ballistic peaks correspond to a frozen spectrum, τ M (q) = 0 for q > 1.
The frozen feature in the momentum space wavefunction suggests that the ballistic peaks are sharply defined, with the finest localization length π/L. The τ M (q)s extracted from the numerics weakly depend on the choice of binning factors B 1 and B 2 . As W increases, satellite peaks with weaker amplitude arise. When W ≈ W c , the τ M (q)s still show freezing behavior for B = 1, 2 but become generically non-zero for larger binning factors. Such features suggests that the distance between a ballistic peak and the nearest satellite peaks is around 2 * 2π/L to 4 * 2π/L. The ballistic peaks start to hybridize with the satellite peaks when W ≥ W c . This corresponds to an unfreezing transition in momentum space, which describes a zero-measure set to an extensive set of Fourier modes. The W c determined here coincides well with the semimetal-metal transition extracted from the density of states calculations.
For larger values of W , an inverse semimetal transition, metal-to-semimetal, takes place. The same multifractal analysis in momentum space also applies.

VI. WANNIER STATES AND BUILDING HUBBARD MODELS
In order to analytically build the Hubbard models analyzed in the main text, we use the Hamiltonian as perviously definedĤ =T +V +Û , cos(Q n r µ + φ µ )c † r c r , U = U r c † r,↑ c † r,↓ c r,↓ c r,↑ .

(S30)
We take Q n = 2πF n−2 /F n and the system size to be L = mF n such that F n defines a supercell of size = F n . The model has a particle-hole symmetry when φ µ = π/2, which we will concentrate on in order to isolate bands near E = 0. The first task is to isolate bands. With L = m , our Brillioun Zone will have m 2 sampled points. Due to the imposed discrete translational symmetry, one can write the single particle HamiltonianĤ sp =T +V aŝ for n = (n x , n y ) and, using Bloch's theorem, the individual Hamiltonians are (in first quantized notation) H n = r,µ it 2 e −i2πnµ/L |n, r n, r +μ| ⊗ σ µ + h.c. + W r,µ cos(Q n r µ + φ µ ) |n, r n, r| , (S32) on a system size of size . In the thermodynamic limit m → ∞, and we can plot dispersions to see if a gap has opened. Once diagonalized, we have a set of states {E n,j }, and looking near E = 0, we can form projectors onto energy states within a minibandP = En,j ∈Miniband |n, j n, j| .
This projector should have an integer multiple of m 2 states within it, and each (n x , n y ) pair should contribute the same number of states. This is a check to determine if we have a "good band." Further, we can get an idea of where the Wannier centers should be by looking at the integrated local density of states ρ Band (r) = r|P |r . (S34) Considering systems with periodic boundary conditions, we consider the position operators {cos(2πx/L), sin(2πx/L), cos(2πŷ/L), sin(2πŷ/L)}, and we construct the projected operators {P cos(2πx/L)P ,P sin(2πx/L)P ,P cos(2πŷ/L)P ,P sin(2πŷ/L)P }. (S35) These operators no longer commute with one another and cannot be simultaneously diagonalized. The solution is approximate joint diagonalization (AJD) achieved by minimizing a cost function: for a set of operators O, we need to find a unitary matrix U such that FIG. S12. The interaction of the center band approaching the transition. Each successive pair of odd represent the opening of a new miniband. The data labeled "3 bands" comes from the mulitorbital models, following the Wannier states continuously connected to the single-orbital model. Note that Q = 2π/ϕ 2 for golden ratio ϕ and Q = 2πFn−2/ with = Fn.

B. Finite energy bands and multi-orbital Hubbard model
It is interesting to note that at finite energy, flat bands develop and in some cases intersect the band near E = 0 which is the focus of this text. As a clear example, note for = 13, see that a flat, Dirac bands appear around the region labeled 1 in Fig. S13(upper left) and intersect the band in the region labeled 3. These flat bands have greatly increased U eff /t eff away from the transition. When the multiorbital Hubbard model appears, the dispersion changes drastically, see Fig. S13(upper right) and has Wannier centers appearing along diagonals within a supercell as seen in Fig. S13(middle left). With the computed data, an entire translationally invariant Hubbard model can be constructed and the result leads to hoppings as indicated in Fig. S13(middle right), from the banded Hamiltonian Fig. S13(lower