Moire Localization in Two Dimensional Quasi-Periodic Systems

TL;DR

This paper investigates a two-dimensional quasi-periodic system with Moire potential, revealing critical scaling behavior and multiple mobility edges, and proposes an optical lattice experimental scheme for realization.

Contribution

It demonstrates that 2D Moire systems share the same scaling exponent as the Aubry-Andre model and exhibit numerous mobility edges, expanding understanding of quasi-periodic localization phenomena.

Findings

Scaling exponent $ u \,\approx\, 1$ in 2D Moire systems

Presence of infinitely many mobility edges

Experimental scheme using arbitrary wavelength lasers

Abstract

We discuss a two-dimensional system under the perturbation of a Moire potential, which takes the same geometry and lattice constant as the underlying lattices but mismatches up to relative rotation. Such a self-dual model belongs to the orthogonal class of a quasi-periodic system whose features have been evasive in previous studies. We find that such systems enjoy the same scaling exponent as the one-dimensional Aubry-Andre model $ \nu\approx 1 $, which saturates the Harris bound $ \nu>2/d=1 $ in two-dimensions. Meanwhile, there exist an infinite number of mobility edges different from the typical one-dimensional situation where only a few or no mobility edges show up. An experimental scheme based on optical lattices is discussed. It allows for using lasers of arbitrary wavelengths and therefore is more applicable than the one-dimensional situations requiring laser wavelengths close to certain incommensurate ratios.

Tables

Table 1. Table 1: The critical behavior for self-dual, orthogonal class of models in a quasi-periodic potential.

Dimension	Mobility edge?	Exponent $\nu$	Harris bound satisfied?
1 Aubry and André (1980)	No	1	No
2	Yes	$\approx 1$	Yes (saturate)
3 Devakul and Huse (2017)	Yes	$\approx 1.6$	Yes (well above)

Equations

$$H_{\text{bi}}=\sum _{⟨i,j⟩}(a_i^{†}{a}_j+b_i^{†}{b}_j)+\sum _{i,j}{V}_{ij}{n}_i^A{n}_j^B,$$

$$H_{\text{eff}}=\sum _{⟨i,j⟩}{a}_i^{†}{a}_j+\sum _i{\mu }_i^A{n}_i^A+(A\leftrightarrow B),$$

$$H$$

$${\mu }_i$$

$$u_i$$

$${\phi }_1$$

$${\theta }_{\text{c}}=\frac{m\pi }{2}\pm \arccos \frac{q^2-p^2}{q^2+p^2},m\in Z$$

$${\mathbf{T}}_1=q\hat{x}-p\hat{y},$$

$${\mathbf{T}}_1=\frac{q+p}{2}\hat{x}+\frac{q-p}{2}\hat{y},$$

$$H$$

$$-\frac{V_d}{2}\sum _{x,y}{\mu }_{x,y}{c}_{x,y}^{†}{c}_{x,y},$$

$${\mu }_{x,y}$$

$$u_{x,y}$$

$${\phi }_1$$

$$-({\varphi }_{x+1,y}+{\varphi }_{x-1,y}+{\varphi }_{x,y+1}+{\varphi }_{x,y-1})-\frac{V_d}{2}{\mu }_{x,y}{\varphi }_{x,y}=E{\varphi }_{x,y},$$

$${\varphi }_{x,y}=e^{2\pi i(x{\phi }_2+y{\phi }_1)}\sum _{m,n=-\infty }^{\infty }{\psi }_{m,n}{e}^{2\pi i\left[n(u_{x,y}+{\phi }_1)+m(v_{x,y}+{\phi }_2)\right]},$$

$$-2{\mu }_{m,n}{\psi }_{m,n}-\frac{V_d}{4}({\psi }_{m+1,n}+{\psi }_{m-1,n}+{\psi }_{m,n+1}+{\psi }_{m,n-1})$$

$$=E{\psi }_{m,n}.$$

$$V_d^{(D)}=4.$$

$${\text{IPR}}_m=\sum _{x,y=1}^L|{\varphi }_{x,y}^{(m)}{|}^4,$$

$$\xi \sim |V_d-V_d^{(c)}{|}^{-\nu }$$

$${\tilde{\alpha }}_0$$

$$A_{a,b}$$

$${\tilde{\alpha }}_0=f(L/\xi )=g\left((V_d-V_d^{(c)})L^{1/\nu }\right),$$

$${\tilde{\mu }}_{x,y}$$

$${\tilde{u}}_{x,y}$$

$$m_1\hat{x}+\gamma m_2\hat{y}=n_1{\hat{x}}^{\prime }+\gamma n_2{\hat{y}}^{\prime },$$

$$\left(\begin{array}{c}{m}_1\\ \gamma m_2\end{array}\right)=\left(\begin{array}{cc}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{array}\right)\left(\begin{array}{c}{n}_1\\ \gamma n_2\end{array}\right).$$

$$\cos \theta =k_1/k_3,$$

$$k_1^2+k_2^2=k_3^2.$$

$$k_1=q^2-p^2,$$

$$S=\frac{1}{\sqrt{q^2+p^2}}\left(\begin{array}{cc}q& p\\ -p& q\end{array}\right)=\frac{q}{\sqrt{q^2+p^2}}\mathbb{I}+\frac{p}{\sqrt{q^2+p^2}}i{\sigma }_y,$$

$$S^{-1}=S^T,$$

$$A=\frac{1}{q^2+p^2}\left(\begin{array}{cc}{q}^2-p^2& -2qp\\ 2qp& q^2-p^2\end{array}\right)=\frac{q^2-p^2}{q^2+p^2}\mathbb{I}-\frac{2qp}{q^2+p^2}i{\sigma }_y$$

$$\left(\begin{array}{c}{n}_1\\ \gamma n_2\end{array}\right)$$

$$\left(\begin{array}{c}{m}_1\\ \gamma m_2\end{array}\right)$$

$${\mathbf{T}}_1=\left(\begin{array}{c}q\\ -p\end{array}\right),$$

Full text

I Introduction

Overlaying lattices with mismatches leads to the so-called Moire pattern, which has long been an intriguing subject in two-dimensional heterostructures Castro Neto et al. (2009); Dean et al. (2013); Hunt et al. (2013). Usually, a commensurate twist angle between two layers are considered, such that the bilayer system still possesses an enlarged periodicity. The interlayer hopping and interactions would then induce a reconstruction of the original Bloch waves in each layer into multiple bands. That provides a new knob to tune the property of systems with various twisting angles, leading to novel examples including the flat-band induced superconductivity for bilayer graphene twisted at magic angle Cao et al. (2018); Yankowitz et al. (2019). More recently, the experimental advancement has made it possible to stabilize the bilayer graphene system at incommensurate large twist angle Yao et al. (2018). In such situations, the combined bilayer system breaks any translation symmetry and the Bloch waves in each layer are destroyed completely. Then, it is of interest to ask what may be the general phenomena expected in these circumstances.

