Energy Spectrum Theory of Incommensurate Systems

Due to the lack of the translational symmetry, calculating the energy spectrum of an incommensurate system has always been a theoretical challenge. Here, we propose a natural approach to generalize the energy band theory to the incommensurate systems without reliance on the commensurate approximation, thus providing a comprehensive energy spectrum theory of the incommensurate systems. Except for a truncation dependent weighting factor, the formulae of this theory are formally almost identical to that of the Bloch electrons, making it particularly suitable for complex incommensurate structures. To illustrate the application of this theory, we give three typical examples: one-dimensional bichromatic and trichromatic incommensurate potential model, as well as a moir\'{e} quasicrystal. Our theory establishes a fundamental framework for understanding the incommensurate systems.

In this work, we propose a natural approach to generalize the energy band theory of the Bloch electrons to the incommensurate systems, by which the energy spec-trum of the incommensurate systems can be conveniently calculated without the need for any commensurate approximation.The formulae of such incommensurate energy spectrum (IES) theory are formally almost identical to those of Bloch electrons, with the only difference being a truncation-dependent weighting factor for each eigenstate.Therefore, it provides a unified method for handling both commensurate and incommensurate potentials, making it an ideal choice for multiple periodic potential models.The IES theory establishes a fundamental theoretical framework for comprehending incommensurate systems.
The Hamiltonian of the BIP model is where V 1,2 are the amplitude of the two periodic potentials.G 1,2 are the magnitude of the reciprocal lattice vectors.We define α = G 2 /G 1 as the ratio between the two periods, which is an irrational number in the incommensurate case.
Using plane wave basis, we get the Schrödinger equation in momentum space Here, c q is the coefficient of the plane wave with φ(x) = q c q e iqx , and ε is the eigenvalue to be determined.The Eq. ( 2) represents a set of algebra equations, which couple only the plane wave states in the set

