Mobility Edges in one-dimensional Models with quasi-periodic disorder
aa r X i v : . [ c ond - m a t . d i s - nn ] J a n Mobility Edges in one-dimensional Models with quasi-periodicdisorder
Qiyun Tang and Yan He College of Physics, Sichuan University, Chengdu, Sichuan 610064, China
Abstract
We study the mobility edges in a variety of one-dimensional tight binding models with slowlyvarying quasi-periodic disorders. It is found that the quasi-periodic disordered models can beapproximated by an ensemble of periodic models. The mobility edges can be determined by theoverlaps of the energy bands of these periodic models. We demonstrate that this method providesan efficient way to find out the precise location of mobility edge in qusi-periodic disordered models.Based on this approximate method, we also propose an index to indicate the degree of localizationof each eigenstate. . INTRODUCTION More than fifty years ago, P. W. Anderson pointed out in a seminal work [1] that thedisordered on-site potential will make the electronic wavefunction become localized at certainsite rather than propagating around the whole system like a plane wave. Since then, theAnderson localization was extensively studied in the condensed matter physics. In three-dimensional (3D) systems, the Anderson localization will occur only when the disorderstrength is larger than certain value. When the disorder is not strong enough, the bandstructure will contain both extended or localized eigenstates. This leads to the concept ofmobility edge which is the energy separating the localized and extended energy levels[2]. Thescaling theory points out that the dimension is an important factor for Anderson localization,and the critical point for the existence of mobility edge is 2D [3, 4], that is, in 2D and 1Dsystems, even if the disorder intensity is infinitesimal small, the mobility edge cannot appear,and all wave functions are localized.For a very long time, physicist only focused on the study of systems with uncorrelateddisorders. Stating from 80’s, low-dimensional quasi-periodic systems have been widely in-vestigated. Due to the highly correlated disorders, it is still possible to have mobility edgeseven in low-dimensional systems[5–11]. A famous example is the so called Aubry-Andre(AA) model [12, 13] which is simply a 1D tight-binding model with incommensurate on-sitepotentials. One can show that the AA model has self-dual symmetry. With the increasing ofthe incommensurate on-site potentials, all the eigenstates of the AA model will undergo thetransition from extended to localized, and the self-dual point is the transition point betweenthe localized and extended states. Therefore, there is no mobility edge in the standard AAmodel, but if the model is further modified, such as the introduction of long-range hopping,then the mobility edge can appear in the improved AA model[14–20].Through the investigation of the AA model and its variants, a new class of system withslowly varying incommensurate potential has been extensively studied by Das Sarma and hisco-workers [21, 22]. They found two mobility edges and proposed a heuristic semi-analyticalmethod to solve the mobility edge in such systems. This further stimulated a lot of moretheoretical works on the mobility edges in 1D quasi-periodic systems [23–25], includinggeneralizing the disorders to the off-diagonal hopping terms [26] or to some topological1D model such as Kitaev chain [27]. Among them, the off-diagonal hopping AA model has2ttracted extensive research interests, because it contains rich and novel physical phenomenasuch as zero-energy topological edge states[28, 29].In the present paper, we will study the generalized AA model with slow varying quasi-periodic disorders in hopping and potential terms. We will first introduce the heuristicsemi-analytical method to determine the mobility edges. This method is based the propertythat the slow varying disorder approach to a constant asymptotically. Although this methodis very useful for some simple models, it also has severe limitation in finding mobility edgesfor models with more complicated disorders. To overcome this difficulty, we approximatethe quasi-disordered model with an emsemble of different periodic models. Then the regionof extended states can be approximately obtained by the overlaps of energy bands of theseperiodic models. We will demonstrate this method in details for several variations of AAmodel and also the Su-Schrieffer-Heeger (SSH) model [30]. Based on this method, we alsopropose a new index to indicate the degree of localization of eigenvectors.
II. A SEMI-ANALYTICAL METHOD TO DETERMINE THE MOBILITY EDGE
We consider the AA model with slowly varying quasi-periodic disorders. H = − N − X i =1 ( t + w i )( c † i c i +1 + c † i +1 c i ) + N X i =1 µ i c † i c i (1)Here N is the total number of lattice sites. For convenience, we assume that t = 1 as theenergy unit. The slowly varying disorder is given by w i = w cos(2 παi v ) , µ i = µ cos(2 πβi v + φ ) (2)Here w and µ is the amplitude of the disorder. α and β are certain irrational numbers,which determines the quasi-period. For example, one can assume that α = β = ( √ − / α is irrational or rational does not affect the results. This is because0 < v < i v is irrational number and make the system quasi-periodic. For later convenience, we will assume that α = β = 0 .
5. Here φ is certain possibleinitial phase angle.Since the exponent v satisfies 0 < v <
1, we have dw i di = − wvπαi v − sin(2 παi v ) , lim i →∞ dw i di = 0 (3)3herefore, because the power 0 < v <
1, the disorder term eventually become a constantfor large enough i .First, we discuss the semi-analytical method to determine the energy region of extendedstates or the locations of mobility edges. As an example, we first consider a simple casewhere β = 0, thus the potential is just a constant µ . If we assume that the eigenstate is ψ = ( a , · · · , a N ) T , then the eigenvector equation in the real space is given by( t + w j ) a j − + ( E − µ ) a j + ( t + w j ) a j +1 = 0 (4)Since w j approaches to a constant as j become very large, we find the following asymptoticeigenvector equation ( t + w ∞ ) a j − + ( E − µ ) a j + ( t + w ∞ ) a j +1 = 0 (5)Here w ∞ = lim j →∞ w j . As usual, the eigenvector is assumed to take the form a j ∝ λ j , thenwe arrive at the following characteristic equation λ + E − µt + w ∞ λ + 1 = 0 (6)which gives two roots λ , = − A ∓ √ A − , A = E − µt + w ∞ (7)In order to find extended states, λ , must be complex numbers with unit modular | λ , | = 1.Therefore, the condition for extended states is (cid:12)(cid:12)(cid:12) E − µt + w ∞ (cid:12)(cid:12)(cid:12) < w < t , the most stringent inequality one can get from Eq.(8) is obtained by taking w ∞ = − w . Therefore,we find the following condition for extended states µ − t − w ) < E < µ + 2( t − w ) (9)which agrees with the mobility edge shown in Figure 1 (a).In Figure 1, the localization of the n -th eigen-wavefunction ψ n is indicated by the inverseparticipation ratio (IPR) [31–33], which is defined asIPR n = N X j =1 (cid:12)(cid:12)(cid:12) a nj (cid:12)(cid:12)(cid:12) , ψ n = ( a n , · · · , a nN ) (10)4ere we assume that ψ n is normalized, i.e P Nj =1 (cid:12)(cid:12)(cid:12) a nj (cid:12)(cid:12)(cid:12) = 1. For the extend states, one expectsthat each component of ψ n is roughly the same order of magnitude | a nj | ∼ /N for all j .Therefore, we find that IPR ∼ P Nj =1 1 N = N , which is very small for the extended states.One the other hand, the non-zero amplitude of localized states will mostly confined to onlya few components, therefore in this case, we expect that the IPR will be order 1.In Figure 1 (a), we plot the eigenenergy of AA mdoel of Eq. (1) as a function of µ . Theparameters used in this calculation is listed in the figure caption. The color of each datapoints represents the IPR value of the corresponding eigenstates. The brighter colors indicatelarger values of IPR which represents a localized states. The darker colors indicate that IPRis close to zero which represents extended states. One can see that the dark colored areais precisely confined between the red and green solid lines which are the mobility edges wehave obtained from the semi-analytic calculations. In order to make a more quantitativelyobservation, in Figure 1 (c), we also plot the IPR as a function of eigenenrgy for a fewselected points µ = 0 , . ,
5. One can see that the IPR drops from the order of magnitudeof 10 − to 10 − when the eigenenergy across some critical values. This clearly signals thetransition from localized states to extended states. These critical values precisely matchesthe mobility edges we have discussed.We can also use the same method to consider a more complicated case like β = α , i.e.both µ and w are quasi-periodic disordered. Following similar steps, we find the conditionfor extended states is given by (cid:12)(cid:12)(cid:12) E − µ ∞ t + w ∞ (cid:12)(cid:12)(cid:12) < III. AN ALTERNATIVE METHOD TO DETERMINE THE MOBILITY EDGE
Now we turn to discuss a new method, which can determine the extended eigenstatesmore efficiently. The basic idea of this method is to approximate the above quasi-periodicdisordered model by an ensemble of periodic models which can be denoted as M a . The5 a) (b) -4 -2 0 2 4 6 8 Eigenenergy -4 -3 -2 -1 I P R (w, )=(0.5,0)(w, )=(0.5,0.5)(w, )=(0.5,5) =0.5 =0v=0.5 (c) -4 -2 0 2 4 6 Eigenenergy -4 -3 -2 -1 I P R (w, )=(0.5,0.9)(w, )=(0.5,1.6)(w, )=(0.5,2.1) =0.5 =0.5v=0.5 (d) Figure 1. (a) Eigenenergy of AA model Eq. (1) as a function of µ with α = 0 . β = 0, w = 0 . v = 0 .
5. (b) Eigenenergy of AA model Eq. (1) as a function of mu with α = 0 . β = 0 . w = 0 . v = 0 .
5. The colors in panels (a) and (b) represent the IPR of each eigenstates.The red and green curves represent the mobility edges. (c) The IPR as a function of eigenenergyfor ( w, µ ) = (0 . , , (0 . , . , (0 . , w, µ ) = (0 . , . , (0 . , . , (0 . , . Hamiltonian of M a is given by H a = − N X i =1 h ( t + w a )( c † i c i +1 + c † i +1 c i ) + µ a c † i c i i , a = 1 , · · · N (12)Here we have assumed periodic boundary condition such that c N +1 = c . For each model M a , since w a and µ a are constant, we can diagonalized the above Hamiltonian in momentum6pace and find extended Bloch wave function with the following energy band E = µ a − t + w a ) cos k (13)Therefore, the region of extended states for the model M a is just the range of the energyband µ a − t + w a ) < E < µ a + 2( t + w a ) (14)One can imagine that for the quasi-periodic disordered model, the extended states willappear for a given energy E , if this E is inside the energy bands of all the above models M a .Therefore, we expect that the energy region of extended state of Eq.(1) is given by E ∈ \ a (cid:16) µ a − t + w a ) , µ a + 2( t + w a ) (cid:17) (15)We can roughly call the above method as “energy matching method”. Now we can apply thisformula to determine the mobility edge of AA model in the following 4 different situations. AA model with disordered hopping
For the case of β = 0, µ a is just a constant µ . Then the intersection of all the aboveintervals is just the smallest possible interval for certain a . It is easy to see that the smallestinterval corresponds to w a = − w . Therefore, we find that the energy region of extendedstate is E ∈ (cid:16) µ − t − w ) , µ + 2( t − w ) (cid:17) (16)which recovers the results we have obtained in the last section. AA model with the same disorders in hopping and potential
Now we turn to a more complicated model with β = α . In other word, µ a and w a havethe same disorder. We have seen that the semi-analytic method is not very useful in thiscase. We will see that our energy matching method provides us a better way to determinethe mobility edges. For a given µ , the smallest interval is again given by taking w a = − w .Since µ a = µw w a , we also have µ a = − µ at the same time. Therefore, the energy region ofextended state should be given by E ∈ (cid:16) − µ − t − w ) , − µ + 2( t − w ) (cid:17) (17)The other extreme values of w a = w and µ a = µ will give us another lower bound E min = µ − t + w ). When µ < w , this lower bound is too low comparing with the one from7q.(17). But when µ > w , this lower bound will be higher than the previous one, whichbecome a new lower bound. In summary, when β = α , we find the condition for extendedstates is given by E ∈ (cid:16) − µ − t − w ) , − µ + 2( t − w ) (cid:17) , for 0 < µ < wE ∈ (cid:16) µ − t + w ) , − µ + 2( t − w ) (cid:17) , for µ > w (18)To verify the above mobility edges, we plot the eigenenergy of AA mdoel with disordersin both hopping and potential in Figure 1 (b). The parameters used in this calculationis listed in the figure caption. Again, the color of each points represents the IPR valueof the eigenstates. One can see that the dark area matches the mobility edges we havedecided above by the energy matching method. An interesting feature of this model is thatthe lower mobility edges is not monotonic. This trend is not easy to get in the traditionalmethod. In Figure 1 (d), the IPR values as a function of eigenenergy is shown for a fewpoints µ = 0 . , . , .
1. It provides a more quantitative evidence of the mobility edges ofEq.(18).
AA model with different disorders in hopping and potential
As a more exotic example, we can have different disorders in hopping and potential. Thiscan be easily achieved by setting φ = 0. In this case, w a and µ a have relatively independentdisorder, but they are not completely independent, and there is a phase difference betweenthem. For the convenience of demonstration, here we choose w a = w cos(2 παa v ) and µ a = µ cos(2 πβa v + π/
2) = µ sin(2 πβa v ), and the phase difference between them is φ = π w a and µ a are no longer the same disorder, we use the auxiliary angle formula to furthersimplify the Eq.(15) E ∈ \ a (cid:16) R − a − t, R + a + 2 t (cid:17) (19)Where R ± a are given by R ± a = R cos( ± παa v − arctan( µ w )) , R = p µ + (2 w ) (20)For a given R , the smallest interval is obtained by taking R − a = R for lower bound and R + a = − R for upper bound. Therefore, we simply find that the region of extended states isgiven by E ∈ (cid:16) R − t, − R + 2 t (cid:17) (21)8 a) (b) -5 -3 -1 1 3 5 Eigenenergy -4 -3 -2 -1 I P R (w, )=(0.5,0)(w, )=(0.5,1)(w, )=(0.5,2) =0.5 =0.5v=0.5 (c) -4 -3 -2 -1 0 1 2 3 4 Eigenenergy -4 -3 -2 -1 I P R ( , )=(1.5,0)( , )=(1.5,0.2)( , )=(1.5,0.4)( , )=(1.5,0.6) =0.5 v=0.5 (d) Figure 2. (a) Eigenenergy of AA model Eq. (1) as a function of µ with α = 0 . β = 0 . w = 0 . v = 0 . φ = π/
2. (b) Eigenenergy of SSH model Eq. (22) as a function of w with λ = 1 . β = 0 . v = 0 .
5. The colors in panels (a) and (b) represent the IPR of each eigenstates. Thered and green curves represent the mobility edges. (c) The IPR as a function of eigenenergy for( w, µ ) = (0 . , , (0 . , , (0 . , λ, µ ) = (1 . , , (1 . , . , (1 . , . , (1 . , . To verify the above mobility edges, we plot the eigenenergy of this model as a function µ inFigure 2 (a). The IPR values are reflected by the colors of each points. One can see that thedark area matches the mobility edges of Eq.(21). In Figure 2 (c), we plot the IPR values asa function of eigenenergy for a few points µ = 0 , ,
2. It provides a more precise numericalevidence of the mobility edges of Eq.(21).
SSH model with disordered hopping
9n addition to AA model, SSH model is also an extensively studied 1D topological model[34–37]. When the slow varying quasi periodic disorder is applied to the SSH model, mobilityedge will also appear [38]. Because this kind of model is more comlicated than AA model,the semi analytical method is relatively difficult, and the advantages of the new method aremore obvious. The Hamiltonian of SSH model is given by H = N X i =1 ( t − w i )( c † i,A c i,B + c † i,B c i,A ) + N − X i =1 ( t + w i )( c † i,B c i +1 ,A + c † i +1 ,A c i,B ) (22)Here w i = w cos(2 παi v ) and assume that t , = t ± λ . We approximate this model by anensemble of models without disorder as H a = N X i =1 ( t − w a )( c † i,A c i,B + c † i,B c i,A ) + N − X i =1 ( t + w a )( c † i,B c i +1 ,A + c † i +1 ,A c i,B ) (23)Here a = 1 , · · · , N . For a fix a , then energy band is given by E = ± p ( t − w a ) + ( t + w a ) + 2( t − w a )( t + w a ) cos k (24)Since the energy band is symmetric about E = 0, we can focus on the upper band, the lowerone is the same. To determine the mobility edges, we have to distinguish the following twocase.(1) When λ < t − w , the energy band is located inside the following region E ∈ (cid:16) λ +2 w a , t (cid:17) . Therefore, for the disordered model, we find the condition of extended states isgiven by E ∈ \ a (cid:16) λ + 2 w a , t (cid:17) (25)Clearly, the intersection is given by the smallest interval corresponding to w a = w . Thusthe condition of extended states for this case is E ∈ (cid:16) λ + 2 w, t (cid:17) (26)(2) On the other hand, if λ > t + w , the energy band is located inside the following region E ∈ (cid:16) t, u + 2 w a (cid:17) . Then, the condition of extended states for the disordered model isgiven by E ∈ \ a (cid:16) t, λ + 2 w a (cid:17) (27)10hen the intersection is given by the smallest interval corresponding to w a = − w . Thus thecondition of extended states for this case is E ∈ (cid:16) t, λ − w (cid:17) (28)To verify the mobility edges of this case, we plot the eigenenergy of the SSH model Eq.(22)as a function w in Figure 2 (b). The parameters used in this calculation is listed in the figurecaption. The color of each point represents the IPR value of the eigenstates. The larger IPRvalues correspond to brighter colors. One can see that the dark area matches the mobilityedges of Eq.(28). In Figure 2 (d), we also plot the IPR values as a function of eigenenergyfor a few points µ = 0 . , . , .
6, which verify the mobility edges more precisely.For completeness, we mention that there is no extended states for t − w < λ < t + w .This can be easily verified numerically (Not shown in the Figure). IV. DISCUSSION AND AN INDEX OF LOCALIZATION
The existence of mobility edge in the models with slowly varying quasi-periodic disordersis due to the coexistence of periodicity and disorders in the same model. In Eq.(3), thedisorder term w i tends to be constant when lattice site i is large enough. This representsthe periodicity of the model which favors an extended wave-function. On the other hand, forsmall lattice site i , the term w i is quite random, which favors the localized wave-function.With these considerations, one would expect that how fast the w i term approaching a stablevalue is likely to affect the degree of localization. When w i approaches to a stable value ata rapid speed, the periodicity of the model is dominant, and most of the eigenstates of themodel should be extended states. When the w i approaches to the stable value at a veryslow speed, most of the eigenstates of the model should be localized states. If it is slow toa certain extent, the mobility edge used to separate the localized state from the extendedstate will disappear.In the previous section, we have approximate the disordered model by an ensemble ofperiodic models M a for a = 1 , · · · , N . The region of extended states of the disordered modelis obtained by the overlap of the energy bands of all these periodic models M a . This methodactually suggest us a new index to represent the degree of localization at energy E . Fora given energy E , we can find the number of models M a whose energy bands contains E .11 a) (b) -4 -2 0 2 4 Eigenenergy -4 -3 -2 -1 I P R =0.5 =0.5=0.5 v=0.5w=0.5 (c) -4 -2 2 4 Eigenenergy -4 -3 -2 -1 I P R =0.5 =2.5v=0.5 w=0.5 (d) Figure 3. (a) Eigenenergy of AA model Eq. (1) as a function of µ with α = 0 . β = 0 . w = 0 . v = 0 . φ = π/
2. (b) Eigenenergy of SSH model Eq. (22) as a function of mu with λ = 1 . β = 0 . v = 0 .
5. The colors in panels (a) and (b) represent the index of localization of eacheigenstates. The red and green curves represent the mobility edges. (c) The IPR as a functionof eigenenergy for ( w, µ ) = (0 . , , (0 . , , (0 . , λ, µ ) = (1 . , , (1 . , . , (1 . , . , (1 . , . If this number equals the lattice site number N , we claim that E belongs to the extendedregion. Therefore, this number actually reflect how extended this energy level E is. Thedifference between N and this number reflect how localized energy level E is. With theabove consideration, one can see that the index of localization I ( E ) can be calculated as12ollows I ( E ) = N − N X a =1 Z π dk δ (cid:16) E − E a ( k ) (cid:17) (29)For the AA model, since the energy bands of model M a is E ( k ) = µ a − t + w a ) cos k , wefind I ( E ) = N − N X a =1 Z π dk δ (cid:16) E − [ µ a − t + w a ) cos k ] (cid:17) (30)Now we show some numerical results about this new localization index. In Figure 3(a) and (b), we plot the eigenenergy of AA model and SSH model as a function of µ and w respectively. The parameters used in these calculations are the same as the one usedin Figure 2 (a) and (b). The only difference is that the colors in Figure 3 represents thelocalization index we have just defined. Comparing to Figure 2, one can see that the newindex can also indicate the transition from localized states to extended states, just as IPRdid. To make a more precise comparison, in Figure 2 (c) and (d), we plot IPR as a functionof eigenerngy for the AA and SSH model for a fixed µ and w respectively. In the same time,we use color to show the value of the localization index. One can see that the smaller valuesof IPR correspond to darker colors or smaller values of I ( E ) and the larger values of IPRcorrespond to brighter colors or larger values of I ( E ). This verifies that the new index canalso indicate the degree of localization. V. CONCLUSION
In this paper, we have presented a so called “energy matching” method to study the mo-bility edges of several 1D tight-binding models with slowly varying quasi-periodic disorders.In this method, the disordered model is approximated by an ensemble of periodic models.Then the region of extended states is determined by the overlap of the energy bands of allthese periodic models. This method provides us a very efficient way to characterize themobility edges in a model with complicated quasi-periodic disorders. Based on this method,we also proposed a new index which can help us to visualize the degree of localization ofeigenstates. 13
CKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China underGrant No. 11874272 and Science Specialty Program of Sichuan University under Grant No.2020SCUNL210. [1] P. W. Anderson, Phys. Rev. , 1492 (1958).[2] N. Mott, J. Phys. C , 3075 (1987).[3] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. , 673 (1979).[4] P. A. Lee and T. V. Ramakrishnan, Rev. Mod. Phys. , 287 (1985).[5] Y. Hashimoto, K. Niizeki, and Y. Okabe, J. Phys. A , 5211 (1992).[6] D. J. Boers, B. Goedeke, D. Hinrichs, and M. Holthaus, A , 063404 (2007).[7] J. Biddle, B. Wang, D. J. Priour Jr, and S. Das Sarma, Phys. Rev. A , 021603 (2009).[8] S. Lellouch and L. Sanchez-Palencia, Phys. Rev. A , 061602 (2014).[9] S. Ganeshan, J. H. Pixley, and S. Das Sarma, Phys. Rev. Lett. , 146601 (2015).[10] H. Yao, H. Khoudli, L. Bresque, and L. Sanchez-Palencia, Phys. Rev. Lett. , 070405 (2019).[11] Y. Liu, X.-P. Jiang, J. Cao, and S. Chen, Phys. Rev. B , 174205 (2020).[12] S. Aubry and G. And´ r , Ann. Israel Phys. Soc. , 18 (1980).[13] P. G. Harper, Phys. Soc. A , 874 (1955).[14] J. Biddle, D. J. Priour, B. Wang, and S. Das Sarma, Phys. Rev. B , 075105 (2011).[15] X. P. Li, J. H. Pixley, D. L. Deng, S. Ganeshan, and S. Das Sarma, Phys. Rev. B , 184204(2016).[16] X. Li, X. P. Li, and S. Das Sarma, Phys. Rev. B , 085119 (2017).[17] X. Deng, S. Ray, S. Sinha, G. V. Shlyapnikov, and L. Santos, Phys. Rev. Lett. , 025301(2019).[18] M. Saha, S. K. Maiti, and A. Purkayastha, Phys. Rev. B , 174201 (2019).[19] Yucheng Wang, Xu Xia, Long Zhang, Hepeng Yao, Shu Chen, Jiangong You, Qi Zhou, andXiong-Jun Liu, Phys. Rev. Lett. , 196604 (2020).[20] Alexander Duthie, Sthitadhi Roy, and David E. Logan, arXiv:2012.01450.
21] S. Das Sarma, S. He, and X. C. Xie, Phys. Rev. Lett. , 2144, (1988).[22] S. Das Sarma, S. He, and X. C. Xie, Phys. Rev. B , 5544, (1990).[23] S. Ganeshan, K. Sun, and S. Das Sarma, Phys. Rev. Lett. , 180403 (2013).[24] Devakul and D. A. Huse, Phys. Rev. B , 214201 (2017).[25] J. Biddle and S. Das Sarma, Phys. Rev. Lett. , 070601 (2010).[26] Tong Liu and Hao Guo, Phys. Rev. B , 104201 (2018).[27] Tong Liu and Hao Guo, Phys. Rev. B , 174207 (2017).[28] Fangli Liu, Somnath Ghosh, and Y. D. Chong, Phys. Rev. B , 014108 (2015).[29] J. C. C. Cestari, A. Foerster, and M. A. Gusm˜ a o, Phys. Rev. B , 205441 (2016).[30] W. Su, J. Schrieffer and A. J. Heeger, Phys. Rev. Lett. , 1698 (1979).[31] D. J. Thouless, Phys. Rep. , 93 (1974).[32] M. Kohmoto, Phys. Rev. Lett. , 1198 (1983).[33] M. Schreiber, J. Phys. C , 2493 (1985).[34] Myles Scollon and Malcolm P. Kennett, Phys. Rev. B , 144204 (2012).[35] Linhu Li, Zhihao Xu and Shu Chen, Phys. Rev. B , 085111 (2014).[36] B. Zhu, R. Lu and S. Chen, Phys. Rev. A , 062102 (2014).[37] Hsiu-Chuan Hsu and Tsung-Wei Chen, Phys. Rev. B , 205425 (2020).[38] Tong Liu and Hao Guo, Physics Letters A , 3287 (2018)., 3287 (2018).