The incommensurately twisted bilayer system resembles quasi-crystals, in the sense that only rotation but not translation symmetry is preserved therein. In the case of quasi-crystals, localization of particles has been a focus of study for decades Sokoloff (1986); Krajci and Hafner (1995); Villa et al. (2006); Park et al. . One could relate the two systems using the following schematic reasoning. Consider a bilayer, two-dimensional system with interlayer interactions,

[TABLE]

where $a_{i}^{\dagger},b_{i}^{\dagger}$ create particles in two layers respectively, and $n_{i}^{A}=a_{i}^{\dagger}a_{i},n_{i}^{B}=b_{i}^{\dagger}b_{i}$. Now, suppose one could take as a starting point the decomposition $\mu_{i}^{A}=\sum_{j}V_{ij}\langle n_{j}^{B}\rangle$, and $\mu_{j}^{B}=\sum_{i}V_{ij}\langle n_{i}^{A}\rangle$, then the above model reduces to

[TABLE]

which describes a decoupled bilayer under onsite chemical potentials $\mu_{i}^{A,B}$ respectively. With reasonable interactions, i.e. not infinite-range interaction with identical strength, one would expect an incommensurate twist angle to result in an incommensurate, quasi-periodic pattern for $\mu_{i}^{A,B}$. That gives an analogous scenario as in quasi-crystals, with the difference that here, it is the chemical potentials rather than lattice site locations that break translation symmetry.

The inter-sample interaction induced slow-relaxation (“quasi many-body localization”) has previously been demonstrated for a one-dimensional translation invariant ladder by Yao et. al. Yao et al. (2016). Specifically, given a highly non-uniform initial density distribution and strong interactions, the particles in one chain would serve as random chemical potential for the other, similar to the case mentioned above. But in contrast, due to the incommensurate twist angles breaking translation symmetry here, it is expected a genuine localization without initial state dependence would occur for incommensurate Moire bilayers. In particular, even if the particles are uniformly distributed in one layer, they could still function as disorder potentials for the other through quasi-periodic $V_{ij}$.

Thus, we are motivated to investigate the “Moire localization” possibly held in a two-dimensional quasi-periodic system. As a first step, we focus on an effective scenario in Eq. (2), namely, a regular single-layer lattice under the perturbation of quasi-periodic potentials. This would serve as a good starting point for further taking into account fluctuations in Eq. (1). Also, although the localization of quasi-periodic systems in one-dimension have been thoroughly studied for decades Nandkishore and Huse (2015); Abanin et al. ; Khemani et al. (2017); Roati et al. (2008); Schreiber et al. (2015); Choi et al. (2016), since the renowned work by Aubry and Andre Aubry and André (1980) in 1980, its generalization to higher dimensions has just started quite recently Bordia et al. (2017); Devakul and Huse (2017); Sutradhar et al. . The theoretical studies on this regard have remained chiefly in the single-particle level, with many crucial aspects still open to discussions. In particular, the scaling exponents for orthogonal class of models in two-dimensions have been left out in Ref. Devakul and Huse (2017); Sutradhar et al. due to various difficulties. But this class of models are most relevant for experiments on genuine two-dimensional quasi-periodic system with real hopping, and also for the theoretical analysis based on Eqs. (1) and (2). Therefore, it is useful to clarify the relevant single-particle physics before considering the full many-body interactions.

Apart from preparing for analysis of interaction-induced localization, generalizing the paradigmatic Aubry-Andre model to higher dimensions itself may yield interesting results. As shown by Devakul and Huse in Ref. Devakul and Huse (2017), for orthogonal class of models in three dimensions, the scaling exponents for (single-particle) localization transitions are the same for both quasi-periodic and purely random potentials. This is in contrast to the situation in one-dimension, where the quasi-periodic Aubry-Andre model gives rise to the scaling exponent strongly violating the Harris bound Harris (1974) formulated for purely random models. Such a violation in one-dimension is indicated in Ref. Khemani et al. (2017) to persist into the many-body localization scenario. Experimentally, a higher-dimensional quasi-periodic potential is actually more suitable for cold atom experiment. This is because the property of the system is most sensitive to the relative rotation angle between the Moire perturbing potential and the main lattice potential, rather than the relative lattice constants. That means one could use lasers of any frequency available in the laboratory to engineer the desirable system. Finally, having a good understanding of the single particle physics for quasi-periodic systems would pave the way for possible perturbative treatment Basko et al. (2006) of many-body localization therein in the future.

In this work, we study a two-dimensional square lattice under the perturbation of Moire-type of potential, which takes exactly the same lattice geometry and constants as the underlying lattice. A commensuration condition for the relative rotation angles between the two lattice potentials is provided. Such a condition also helps identifying “better” rotation angles that are further from a commensurate one, resembling the “better” ratio of lattice constants in a one-dimensional quasi-periodic system, that is usually taken to be the golden ratio. Two important features of such a system is revealed. First, for the inverse-participation-ratios in the eigenstate-disorder plane, in contrast to the one-dimensional situation where only a few step-like jumps exist, our two-dimensional model exhibits a rapid and continuous change. Such a character indicates the unusual mobility edges in two dimensions. Second, the critical exponent for localization length is extracted from the multifractal analysis, with the value $\nu\approx 1$ saturating the Harris bound $\nu>2/d=1$ for two-dimensions. Finally, we provide an experimental scheme to realize such a Moire type of model, where lasers of arbitrary wavelengths could be applied.

II Models and the Commensuration Conditions

Consider a square lattice under perturbations of Moire-type of potentials with the same geometry as the underlying major one. There are three types of controlling parameters for such a perturbing potential. (a) Stretching: it could possess different lattice constants. (b) Rotation: it may be rotated by an angle $\theta$ with respect to the underlying lattice. (c) Translation: there may be a global relative translation between these two lattices. Both (a) and (b) could generate an incommensurate potential, but as in the experiment Ref. Bordia et al. (2017), without interactions, (a) alone would result in a separable $V(x,y)=V(x)+V(y)$ such that the model is reduced to two orthogonal lower-dimensional ones. Meanwhile, (c) does not change the universal properties of the system in thermodynamic limits, but only affects microscopic details for a finite size system, i.e. the wave functions in the boundary. Thus, the global rotation angle of perturbing potential is the most important factor characterizing a genuine two-dimensional quasi-periodic system. It was also found in previous works that the system’s properties are most sensitive to rotation angles of the perturbing potentials Devakul and Huse (2017); Sutradhar et al. .

Thus, we focus on the following model with two overlapping square lattices, where the stronger one gives rise to the tight-binding approximation, and the weaker one generates the Moire type of onsite chemical potentials,

[TABLE]