arXiv:2309.01367v1 [cond-mat.mes-hall] 4 Sep 2023
To obtain the energy spectrum, some kind of truncation of (m, n) should be given first, so that the Hamiltonian H(q) becomes a matrix of finite dimensions.Then, for a given q, we can directly calculate the eigenstates of the Hamiltonian matrix H(q).Now, the key question becomes how to determine range of values for q, so that the entire energy spectra of the incommensurate system can be obtained without any omissions or duplications.
To solve this question, let us first review the case of Bloch electrons, i.e.V 2 = 0.In this case, we can also get a set of coupled equations as Eq. ( 2), which are exactly the so-called central equations of the Bloch electrons [96].Here, the coupled plane waves are Q B q = {k|k = q + mG 1 : m ∈ Z} now.In principle, q ∈ (−∞, +∞).But, to calculate the eigenstates of H(q), some momentum q are equivalent or duplicated.It is because that the two momenta q and q + mG 1 (m ∈ Z) actually correspond to the same set of central equations.In other words, H(q) and H(q + mG 1 ) are the same matrix, thereby yielding identical energy spectra.In this sense, we say that all the momenta in Q B q are equivalent.From such equivalence relations, we can get two important facts: 1.For any arbitrary momentum q, there exists an integer m 0 that yields an equivalent momentum q 0 = q + m 0 G 1 within the range [0, G 1 ).This implies that the complete energy spectra can be obtained by considering only momenta within [0, G 1 ).
2. Any two momenta in the interval [0, G 1 ) are unequivalent.The implication is that the energy spectra are not computed repetitively while q traverses the interval [0, G 1 ).
Thus, the conclusion is that, by traversing the interval [0, G 1 ) with q, we can calculate the eigenstates of all H(q) and obtain the complete energy spectra of the periodic system without repetition.Obviously, this is exactly the standard procedure of the energy band theory of the Bloch electrons, and the interval [0, G 1 ) is just the first Brillouin Zone (FBZ).
The similar idea can be generalized to the incommensurate case.As indicated in Eq. ( 2) and ( 3), all the momenta in Q q = {k|k = q + mG 1 + nG 2 : m, n ∈ Z} are equivalent.Such equivalence relation can also give us several important messages: 1.With a specified n, we can always find an equivalent momentum q n = (q+nG 2 )+m n G 1 within the range of [0, G 1 ) for any given q, where m n is an integer.
2. Suppose a truncation of n is provided, and N E represents the total number of all the allowed n.
Then, for any q, there are N E equivalent momenta in the interval [0, G 1 ), which do not overlap with each other due to the incommensurability.denotes value of n (k), so that all equivalent momenta k = q + mG 1 + nG 2 for a given q can be represented as discrete dots on the graph.Moreover, for a specified n, the equivalent momenta (q + nG 2 ) + mG 1 with different values of m lie on a horizontal line.Within each horizontal line, we see that there always exists one equivalent momentum within the range [0, G 1 ) (red dot), which is just the q n mentioned above.A natural truncation scheme of (m, n) is: Here, n c and k c are two truncation constants, where k c reflects the truncated energy of the plane waves and n c determines the minimum interval between all the coupled momenta [97].The truncation of n (k) is represented as the two black horizontal (orange vertical) dashed lines in Fig. 1 (a).Once n c is given, the total number of allowed n (red dots) are N E = 2n c + 1.It can be proved that all the red dots (q n ) never overlap with each other in the incommensurate case [98].Therefore, for any q, we can always find N E distinct equivalent momenta in the interval [0, G 1 ).
The two points above suggest that a correct and convenient way to calculate the energy spectra of an incommensurate system is to let q traverse the interval [0, G 1 ), calculate the eigenstates of each H(q) and assign a weight factor of 1/N E to each eigenstate.
Based on this idea, the density of states (DOS) and expectation value of an observable in an incommensurate system can be calculated by the following formulae: where ε qi and |ϕ qi ⟩ denote the ith eigenvalue and eigenstate of H(q), and p qi is the Boltzmann factor.Here, we name the interval [0, G 1 ) primary Brillouin zone (PBZ) of the incommensurate systems.Note that we can also choose G 2 as the PBZ, which will give the same results with proper truncation [99].The two formulae above are formally almost identical to that of Bloch electrons, which are the central results of the proposed IES theory.
The IES method can provide accurate enough numerical results as long as the cutoff k c and n c is large enough, since the plane wave basis is used.In practical calculations, we need to test different truncation values to ensure the convergence of results.An interesting case is n c → ∞, where the equivalent momenta q n of an incommensurate model will uniformly and densely fill the whole PBZ.It implies that, once n c is sufficiently large, the entire incommensurate spectra can be obtained by just calculating the eigenstates of H(q) with one any chosen q [75,95].But the cost is the dimension of the Hamiltonian matrix is extremely huge.
The IES method, i.e.Eq.( 4) and ( 5), is also valid for the commensurate case.The only difference is that the N E for the commensurate case has to be determined by directly counting the distinct q n in PBZ, since the the equivalent momenta in the PBZ (q n ) now may overlap with one another [100].In this sense, the incommensurate and commensurate systems exhibit no fundamental distinction, which thus can be treated under a unified theoretical framework.
Results of the BIP model.-Wethen apply this IES method to three specific examples for demonstration.The first one is the BIP model, where we set α = ( √ 5 − 1)/2.It is a typical incommensurate case and has been intensively studied with various methods [26,[28][29][30][31].
The numerical results are given in Fig. 1.First of all, like the Bloch electrons, we can plot energy spectra ε qi as a function of q [see Fig. 1 (b)], where i become the "band" index and q ∈ PBZ.Such energy spectrum diagram clearly indicates the existence of the energy gaps.However, the incommensurate energy spectrum diagram here is different from the energy bands of Bloch electrons, because its shape depends on the N E , i.e. the truncation n c .Note that, in the PBZ, there are N E identical spectra at the N E equivalent momenta, which is quite unlike that of Bloch electrons.The other obvious characteristic is that there exist a kind of momentum edge states (orange lines) in the energy gaps.Formally, Eq. ( 2) can be viewed as a tight binding (TB) model in momentum space, where a plane wave is equivalent to a discrete site in momentum space.Thus, when the truncation n c is given, it actually results in open boundaries in momentum space.Like the TB model in real space, the open boundary will give rise to edge states, i.e. the so-called momentum edge states in Fig. 1 (b).We can distinguish the momentum edge states by examining the distribution of the wave functions in momentum space [102].Since the momentum edge states rely on the truncation n c , we think that they are artificially induced states that need to be eliminated.
To demonstrate the correctness of the IES method, we compare it with the commensurate approximation.The DOS calculated with the two methods are plotted in Fig. 1 (c), where the blue (red) lines represent the IES method (commensurate approximation with α ≈ 55/89).Both methods yield the same DOS, which clearly indicates that the IES method can correctly describe the en-ergy spectrum of the incommensurate system.
The IES method can accurately capture the wave function characteristics of the incommensurate system as well.A direct example is the single particle mobility edge of the BIP model, which refers to the quantum localization of the wave functions in an incommensurate system.With the IES method, the inverse participation ratio in momentum space (IPRM) can be utilized to quantify the degree of localization [103], If the wave function is localized in real space, it should be extended in momentum space, and thus the IPRM will approach to zero, i.e.Clearly, when V 2 is small (large), the wave function is extended (localized).Moreover, we can see that the localization transition is energy dependent (the black oblique line), which indicates that the single particle mobility edge indeed exists.The results from the IES method are completely consistent with the known conclusions [29,30].Very interestingly, the IES method offers an improved approach to calculate the famous Hofstadter butterfly spectrum.It is known that in the TB limit (with deep V 1 potential), the BIP model will asymptotically approach the AAH model, which can be exactly mapped into a 2D Hofstadter model, i.e. 2D square lattice in a perpendicular magnetic field [104].Thus, when the energy spectrum of the BIP model is plotted as a function of α, we indeed get a butterfly-like spectrum, as shown in Fig. 1 (e).Since α is proportional to the magnetic flux ϕ through each plaquette, it means that the IES theory provides a way to calculate the Hofstadter butterfly spectrum under any magnetic field, no matter whether ϕ is rational or irrational.
Results of the TIP model.-TheIES method offers a convenient approach for calculating complex incommensurate systems.One example is the 1D trichromatic incommensurate potential (TIP) model, where three applied incommensurate potentials greatly hinder the commensurate approximation.The Hamiltonian is now where α 2 = G 2 /G 1 and α 3 = G 3 /G 1 are two irrational numbers [101].The coupled momenta are now  can be determined by directly enumerating the distinct equivalent momenta in PBZ.Then, energy spectrum and other properties can be calculated with Eq. ( 4) and (5) in the same way.
The numerical results are given in Fig. 2. Fig. 2 (a) shows the energy spectrum.With three incommensurate potentials, the number of equivalent momenta N E is markedly greater than that in the BIP model, resulting in a highly dense energy spectrum.The presence of momentum edge states is also observed in this case (orange lines).We have calculated the corresponding DOS (blue lines) in Fig. 2 (b), which coincides well with that of an approximate commensurate structure (red lines).In Fig. 2 (d), we present the IPRM to show the localization feature of the eigenstates.We also can get a butterfly-like spectrum by plotting the energy spectra as a function of α 3 with a fixed α 2 , see Fig. 2 (e).
Results of moiré quasicrystal.-Thetwo dimensional (2D) incommensurate systems have drawn great interest very recently.Here, we use the IES method to calculate a special 2D incommensurate system, i.e. moiré quasicrystal, which has eightfold rotational symmetry but no translation symmetry.Very similar to the BIP model, such moiré quasicrystal can be achieved by applying two square periodic potentials V (θ 1 = 0) and V (θ 2 = π/4) with a twist angle θ = π/4, where G a and G b are the two reciprocal lattice vectors of a square lattice, and R(θ) is the rotation matrix.The moiré quasicrystal potential is plotted in Fig. 3 (a), revealing its evident eightfold rotational symmetry.Note that such moiré quasicrystal has already been realized in cold atom systems [23][105].For a given q, the coupled momenta are now The equivalent momenta and the truncation are illustrated in Fig. 3 (b).Here, we choose the FBZ of V (θ 1 = 0) as the PBZ, see the gray square.All the equivalent momenta are plotted as dots in Fig. 3 (b), with a truncation: |n 1,2 | ≤ n c , |k| ≤ k c .k c denotes the energy cutoff of the plane waves, illustrated as the orange circle.For the moiré quasicrystal, all the equivalent momenta in PBZ will never overlap with one another (red dots) [106], which implies N E = (2n c + 1) 2 .The calculated DOS is given in Fig. 3 (c).In Fig. 3 (d), we plot the IPRM of the moiré quasicrystal.The IES theory does provide a convenient method for handling the moiré quasicrystal.All the calculation details of the three examples are given in the supplementary materials [107].
Summary.-In summary, we establish a general energy spectrum theory for the incommensurate systems without relying on the commensurate approximation.In addition to the examples above, it can be readily applied to the TB models and other artificial incommensurate systems, e.g.cold atom systems and optical crystals.More importantly, it implies that all the formulae of the energy band theory can be directly transplanted into the incommensurate systems via the IES theory.Thus, our theory actually gives a fundamental framework for comprehending the incommensurate systems.

III. THE PHYSICAL SIGNIFICANCE OF THE TRUNCATION
To calculate the incommensurate energy spectrum, we need to introduce a truncation of (m, n) to get a finitedimensional Hamiltonian matrix H(q).As discussed in the BIP model, the used truncation is |n| < n c and |k| < k c , where n c and k c are two provided truncation constants.The physical meaning of k c is very clear, which represents the cutoff of the energy of plane waves.Meanwhile, n c also has a clear physical significance, which indeed reflects the minimal interval between the equivalent momenta for given q.Actually, due to the incommensurability, the combination of G 1 and G 2 , i.e. mG 1 + nG 2 , can give rise to arbitrarily small momentum, which is essentially different from the commensurate case.When |n| < n c , no matter what the value of m is, |mG 1 + nG 2 | always has a minimum value, which depends on the value of n c .Such physical pictures of the truncation is valid for all the incommensurate systems.

IV. PRIMARY BRILLOUIN ZONE
In principle, it is better to use a deep potential as the PBZ, and a shallow potential as a perturbation.The advantage lies in the fact that shallow perturbation can lead to a small truncation n c , as demonstrated in the main text.The converse is also valid, but the truncation should be larger.In Fig. S1 (c) and (d), we use the interval [0, G 2 ) as the PBZ and calculate the DOS and energy spectrum.We see that, with a larger truncation n c , we get the same DOS.

V. THE BIP MODEL
Here, we give some calculation details of the BIP model.Fig. S2 plots the calculated DOS with different truncation n c and k c , which illustrates the convergence of the results.The DOS results show that we can get a convergent numerical results even with truncation n c = 4.And the results from the IES theory are in good agreement with the commensurate approximation as well.To calculate the DOS, we uniformly choose 1000 q points in the PBZ.For each q point, the H(q) is a 137 × 137 matrix with n c = 8 and k c = 4G 1 .

VI. THE TIP MODEL
Here, we give some calculation details of the TIP model.Fig. S3 (a) is the same as the Fig. 2 of the main text.Then, we plot the enlarged energy spectrum diagram (region in the red box) in Fig. S3 (b).The DOS with various truncations are shown in Fig. S3 (c) to illustrate the numerical convergence, where Fig. S3 (c) and (d) are the enlarged DOS plots in the gray regions.The DOS of a commensurate approximation (G 1 : G 2 : G 3 ≈ 65 : 49 : 37) is also given as a comparison .First, we see that a convergent DOS can be achieved even with n c = 3.Meanwhile, Fig. S3 (d) and (e) indicate that the IES theory may give better results than that of the common commensurate approximation.It is because that, with the IES theory, the predicted position of peaks always remains unchanged with increasing n c .In Fig. S3, H(q) is a 393 × 393 matrix when n c = 3.

VII. MOIRE QUASICRYSTAL
First, we would like to show that the moire quasicrystal model here is exactly the same as that in Ref. 23 of the main text.In Ref. 23, the potential is Ga and G 4 = Gb , we then obtain the moire quasicrystal potential in the main text, except for a factor 1 2 for V 0 and a constant potential shift.In Fig. S4, we plot the calculated DOS with different truncation to illustrate the convergence.The results of a commensurate approximation structure is given as well.Here, a commensurate approximation structure is described by two integers (m, n) with where θ mn is the twisted angle [24].For the moire quasicrystal, (m, n) = (12, 29) corresponds to θ mn ≈ 44.96

VIII. THE COMMENSURATE CASES
Here, we use the bichromatic potential model as an example to illustrate that the IES formulae are also valid for the commensurate case.Specifically, we consider a commensurate situation: G 1 : G 2 = 3 : 1.Based on the Eq. ( 2) of the main text, we can also plot the equivalent momenta for a given q with G 1 as the PBZ, see Fig. S5 (a).Clearly, distinct from the incommensurate case, the equivalent momenta will overlap in such commensurate case.Thus, we eventually get three distinct equivalent momenta in the PBZ, i.e.N E = 3, which is irrelevant to the n c .Note that [0, G 1 /3) is just the FBZ of such commensurate case.Therefore, Eq. ( 4) and ( 5  that the IES theory actually is applicable to the commensurate cases.So, the IES formulae in fact provides a unified theoretical framework to treat the multi periodic potential models, no matter whether it is incommensurate or not.In practical calculations, the only issue is to correctly determine N E by counting the distinct equivalent momenta in the PBZ.Note that in the commensurate cases, n c has no effect and k c is the only meaningful truncation.

Fig. 1 (
Fig.1(a) can help us intuitively understand the two points above.The vertical (horizontal) axis of Fig.1(a) denotes value of n (k), so that all equivalent momenta k = q + mG 1 + nG 2 for a given q can be represented as discrete dots on the graph.Moreover, for a specified n, the equivalent momenta (q + nG 2 ) + mG 1 with different values of m lie on a horizontal line.Within each horizontal line, we see that there always exists one equivalent momentum within the range [0, G 1 ) (red dot), which is just the q n mentioned above.A natural truncation scheme of (m, n) is:(1) |n| ≤ n c ; (2) |k| ≤ k c .Here, n c and k c are two truncation constants, where k c reflects the truncated energy of the plane waves and n c determines the minimum interval between all the coupled momenta [97].The truncation of n (k) is represented as the two black horizontal (orange vertical) dashed lines in Fig.1 (a).Once n c is given, the total number of allowed n (red dots) are N E = 2n c + 1.It can be proved that all the red dots (q n ) never overlap with each other in the incommensurate case[98].Therefore, for any q, we
log 10 (IPRM) → −∞, as k c increases.Conversely, for an extended wave function, log 10 (IPRM) tends to zero.Fig. 1 (e) plots the IPRM of all the eigenstates of the BIP model as a function of V 2 , where the log 10 (IPRM) is represented by the color.

n 3 ∈
Z}.With the IES theory, we can still choose [0, G 1 ) as the PBZ and use the truncation: |n 2,3 | ≤ n c and |k| ≤ k c .Similarly, we also have N E equivalent momenta in the PBZ, which
This work was supported by the National Key Research and Development Program of China (No. 2022YFA1403501), and the National Natural Science Foundation of China(Grants No. 12141401, No. 11874160).

FIG
FIG. S3.The TIP model.(a) is the energy spectrum diagram, the same as the Fig. 2 of the main text.(b) and (c) are the enlarged energy spectrum plot for the region of the red box in (a).(c) is the DOS calculated with different truncation nc, and the results of a commensurate approximation is also given.(d) and (e) are enlarged DOS plots for the gray regions in (c).Parameters: V1 = 8E0, V2 = V3 = 0.03E0, kc = 4G1.The two irrational numbers are α2 = λ −1 , α3 = λ −2 , where λ = 1.3247 • • • is the root of the equation x 3 − x − 1 = 0 [101].