Here $i=(x_{i},y_{i})$ denotes the square lattice sites, with $x_{i},y_{i}\in\mathbb{Z}$. The Moire superlattice potential $\mu_{i}$, compared with the main lattice, is rotated by an angle $\theta$ at the origin $(x,y)=(0,0)$ and then translated by $(x,y)=(a,b)$. The stretching (a) would be neglected in most parts below, because it does not lead to extra features and only complicates the analysis. For instance, either stretching the lattice constant by a factor $\sqrt{2}$, or a rotation of $\pi/4$ would produce an incommensurate potential, while combining these two only gives a usual commensurate, staggered lattice. We would chiefly investigate the influence of different rotation angles $\theta$, and average over various translations $(a,b)$ (or equivalently, $\varphi_{1},\varphi_{2}$) when extracting universal properties from finite size numerical results.

The first question to ask is under what condition would the rotation angle $\theta$ result in a commensurate perturbing potential $\mu_{i}$. Similar to the analysis in honeycomb lattices Shallcross et al. (2010); Lopes dos Santos et al. (2012); Yao et al. (2018), the commensuration condition is prescribed by solving a Diophantine equation specified by lattice geometries. For square lattices, the solution is readily provided by the Pontryagin’s triples Steuding (2005) (see Appendix A for derivations)

[TABLE]

where $q,p$ are positive, coprime integers. This is the commensurate condition for Moire rotation angles of rectangular lattices with rational aspect ratios $\gamma$. Specific to the square lattice $\gamma=1$, due to four-fold rotation symmetry and mirror symmetry along $\hat{x}\pm\hat{y}$ directions, one only needs to consider $\theta_{c}=\arccos\frac{q^{2}-p^{2}}{q^{2}+p^{2}}\in(0,\pi/4)$ and $\theta_{c}^{\prime}=\pi/2-\theta_{c}=\arccos\frac{2qp}{q^{2}+p^{2}}\in(0,\pi/4)$, with coprime $q>p$. The resulting Moire pattern has periodic structures, with the new Moire-Bravais lattice vectors $\boldsymbol{T}_{1},\boldsymbol{T}_{2}$ given below. When $p$ or $q$ is an even number (they cannot both be even due to the coprime requirement),

[TABLE]

and when $(q,p)$ are both odd numbers,

[TABLE]

In both cases, the two vectors have identical length $|\boldsymbol{T}_{1}|=|\boldsymbol{T}_{2}|=\sqrt{q^{2}+p^{2}}$ (when $q$ or $p$ is even) or $\sqrt{(q^{2}+p^{2})/2}$ (when $q,p$ are both odd), and are orthogonal to each other $\boldsymbol{T}_{1}\cdot\boldsymbol{T}_{2}=0$. An example of commensurate Moire pattern is drawn in Fig. 1(a), where $q=11,p=1$, and the Bravais vectors are $\boldsymbol{T}_{1}=6\hat{x}+5\hat{y}$, $\boldsymbol{T}_{2}=5\hat{x}-6\hat{y}$. In contrast, an incommensurate example is shown in Fig. 1(b) where only the four-fold rotation, but not any translation symmetry, exists.

Here, we are interested in localization induced by incommensurate potentials, so it is worthwhile to choose among different incommensurate angles for further numerical analysis. Strictly speaking, any angle other than those given by Eq. (4), i.e. any $\theta=\arccos(A)$ with $A$ (or $A$ module $\pi/2$) being irrational numbers, would render incommensurate systems belonging to the same universal class. But when one is trying to extract information from numerical results in a finite-size system, there is a subtle question of “how incommensurate” the potentials are. Mathematically, it means if one takes a series of commensurate angles $\{\theta_{n},n=1,2,3\dots\}$, whose $n\rightarrow\infty$ gives an incommensurate angle $\theta_{\infty}$, then a “more incommensurate” angle would require each $\cos\theta_{n}$ being given by larger $(q,p)$ in Eq. (4). Equivalently put, since the length of Moire-Bravais vectors $\sim\sqrt{q^{2}+p^{2}}$, that means given a system size (and therefore an upper limit of $(q,p)$, since for larger $(q,p)$, the system cannot cover even one Moire period), “more incommensurate” $\theta_{\infty}$’s are farther away from commensurate angles. This is in the same spirit as in one-dimension, where the “golden-ratio” $(\sqrt{5}-1)/2$ — the “most irrational number” farthest away from rational ones — is usually adopted in numerics as the relative ratios of lattice constants between main and perturbing lattices.

To illustrate this point more clearly, we plot the commensurate angles for maximal $(q,p)\leq 20$ in Fig. 2. We readily note that for all small angles $\theta\rightarrow 0$ or certain large angles (i.e. $\theta=\pi/5$ or $\theta=\pi/8$), there appears less commensurate angles nearby. For small angles, there is a practical difference between bilayer heterostructure and optical lattice systems. For bilayer ones, there is a competition between interlayer attractive van der Waals force and elastic force. The former one tends to align the two layers by distorting lattices so as to minimize interlayer distance, while the latter one would favor keeping lattice shapes and therefore the interlayer twisting. For small twist angles $\theta\rightarrow 0$, the van der Waals force overwhelms and the system could always be taken as a commensurate one Woods et al. (2014), with the enlarged periodicity approximately given by $1/\theta$. Although such an issue does not exist in optical lattices where lattice shapes are fixed, for generality of the results in practice, we consider the large twist angle $\pi/5$ for most parts in the following.

It is worth emphasizing that even if one chooses an irrational lattice constant $\beta$, there still exists the need to pick out twist angles further away from “almost-commensurate” one so as to minimize finite size effects. The process there would consists of first finding a sequence of rational lattice constants close to $\beta$, and then for each member in this sequence, find the density of commensurate angles. We would not digress to such a situation in this work.

III Moire Localization

For the one-dimensional Aubry-Andre model Aubry and André (1980), its self-duality gives rise to the unique localization transition (lack of mobility edge) and the uniform localization length. But if the incommensurate potential is given by relative twists rather than stretch, the two-dimensional eigenstates cannot be decomposed into orthogonal one-dimensional components. Then we would see that such a higher-dimensional generalization breaks the unique localization length and transition, in accordance with the three-dimensional results in Ref. Devakul and Huse (2017). It can be traced back to that the Thouless formula adopted to prove these two points Aubry and André (1980), which we review in Appendix D, no longer holds for a higher-dimensional system. We would also illustrate the mobility edge clearly in numerics later on.

III.1 What does self duality imply in the Moire model?

For later convenience, we rewrite Eq. (3) as

[TABLE]

which is equivalent to Eq. (3) up to a constant. That also makes the model particle-hole symmetric after averaging over the phases $\varphi_{1},\varphi_{2}$. The single-particle eigenfunction $\phi_{x,y}=\langle x,y|\phi\rangle=\langle 0|c_{x,y}|\phi\rangle$ satisfies

[TABLE]

where $H|\phi\rangle=E|\phi\rangle$ and $|0\rangle$ is vacuum. Consider the Fourier transformation

[TABLE]

where $m,n\in\mathbb{Z}$. Then, for incommensurate rotations where $(m,n)\sin\theta$ or $(m,n)\cos\theta$ are never integers and therefore different Fourier modes do not couple, the wave functions $\psi_{m,n}$ satisfy the dual equation

[TABLE]

Compare Eqs. (8) and (10), we see that the self-duality maps $\frac{V_{d}}{2}\leftrightarrow\frac{8}{V_{d}},$ and $E\leftrightarrow\frac{4E}{V_{d}}$. Then the invariant disorder strength under duality mapping can be extracted as

[TABLE]

The dual Fourier transformation Eq. (9) has the property that if $|\psi_{m,n}|^{2}$ is a localized wave function, in the sense that $\sum_{m,n}|\psi_{m,n}|^{2}$ and $|\psi_{m,n}||_{m\text{ or }n\rightarrow\infty}$ is finite ($|\psi_{m,n}|$ is exponentially localized around certain $(m_{0},n_{0})$), we would have a delocalized $\phi_{x,y}$ as $\sum_{x,y}|\phi_{x,y}|^{2}\rightarrow\infty$ 111 Specifically, $\sum_{x,y}|\phi_{x,y}|^{2}=\sum_{x,y}\sum_{mnkl}\psi_{m,n}^{*}\psi_{k,l}A_{mnkl}^{x}B_{mnkl}^{y}C_{mnkl}$, where $A_{mnkl}^{x}=e^{2\pi ix[(l-n)\cos\theta+(k-m)\sin\theta]}$, $B_{mnkl}^{y}=e^{2\pi iy[-(l-n)\sin\theta+(k-m)\cos\theta]}$, and $C_{mnkl}=e^{2\pi i[(l-n)\varphi_{1}+(k-m)\varphi_{2}]}$. For $m\neq k$ or $n\neq l$, one could use $\sum_{x=-L}^{L}e^{2\pi i\beta x}=\sin[\pi\beta(2L+1)]/\sin[\pi\beta]$ and show that $D_{mnkl}\equiv C_{mnkl}\sum_{x,y}A_{mnkl}^{x}B_{mnkl}^{y}$ is finite as $L\rightarrow\infty$. Then, $\sum_{x,y}|\phi_{x,y}|^{2}=\sum_{m\neq k\text{ or }n\neq l}\psi_{m,n}^{*}\psi_{k,l}D_{mnkl}+\sum_{mn}|\psi_{m,n}|^{2}\sum_{x,y}1\rightarrow\infty$, where the first term is finite due to the localization of $\psi_{m,n}$. . The boundary condition (or the normalization condition) for the wave function is that the maximal value of $|\psi_{m_{0},n_{0}}|$ at the localized site $(m_{0},n_{0})$ is taken to be a constant. As such, the duality mapping implies a one-to-one correspondence between a localized eigenfunction at disorder $V_{d}$ and energy $E$, and a delocalized one at $16/V_{d}$ and energy $4E/V_{d}$. It does not specify which one is the localized solution.

So far, the analysis is completely in parallel with that for the one-dimensional Aubry-Andre model. However, to prove all eigenstates are localized for $V_{d}>V_{d}^{(D)}=4$ (and delocalized when $V_{d}<V_{d}^{(D)}$), that is, a unique localization transition for all eigenstates, Ref. Aubry and André (1980) invoked the Thouless formula Thouless (1972). Such a formula relates the localization length $\xi_{E}$ at energy $E$ with the density of states $D(\varepsilon)$ as $1/\xi_{E}=\int_{-\infty}^{\infty}(\ln|E-\varepsilon|)D(\varepsilon)d\varepsilon$. Together with the duality mapping for density of states, one could show that all eigenstates are localized at the $V_{d}>V_{d}^{(D)}$ side with uniform localization length $1/\xi=\ln(V_{d}/4)$ (see Appendix D). However, the Thouless formula, which requires the nearest-neighbor hopping in one-dimension (i.e. Eq. (39) and (40)), no longer holds in higher dimensions. Then it is expected the mobility edge would generally exist for a higher dimensional model.

In the following subsections, we perform exact diagonalization to further reveal the nature of the localization transitions in our Moire model.

III.2 Level statistics

We first compute the energy level statistics for the ratio of neighboring gaps Oganesyan and Huse (2007); Atas et al. (2013). Arranging the energy levels $\{E_{m}\}$ in the order $E_{m}<E_{m+1}$, and defining the gap between neighboring levels in a finite-size system as $\delta_{m}=E_{m+1}-E_{m}$, the ratio reads $r_{m}=\min(\delta_{m},\delta_{m+1})/\max(\delta_{m},\delta_{m+1})$. The distribution function of $r_{m}$’s, $P(r)$, would approach a Poisson limit $P(r)=2/(1+r)^{2}$ with $\langle r\rangle=\int_{0}^{1}dr(rP(r))\approx 0.386$ for a fully localized system with complete sets of integrals of motion. In contrast, for delocalized systems, the neighboring energy level repulsions would lead to the Gaussian orthogonal ensemble (GOE) with $\langle r\rangle\approx 0.5307$. Ref. Devakul and Huse (2017) computed such a quantity for a fixed system size in three-dimensions, and found an “intermediate” regime near the dual-invariant disorder strength $V_{d}^{(D)}$ with $\langle r\rangle$ approaching the Gaussian ensembles. It is of interest to see whether such an intermediate phase, which would indicate the existence of mobility edge, persist to our two-dimensional situation. Also, we check the change of $\langle r\rangle$ as the system size increases so as to indicate the behavior in the thermodynamic limit.

The qualitative behavior of $\langle r\rangle$ as a function of $V_{d}$ can be expected as follows. At $V_{d}=0$, the translationally-invariant limit, momentums are conserved and $E_{\boldsymbol{k}}=-2(\cos k_{x}+\cos k_{y})$. Thus, it can be regarded as a localized system in momentum space. Deviating away from $V_{d}=0$, the momentum conservation is immediately broken by the Moire perturbing potential, while real-space localization has not been established. Then, due to the lack of complete local integrals of motion, the level repulsion would lead to the GOE distribution for $\langle r\rangle$. As $V_{d}$ increases, more and more eigenstates are localized, and when the full localization for all eigenstates occur, $\langle r\rangle$ would once again approach the Poisson limit with real-space positions as good quantum numbers to denote localized states.

The results for $\theta=\pi/5$ shown in Fig. 3 match the intuitive reasoning given above. Here, we avoid the subtle regime with small $V_{d}\ll 1$. This is because at $V_{d}=0$, the system possesses many symmetries, such as translation, reflection, and four-fold rotation symmetries. That renders the Hamiltonian block-diagonalized into different symmetry sectors and a meaningful level statistics requires the decomposition of Hilbert space into these sectors Santos and Rigol (2010a, b); Sandvik (2010); D’Alessio and Rigol (2014). Close to $V_{d}=0$, those symmetries are only weakly broken and a severe finite-size effect is expected. We readily see from Fig. 3 that at small $V_{d}$, $\langle r\rangle$ decreases with irregular behaviors. On the other hand, the transition near $V_{d}^{(D)}=4$ is well-behaved. At $V_{d}=V_{d}^{(D)}=4$, $\langle r\rangle$ keeps approaching the GOE limit as system sizes $L$ increase. Passing $V_{d}=4$, $\langle r\rangle$ starts to drop. We note that for relatively smaller system size (i.e. $L=60,80,100$), the curves cross at around $V_{d}\approx 4.4$ with the increase of $L$. However, for larger $L$’s such a “crossing” appears to be replaced by a convergence of curves for different $L$’s. We have also verified several other twisting angles and similar results emerge. Thus, the level statistics indicates the existence of mobility edge; only when $V_{d}\gtrapprox 6.5$ would almost all eigenstates be localized as $\langle r\rangle$ approaches the Poisson limit.

The average value $\langle r\rangle$ gives an averaged behavior for all eigenstates. We next look at the eigenstate-resolved inverse participation ratio for more details regarding the mobility edge.

III.3 Inverse Participation Ratio and Mobility Edge

The inverse participation ratio (IPR) is defined as

[TABLE]

for a normalized eigenstate $m$: $\sum_{x,y}|\phi_{x,y}^{(m)}|^{2}=1$. The IPR crosses from 0 in the strongly delocalized (i.e. take $|\phi_{\boldsymbol{r}}|=1/L$) regime to 1 in the localized situation (i.e. take $|\phi_{\boldsymbol{r}}|=\delta_{\boldsymbol{r},\boldsymbol{r}^{\prime}}$). See results in Fig. 4. As mentioned previously, we average over various $\varphi_{1},\varphi_{2}$’s in Eq. (7) so as to reduce finite size effects concerning microscopic details.

From the results, there does seem to be mobility edges as we see that the transition of IPR for different eigenstates appears at different critical $V_{d}^{(c)}$’s. Since the localization lengths do not need to be uniform, it is worthwhile to check whether some eigenstates have much longer localization lengths than others such that they appear like an extended one for a finite-size sample. For this purpose, we compute systems of different sizes for the same parameter in Fig. 4. Should such a scenario occur, those localized states with larger localization legnths would take larger IPR values as one increases the system size (brighter colors should occupy larger areas). But there appears no notable finite-size effects as Fig. 4 (a) and (b) for sizes $L=40,80$ look almost identical. Thus, the IPR results indeed suggest the existence of mobility edges for the two-dimensional generalization of Aubry-Andre model.

Further, we note that the IPR feature here is quite different from what is typically seen in its one-dimensional counterparts. Here, the transition of IPR for different eigenstates changes continuously, constrained only by the particle-hole symmetry (after $\varphi_{1,2}$ averages, or in $L\rightarrow\infty$ limits) around $E_{m_{0}}=0$ for the middle state $m_{0}=L^{2}/2$. That is, there appear to be no “plateaus” of constant $V_{d}^{(c)}$ for certain energy windows. This is to be compared with the one-dimensional scenario where only a few $V_{d}^{(c)}$ exits, related by a “step”-like transition at certain $E_{m}$’s (see, e.g. Ref Li et al. (2017)).

We trace such a difference back to the lack of gaps in the density of states for our two-dimensional models. In the one-dimensional system with mobility edge, large gaps (comparable with band width) typically exist in the density of states. Usually, an abrupt change of $V_{d}^{(c)}$ occurs when the eigenstates go from one gapped band to another, while $V_{d}^{(c)}$ remains almost constant within one band Li et al. (2017). In our two-dimensional case, we note that for small size systems, such as the $\theta=\pi/5$ example above for $L\leq 20$, gaps in the density of states do appear, together with “plateaus” of $V_{d}^{(c)}$ for $E_{m}$ within one band. However, as the system size $L$ increases, those gaps in the density of states shrink, and eventually vanish together with the “plateaus” of $V_{d}^{(c)}$. Then, a smooth change of IPR affected by density of states at $E_{m}$ takes over, with typical features shown in Fig. 4. After this point, the IPR configurations no longer change notably with the increase of $L$. Such a scenario shows up for all cases we have tried, including various rotation angles $\theta$ and stretched lattice constants for perturbing potentials. Thus, we expect that for a generic two-dimensional, quasi-periodic system in the thermodynamic limit, no gapped structure with “plateaus” of $V_{d}^{(c)}$ should appear. The lack of gaps in the density of states is also shown for three-dimensional quasi-periodic models in Ref. Devakul and Huse (2017). It is interesting to note that such a “mini-band” type of finite-size effects also exists in certain strongly interacting systems Papić et al. (2015). As such, for all results discussed in this work, we have made sure that the system size is well above the limit for band gaps to appear.

III.4 Critical Exponent

Near the critical point $V_{d}^{(c)}$, the localization length is expected to diverge

[TABLE]

with the scaling exponent $\nu$. For the two-dimensional quasi-periodic model in the orthogonal class, such an exponent has been evasive in the previous studies Devakul and Huse (2017); Sutradhar et al. . Here, by sampling a large number of systems with different $\varphi_{1,2}$’s based on the “multifractal” analysis, we obtain such an exponent as shown in Fig. 5. One advancement of our work is that the sizes $L$ considered are systematically larger than those in previous works. For instance, the system sizes in Ref. Devakul and Huse (2017) ranges over $L=40\sim 120$, while Ref. Sutradhar et al. uses $L=10\sim 90$ in the exact diagonalizations for two-dimensional systems. We noticed that in our example of $\theta=\pi/5$, for size $L\sim 20,40,60$, there are severe finite-size effects as $\tilde{\alpha}_{0}$ defined below in Eq. (14) shows irregular crossovers for different system sizes. As we tested, such finite size effects are much more severe than the one-dimensional Aubry-Andre model or a factorizable two-dimensional one ($\theta=0$ but with incommensurate stretching) where systems with similar sizes $L$ already provide good scaling behaviors. In the following, let us describe the procedure of extracting $\nu$ in more detail.

The method of multifractal finite size scaling (MFSS) was introduced in Refs. Rodriguez et al. (2009, 2010, 2011). The steps can be summarized as follows. (1) Partition the system of size $L\times L$ into boxes of size $l\times l$, each containing $(L/l)^{2}$ lattice sites. In our case, $L/l=10$, so there are totally 100 boxes. (2) Obtain the real-space eigenstate $\phi_{x,y}$ with energy close to certain value, i.e. $E=0$ in our case. (3) Compute the box-averaged value $\tilde{\alpha}_{0}$, where we have avoided boxes on the boundary so as to reduce finite-size effects,

[TABLE]

Here, $A_{a,b}$ is the wave function amplitudes within the box indexed by $(a,b)$. For each system size $L$, we need to diagonalize a large number of Hamiltonians with different $\varphi_{1},\varphi_{2}$. Each of such a diagonalization (using Lanczos method) would render one eigenstate $\phi_{x,y}$ closest to $E=0$ and its corresponding $\tilde{\alpha}_{0}$. Those $\tilde{\alpha}_{0}$’s for the same system sizes and $V_{d}$’s are to be averaged over to provide one data point in Fig. 5. As emphasized in Ref. Rodriguez et al. (2011), only one sample state $\phi_{x,y}$ should be used in each diagonalization.

The qualitative behavior of $\tilde{\alpha}_{0}$ can be expected as follows. When the system is completely delocalized, i.e. wave function amplitudes are homogeneously distributed among all sites/boxes, $|\phi_{x,y}|^{2}\sim 1/L^{d}$, $A_{a,b}\sim(l/L)^{d}$, so $\tilde{\alpha}_{0}\rightarrow d=2$ where $d$ is the spatial dimension. In contrast, when approaching strongly localized regime with most of the boxes empty $A_{a,b}\rightarrow 0$, $\tilde{\alpha}_{0}\rightarrow\infty$. Thus, near the critical regime when the system crosses from metallic to insulating behavior at certain energy $E$, $\tilde{\alpha}_{0}$ is a monotonically increasing function of $V_{d}$.

The single-parameter scaling form can be written as

[TABLE]

which suggests the data collapse of $\tilde{\alpha}_{0}$ in terms of $L/\xi\propto|V_{d}-V_{d}^{(c)}|^{\nu}L$. Near the critical point, we expand $\tilde{\alpha}_{0}\approx\sum_{n=1}^{N}a_{n}w^{n}$, where $w=(V_{d}-V_{d}^{(c)})L^{1/\nu}$, and fit the data for different $V_{d}$ and $L$ with the polynomial. The quality of fit is evaluated by the $\chi^{2}$ statistics and the $p$-value, and the uncertainty range is generated by a Monte-Carlo type of procedure using synthetic data sets having the same mean and standard deviation as the original data. Detailed fitting procedures are described in Appendix E.1.

The results are shown in Fig. 5, where we see that the scaling exponent $\nu\approx 1$ appears to saturate the Harris bound $\nu>2/d=1$ in two dimensions. To verify the generality of the scaling exponents, we computed another twist angle $\theta=4\pi/9$ in Appendix E.2, where $\nu\approx 1$ is similarly obtained. It is also of interest to see that despite the existence of mobility edge in the two dimensional model, its critical exponent $\nu$ is similar to the one-dimensional Aubry-Andre one with analytically obtainable $\nu=1$, where mobility edge is absent. As such, we can summarize the critical behavior for the self-dual models in various dimensions in Table 1.

IV Experimental Signatures

In principle, one could surely produce exactly the model in Eq. (3) by superposing a standard square optical lattice together with the Moire perturbing potential given by a digital mirror device, like in the experiment Ref. Choi et al. (2016). Even without using the digital mirror device, a standard superlattice technique could realize a generalized version of the previous model. For the generality and broadness of our results, we would focus on such a generalized model in the following.

In the superlattice scheme, we need two identical and mutually independent lattice potentials (without interference) overlapping with each other. Note that two square lattices involve 4 pairs of laser beams of the same wavelengths. Since there are only 3 perpendicular directions for polarizations, at least two pairs of the beams would end up interfering with each other and therefore the two sets of square lattice will not be mutually independent. To circumvent this difficulty, we consider a generalization to two identical rectangular lattices, where laser beams of very different frequencies will not interfere with each other when their effects are averaged over time. That means the Moire potential becomes

[TABLE]

with $x,y\in\mathbb{Z}$ still denoting site indices. The self-duality mapping can be similarly performed by replacing $u_{x,y},v_{x,y}\rightarrow\tilde{u}_{x,y},\tilde{v}_{x,y}/\gamma$ in Eq. (9), and therefore the self-dual point is still $V_{d}^{(D)}=4$. In the following, we choose the aspect ratio of $\gamma=2$ for a plaquette of the rectangular lattices. This is to emphasize that there is no incommensurability arising from the lattice constants of two overlapping lattices, and also, such an aspect ratio can be achieved by using lasers with wavelengths $\lambda_{1}=1064nm$ and $\lambda_{2}=532nm$ as in many experiments. A schematic plot for the lasers and their polarization directions are shown in Fig. 6. Note that the commensurate condition Eq. (4) still holds for rectangular lattices with rational aspect ratios, as discussed in Appendix A.

Consider a spin-polarized fermion system. The initial state can be chosen similar to the previous experiments as a charge-density-wave type Schreiber et al. (2015), where even sites $x=0,2,4,\dots$ along a direction are empty, while odd sites $x=1,3,5,\dots$ are occupied. Fig. 7 shows the numerical simulation for the time evolution, where $N_{\text{odd}}$ are particles at sites $x=1,3,5\dots$ for arbitrary $y$, and $N$ is the total particle number. We readily see that close to self-duality point $V_{d}^{(D)}=4$, i.e. $V_{d}=3$, there already appear a significant portion of localized states such that $N_{\text{odd}}/N$ remain slightly above $0.5$, which signals the mobility edge in such a system.

V Conclusion

We have provided, as the first step, a study on effective models of Moire localization in two-dimensional quasi-periodic systems. The commensuration conditions discussed in Sec. II can be used as a guide to engineer bilayer or multi-layer systems made of simple rectangular/square lattices, which is of primary interest to cold atom experiments.

Such a model also fills up the gap for the understanding of two-dimensional quasi-periodic system in the orthogonal class (time-reversal invariant). We find that the mobility edge here, shown by the IPR, appears quite differently from that in a one-dimensional quasi-periodic system studied before. In our case, there exists a rapid and continuous change without “plateaus” for IPR in the eigenstate-disorder plane. We trace such a difference to the lack of band gaps for the two-dimensional spectrum. Note that such a feature is most notably revealed in the IPR calculation rather than the level statistics, as the latter one would inevitably involve an average over many nearby eigenstates. On the other hand, level statistics reveals a cross-point for relatively small system sizes but appear to converge for larger $L$’s, with the region well above duality point $V_{d}=4$ still far away from the Poisson limit. Such a fact also confirms the existence of mobility edge.

We have also found the scaling exponent for such a model, which has been absent in previous studies Devakul and Huse (2017). It also resolves the earlier concern in Ref. Sutradhar et al. that a single-parameter scaling may be insufficient for this class of model, as we found that choosing sufficiently large system sizes would yield satisfactory scaling behaviors. The exponent for our two-dimensional model, $\nu\approx 1$, turns out to be the same as its one dimensional cousin, while the difference is that in two-dimensions such an exponent saturates the Harris bound $\nu>2/d=1$ formulated for a purely random scenario. Also, the exponent in our orthogonal class is smaller than the two-dimensional systems in symplectic classes $\nu\approx 1.38$ Devakul and Huse (2017).

Finally, we have provided a scheme for cold atom experiments relying solely on standard superlattice techniques, without using the digital mirror device. Since the system is only sensitive to the relative rotation angles rather than lattice constants between main and perturbing potentials, our results indicate that the two dimensional system is actually more viable for cold atom experiments as there are no requirements for particular laser wavelengths.

VI Acknowledgment

This work is supported by AFOSR Grant No. FA9550-16-1-0006, MURI-ARO Grant No. W911NF-17-1-0323, and NSF of China Overseas Scholar Collaborative Program Grant No. 11429402 sponsored by Peking University. This work used the Bridges system, which is supported by NSF award number PHY190004P at the Pittsburgh Supercomputing Center (PSC).

Appendix A Commensurate rotations for rectangular Lattices

Take an arbitrary lattice site as the center of rotation, then the commensurate rotations should satisfy

[TABLE]

where $\hat{x}^{\prime}=\hat{x}\cos\theta+\hat{y}\sin\theta$, $\hat{y}^{\prime}=-\hat{x}\sin\theta+\hat{y}\cos\theta$, and $\gamma$ is the aspect ratio for the rectangular lattices. The following derivations are valid for a rational $\gamma$, i.e. $\gamma=2$ for the experiment part in the main text. Written explicitly,

[TABLE]

To have integer solutions $(m_{1},m_{2}),(n_{1},n_{2})$, it is necessary to have

[TABLE]

The integers $k_{1},k_{2},k_{3}$ satisfy

[TABLE]

The above diophantine equation corresponds to finding Pontryagin triples whose solutions are well-known Dickson (2005); Steuding (2005); Shallcross et al. (2010). One can enumerate them by the following procedure. Note that Eqs.  (19) and (20) mean finding the rational points on the unit circle $\cos^{2}\theta+\sin^{2}\theta=1$. One such point can be pinned down as $(0,1)$, and the remaining ones can be found by intersecting the circle with the line passing through rational points $(0,1)$ and $(q/p,0)$, where $q,p\in\mathbb{Z}$. Combining (20) and $k_{1}/k_{3}=1-(q/p)k_{2}/k_{3}$, we can set

[TABLE]

Then we have Eq. (4) in the main text. Note here $k_{1},k_{2}$ are interchangeable, which correspond to angles $\theta$ and $\pi/2-\theta$.

Note that the commensuration condition is the same for arbitrary aspect ratios $\gamma$. But the Bravais lattice vectors are different, as we discuss below.

Appendix B Moire Bravais lattice vectors at commensurate rotation angles

Define the matrices

[TABLE]

where T is transpose, and $A$ is the rotation matrix in Eq. (18) parametrized by Eq. (19) and (21). They satisfy the relation $SAS=\mathbb{I}$. Then the vectors $\boldsymbol{m},\boldsymbol{n}$ can be written as

[TABLE]

where $\alpha,\beta$ are some rational numbers. Thus, the Moire Bravais vectors for a square lattice $\gamma=1$ can be chosen as

[TABLE]

We aim to find the shortest vectors $\boldsymbol{T}_{1,2}$ as it is possible that some linear combinations could yield

[TABLE]

where $N_{i}\in\mathbb{Z}$ and $\boldsymbol{T}^{\prime}_{i}$ is also a Moire Bravais vector. Note here $(\alpha_{i},\beta_{i})$ are coprime numbers, for otherwise the common factors can be canceled by $N_{i}$. It turns out that the only exception is when $(q,p)$ are both odd numbers, in which case

[TABLE]

For a rectangular lattice with rational aspect ratio $\gamma$, one should apply the general conditions (24) to find the linear independent vectors $\boldsymbol{T}=(n_{1},n_{2})$ such that $\alpha q+\beta p$ and $(-\alpha p+\beta q)/\gamma$ are integers, and $q,p$ are coprime. For our purpose, we give the lattice vectors for $\gamma=2$: when $q,p$ are both odd, $\boldsymbol{T}_{1}=(q+p,(q-p)/2)^{T},\boldsymbol{T}_{2}=(-(q-p),(q+p)/2)^{T}$; when $q$ is odd and $p$ is even, $\boldsymbol{T}_{1}=(q,-p/2)^{T},\boldsymbol{T}_{2}=(2p,q)^{T}$; finally, when $q$ is even and $p$ is odd, $\boldsymbol{T}_{1}=(2q,-p)^{T},\boldsymbol{T}_{2}=(p,q/2)$.

Appendix C Numerical algorithm for free fermion dynamics

Consider

[TABLE]

with $c_{i}$ the fermion operator at site $i$ and ${\cal H}$ an $N\times N$ matrix. For an observable that can also be expressed in the bilinear form (such as the density)

[TABLE]

note $[AB,C]=A\{B,C\}-\{A,C\}B$,

[TABLE]

and therefore

[TABLE]

Now, if the Hamiltonian matrix can be diagonalized by the matrix ${\cal T}$,

[TABLE]

we have

[TABLE]

For our purposes, consider the initial state where half of the system is filled, $|\psi_{\text{ini}}\rangle=c_{1}^{\dagger}\dots c_{N/2}^{\dagger}|0\rangle$. Then

[TABLE]

where $\alpha,\beta$ are summed over all eigenstates, and $n_{i}^{(0)}$ is the initial particle number at site $i$. The above formula applies to both free bosons and fermions when the initial state is a Fock one, with the restrictions that for fermions, $n_{i}^{(0)}\leq 1$.

Appendix D Review of Aubry-André model in one-dimension

To be self-contained and to help trace the origin for mobility edge in two-dimensional models, we review some features of the one dimensional Aubry-Andre model. There are two necessary conditions for having uniform localization lengths: (1) A one-dimensional model with nearest-neighbor hopping; (2) The duality between localization/delocalization in real/momentum spaces. Generalizing to higher-dimensions clearly violate condition (1) and therefore do not exhibit unique localization length. The following derivation is based on the discussions in Refs. Herbert and Jones (1971); Thouless (1972); Aubry and André (1980); Goldsheid and Khoruzhenko (2005).

The Aubry-Andre model is

[TABLE]

Performing a similar analysis as in Sec. III.1, we could similarly obtain the duality mapping with critical disorder strength $V_{d}^{(D)}$. However, due to the one-dimensional nearest neighbor feature, one could related the localization length with density of states through the Thouless formula as shown below.

Consider the real space matrix for Green’s function

[TABLE]

where $m,n=1,\dots,N$ denote lattice sites, $|0\rangle$ is the vacuum. For an open boundary system, the matrix

[TABLE]

Due to the one-dimensional nearest-neighbor feature, the corner element for the inverse matrix can be easily acquired as the cofactor

[TABLE]

where $E_{\alpha}$ are eigenvalues of $H$. On the other hand, according to the definition, one can first expand the Green’s function in the energy eigenbasis $G(E)=(E-H)^{-1}=\sum_{\alpha=1}^{N}(E-E_{\alpha})^{-1}|\alpha\rangle\langle\alpha|$, and then the real space matrix element is

[TABLE]

where $\psi_{\alpha}(m)=\langle 0|a_{m}|E_{\alpha}\rangle$ is the real-space eigenfunctions. Simultaneously using the above two equations to compute the residual of $G_{1N}(E)$ at a certain $E_{\beta}$,

[TABLE]

Suppose the eigenstate $|E_{\beta}\rangle$ is localized centering at some site $m_{0}$, with the ansats

[TABLE]

where $\lambda$ is the exponent of interest, and $A_{m}$ is the multiplication of some constant and non-decaying oscillating factors at site $m$. Combing the above two equations and taking logarithm, we have

[TABLE]

Taking the $N\rightarrow\infty$ limit, the $\ln(A_{1}A_{N})$ term vanishes. The second term can be converted into an integration. Then

[TABLE]

where $dN(\varepsilon)=D(\varepsilon)d\varepsilon$ and $D(\varepsilon)$ is the density of state at $\varepsilon$. This is the Thouless formula for the localization length $\xi=1/\lambda$. Usually, the hopping constant $t=1$, so the first term drops out. This formula can be easily generalized to inhomogeneous nearest neighbor hopping situations Thouless (1972). But note that once the hopping goes beyond nearest neighbor, i.e. for second-neighbor hopping or in higher dimensions, Eq. (39) no longer takes the same form, and Eq. (40) cannot be obtained. Then the Thouless formula would not hold generally.

Due to the duality mapping, one can obtain the relation

[TABLE]

for an incommensurate system. Then, applying Thouless formula,

[TABLE]

Since $\lambda\geq 0$ by definition, for $V_{d}>4t$, $\lambda_{t,V_{d}}(E_{\beta})\geq\ln(V_{d}/4t)>0$. That means all eigenstates $|E_{\beta}\rangle$ are exponentially localized. From the transformation between eigenfunctions in two dual representations, all eigenfunctions in the dual space are then delocalized (localization/delocalization are opposite in the two spaces related by Fourier transformation), meaning $\lambda_{t,8t^{2}/V_{d}}(4tE_{\beta}/V_{d})=0$. Then, we have the same critical exponent (inverse of localization length $\xi$)

[TABLE]

for all eigenstates when $V_{d}/4t>1$.

Appendix E Additional details for scaling exponents

E.1 Fitting procedure

The procedure for fitting the scaling of $\tilde{\alpha}_{0}(V,L)$ consists of three steps.

First, the average $\langle\tilde{\alpha}_{0}\rangle(V,L)$ and standard deviation $\sigma(V,L)$ of $\tilde{\alpha}_{0}$ at each point $(V,L)$ is determined. There are data of the order $10^{5}$ at each $(V,L)$ calculated by using different $\varphi_{1},\varphi_{2}$. Note that although the average value is straightforward to understand, the “standard deviation” must be taken with care. In all regimes, the $\tilde{\alpha}_{0}$ takes values according to a broad distribution function (see, i.e. Ref. Rodriguez et al. (2011) Fig. 3). We reproduced such a character for our model in Fig. 8. It is in fact the maximal point of the distribution function that gives the scaling behavior. Only when averaged over large numbers of $\varphi_{1,2}$, the mean value of $\tilde{\alpha}_{0}$ would approach the maximal probability point. Thus, at each $(V,L)$, we divide all data into 10 bins, obtain the respective averaged $\tilde{\alpha}_{0}$, and compute the standard deviation using the 10 bin-averaged $\tilde{\alpha}_{0}$’s.

Second, we fit $\langle\tilde{\alpha}_{0}\rangle(V,L)$ into the scaling function $g(L/\xi)$ mentioned in the main text by minimizing $\chi^{2}$, $\chi^{2}=\sum_{V,L}(\langle\tilde{\alpha}_{0}\rangle(V,L)-g(V,L))^{2}/\sigma^{2}(V,L)$. For a reasonable fit, the $\chi^{2}$ value should be similar or less than the number of degrees of freedom, $k=N_{D}-N_{p}$, where $N_{D}$ is the number of $\langle\tilde{\alpha}_{0}\rangle(V,L)$ to be fit for different $(V,L)$, and $N_{p}$ is the number of fitting parameters. Whether or not a fit is acceptable is decided by the $p$-value,

[TABLE]

where $\Gamma(a,z)=\int_{a}^{\infty}t^{z-1}e^{-t}dt$ is the generalized Euler-Gamma function, with $\Gamma(z)=\Gamma(0,z)$. A larger $p$-value signals a better fit, and we take $p\geq 0.1$ as the acceptable criterion. In addition to the $p$-value, we also make sure the fitting stability by checking that increasing the expansion order would result in $V_{d}^{(c)},\nu$ within the uncertainty range (to be discussed below). Under such criterions, the expansion order is kept as low as possible.

Finally, once a best fit $g(V,L)$ is obtained for $\langle\tilde{\alpha}_{0}\rangle(V,L)$ and $\sigma(V,L)$, we generate 1000 synthetic data $\alpha^{(syn)}_{i}(V,L),i=1,\dots,1000$, at each point $(V,L)$. They distribute according to the same mean and standard deviations discussed previously. Then, we fit each set of $\{\alpha^{(syn)}_{i}(V,L)|\forall V,L\}$ with the same expansion as in the best fit $g(V,L)$, and obtain a histogram of $V_{d}^{(c)},\nu$ corresponding to $i=1,\dots,1000$. After discarding the ending (maximal and minimal) $2.5\%$ of data, the uncertainty range is obtained. The histogram of $V_{d}^{(c)},\nu$ for $\theta=\pi/5$ is shown in Fig. 9.

E.2 Scaling at another angle

We perform the same analysis for the angles $\theta=4\pi/9$ at $E=0$. But here we used a wider range of system lengths $L$, and it is clear from Fig. 10 that there is a shift in the “crossing” point with the increase of $L$. Such a “shift” diminishes as $L$ becomes larger, which means there is an additional scaling-irrelevant term entering $\tilde{\alpha}_{0}$ Rodriguez et al. (2011). Therefore, we include such a term in the scaling function,

[TABLE]

Compared with Eq. (15), the $L^{-|y|}$ term describes the first order expansion of scaling-irrelevant term. With the fitting parameter $-|y|<0$, this term vanishes in the thermodynamic limit. Then we expand $g_{0,1}$ similarly as before, $g_{0}=\sum_{m=0}^{N_{0}}a_{m}w^{m}$, $g_{1}=\sum_{m=0}^{N_{1}}b_{m}w^{m}$, where $w=(V_{d}-V_{d}^{(c)})L^{1/\nu}$. The fitting parameters are $a_{m},b_{m},y,\nu,V_{d}^{(c)}$, with the total number $N_{p}=N_{0}+N_{1}+5$.