Local formula for the Z 2 invariant of topological insulators
LLocal formula for the Z invariant of topological insulators Zhi Li
1, 2 and Roger S. K. Mong
1, 2 Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, United States Pittsburgh Quantum Institute, Pittsburgh, Pennsylvania 15260, United States
We proposed a formula for the Z invariant for topological insulators, which remains valid withouttranslational invariance. Our formula is a local expression, in the sense that the contributions mainlycome from quantities near a point. Using almost commute matrices, we proposed a method toapproximate this invariant with local information. The validity of the formula and the approximationmethod is proved. I. INTRODUCTION
One of the most important progresses of condensedmatter physics in recent years is the realization ofmany topological phases of matter beyond the Landau-Ginzburg paradigm. While the general classification oftopological phases is still in progress, the classificationfor gapped non-interacting fermions is well-established[1–4] and shows beautiful connections to K -theory andsymmetric spaces. According to the action of several dis-crete symmetries, systems are classified into 10 classes.In each class, systems are labeled by a topological invari-ant valued in Z or Z . The pattern appearing for variousdimensions can be naturally explained by the Bott peri-odicity [5, 6] and can be arranged into a periodic table.A topological insulator, first proposed by Kane andMele in Ref. 7, is a nontrivial system in two-dimension(2D) with time-reversal symmetry which squared to − Z number whichwe call Kane-Mele invariant.For systems with translational invariance, one canget analytical formulas for the topological invariants byworking in momentum space and considering essentiallysome vector bundles (with symmetries) [2–4] over theBrillouin zone. For example, see Refs. 7, 9–13 for var-ious formulas for the Z Kane-Mele invariant.While the classification is believed to be robust againstdisorder [3, 14, 15], analytical formulas are more difficultto find. Nevertheless, one can still get useful results fromnoncommutative geometry/topology considerations [16–19], which may manifest itself as a (Fredholm, mod 2Fredholm, Bott, etc) index [20–29] (although some ofthem are abstract definitions and do not tell us how tocalculate them efficiently); or from physical considera-tions such as scattering theory [30]. A nice example isthe following formula [31, 32] for two-dimensional Cherninsulators (class A): ν ( P ) = 12 πi (cid:88) j ∈ A (cid:88) k ∈ B (cid:88) l ∈ C ( P jk P kl P lj − P jl P lk P kj ) , (1)where P is the orthogonal projection operator onto filledbands, or equivalently the ground-state correlation ma- trix. { A, B, C } is a partition of the plane into three parts,as in Fig. 1(b). This formula reveals the local nature ofthe Chern number: assuming P ij decays fast enough as | i − j |→ ∞ , then v ( P ) can be well approximated by onlysumming over j, k, l near the intersection point. For ex-ample, truncate the plane with a circle as in Fig. 1(c),then the same summation (with A, B, and C now finite)provides a good estimation. Φ (a) A BC l l l ⊗ (b) A BCD ⊗ (c) FIG. 1. (a) Insert a flux in the hole. (b) Divide the plane intothree regions A, B, C. The intersection point is where a fluxwill be inserted. (c) Truncate the plane with a circle. Denotethe region outside the circle by D. The intersection point iswhere a flux will be inserted.
In this paper, we propose a formula for the Z invariantfor topological insulators in two dimensions, which re-mains valid with disorder. Importantly, our approach is purely topological , in the sense that we discard many geo-metrical information/choices such as distances and angles[see Eq. (16)]. Moreover, we only require a mobility gapinstead of a spectral gap. Similar to Eq. (1), the input ofour formula is the projection P . Also similar to Eq. (1),our formula is essentially a local expression, in the sensethat the contribution mainly comes from quantities neara point. As a result, one can expect to calculate it withsufficient precision by a truncation near that point.This paper is organized as follows. In Sec. II, we ex-plain the physics intuition and give a physical derivationof our formula. In Sec. III, we formally state our formulaand show that it is well-defined. In Sec. IV, based onthe theory of almost commuting matrices, we introducea method to numerically calculate the invariant from afinite-size system. We present some numerical results inSec. V. In Sec. VI, we investigate the properties of ourformula and sketch the proof of our main proposition. Tokeep the paper more accessible, some technical details aregathered into the Appendix A. a r X i v : . [ c ond - m a t . m e s - h a ll ] N ov II. INTUITION–FLUX INSERTION ANDTOPOLOGICAL INVARIANT
In this section, we put Chern insulator/topological in-sulator on a punctured plane and insert fluxes at theorigin [see Fig. 1(a)]. We will explain how the physicsof flux insertion is related to topological invariant. Thissection aims to explain our intuition and provide a phys-ical derivation of our formula, hence, some statementshere may not very rigorous. We will establish our resultscarefully in the following sections.Recall the simple case, Chern insulators, which canbe realized in integer anomalous quantum Hall systems.In this case, we have the well-known Thouless chargepump [33]: when a flux unit is adiabatically inserted, itinduces an annular electric field, which in turn produces aradial electric current due to the Hall effect. As a result,electrons are pushed away from (or close to, dependingon the sign of the current) the origin for “one unit”. InFig. 2 we draw the band structure for boundary states(near the puncture). Diagrammatically, when a flux unitis inserted, every occupied state moves toward top rightto a lower level. insert 1-flux k k k k E F E F k k k k FIG. 2. Band for boundary states of a Chern insulator. • means filled, ◦ means empty. After a unit flux insertion, everyfilled state moves towards top right to the next level. In thisprocess, the label k i is tight to the electron, not the level. The many-body state after the flux insertion is notthe ground state, because there is an filled state abovethe Fermi level. Compared to the ground states, we cansee that the new ground state has one less electron ( k electron in Fig. 2) than the old one (note that we aredoing ∞ − ∞ , see comments below). The difference ofnumber of electrons in ground states is exactly the Chernnumber. This is the idea behind Ref. 20:Chern number = Ind( P, P (cid:48) )= dim Ker( P − P (cid:48) − − dim Ker( P − P (cid:48) + 1) , (2)where P/P (cid:48) is the projection operator onto filled statesbefore/after the flux insertion, Ind is the relative indexfor a pair of projections, which intuitively counts the dif-ference of their ranks (dimension of eigenvalue 1 sub-space, number of filled levels in physics). Since the rankis just Tr( P ) and Tr( P (cid:48) ), one may expectInd( P, P (cid:48) ) ∼ Tr( P − P (cid:48) ) . (3)This formula is indeed correct if Tr( P − P (cid:48) ) is well-defined—if ( P − P (cid:48) ) is trace class [34]. This is not the case for nontrivial Chern insulator though: Ind( P, P (cid:48) ) isstill well-defined [20], but one needs a more complicatedformula to evaluate it, which is essentially Eq. (1).Now we turn to topological insulators. In this case, weadiabatically insert a -flux quanta. As shown in Fig. 3,what happens is: energies for left-movers increase, whileenergies for right-movers decrease. If we assume (withoutloss of generality) the Fermi level is right below an emptystate, the ground state after the flux insertion will haveone more electron than before. We want to count thenumber of extra electrons to determine the Kane-Meleinvariant ν KM (mod 2). E F insert 1/2-flux E F k k k k k k k k k k k k k FIG. 3. Band for boundary states of a topological insulator.After a half flux insertion, we get one more state under theFermi level, which is geometrically near the vertex (flux).
To do this, we first count the number of electrons in afinite disk with radius r (the system is still on an in-finite plane, we just draw a virtual circle to define adisk). Due to time reversal symmetry, topological in-sulator have zero total Hall conductance, so the numberof electrons inside the disk remains unchanged under adi-abatic flux insertion. However, there is a vertex state ( k R in Fig. 3) that is left empty, so in the new ground statethe number of electrons in disk is increased by 1:∆ (cid:104) No. electrons in a (large) disk (cid:105) = 1 = ν KM (mod 2) , (4)where (cid:104)(cid:105) means ground state expectation value (noteagain that ground states before and after the flux in-sertion are different).Since P is the projection onto filled states, H =1 − P can be regarded as a spectral-flattened Hamiltonian(filled= 0, empty= 1). Denote H to be the Hamiltonianafter flux-insertion, consider Q = 1 − H and correspond-ing projection matrix Q (see Sec. III for details). We have (cid:104) N r (cid:105) (before) = (cid:104) (cid:88) | i | No. electrons in a (large) disk (cid:105) = lim r →∞ Tr( Q r − Q r ) . (7)However, there will always be an electron go into (or outof) the disk adiabatically, which compensates the lost (orextra) state, so ∆ (cid:104) No. electrons in disk (cid:105) in the left-handside is always 0 in this case. This can also be seen fromthe (large) gauge equivalence between the two systemsbefore and after a unit flux insertion. For the right-handside, since Q = P (cid:48) in this case is already a projection, Q = Q , so the r.h.s is 0. The difference between topo-logical insulators and Chern insulators is as follows: inthe former case the number of electrons go through theboundary r is 0 in average (because of zero Hall conduc-tance), while in the latter it is nonzero and is essentiallythe Chern number.As a side note, one may also consider the insertionof a unit flux and consider the difference between twoground states. A direct application of the relative indexEq. (2) gives 0. However, one may note that two terms(dimension of the kernel) in Eq. (2) come from left moversand right movers separately and one can therefore definethe Z index with one kernel. This is the idea behindRef. 27. III. FORMULA FOR INFINITE SYSTEM In this section, we will carefully define the quantities inour main formula Eq. (6) and show its well-definedness.The input of our formula will be the single-body pro-jection operator P , which is related to the spectrum-flattened Hamiltonian H = 1 − P . For gapped states, P decays at least exponentially [35]: P x , y < C e − C | x − y | . (8)According to the Peierls substitution [36], if we insert a1/2-flux at a vertex, the new (single-body) Hamiltoniancan be written as H = 1 − Q , where Q x , y = s x , y P x , y . (9)Here s x , y are phases such that for any loop l = ( x , x , · · · , x n = x ), we have s l def = n (cid:89) i =1 s x i , x i +1 = (cid:40) − , if the vertex is in the loop1 , otherwise . (10)These phases can be chosen as follows: we divide theplane into three regions, as in Fig. 1(b).Let n x , y to be the number intersections of the straightline segment ( x , y ) with three boundaries l , l , l , set s x , y = ( − n x , y . (11)We call this gauge “insert half fluxes along the bound-aries”. While P is a projection, Q no longer is. Actually,we have (12)( Q − Q ) x , y = (cid:88) z s x , z s z , y P x , z P z , y − s x , y P x , y = − s x , y (cid:88) z (cid:48) P x , z P z , y . where (cid:80) (cid:48) means sum under constraint s xyz = − 1. De-note V = Q − Q . Since matrix elements of P decaysexponentially, V is mainly supported around the vertex(hence the notation V ) due to the constraint. To be spe-cific, we have the following: Proposition III.1. ∃ C (cid:48) , C (cid:48) , such that | V x , y | < C (cid:48) e − C r where r = max {| x | , | y |} . Proof. Let us calculate V x , y : | V x , y | = | (cid:88) z (cid:48) P x , z P z , y | < C (cid:88) z (cid:48) e − C ( | x − z | + | z − y | ) . (13)From geometry, it is obvious that | x − z | + | z − y | > r if s xyz = − 1, so the summation can be controlled by2 C e − C r (cid:88) z e − C ( | x − z | + | z − y | ) < C (cid:48) e − C (cid:48) r . (14)(This is a pretty crude estimation but is enough.)In the following, we will refer call property as “expo-nential decay property” (EDP). Intuitively, Q − Q sat-isfies EDP means the deviation of Q from a projectionmainly comes from states near the vertex point. If wespectral flatten Q to Q (for eigenvalues λ ≤ , convert itto 0, otherwise convert it to 1), we anticipate that Q − Q is mainly supported near the vertex. Actually Q − Q alsoobeys EDP, but we do not need this result. We only needthe following: Proposition III.2. Q − Q is trace class. Proof. | x − ¯ x |≤ | x − x | for ∀ x ∈ R , so | Q − Q |≤ | Q − Q | = 2 | V | as an operator (note that they commute). Since V obeys EDP, V must be trace class (see the corollaryafter Lemma 2 in appendix A 1), so is Q − Q .Therefore, it is legal to define a “trace over vertexstates” as Tr v ( Q ) = Tr( Q − Q ) . (15)Note that in the definition of Q we can arbitrarily choosethe chemical potential µ ∈ (0 , v ( Q ) should benaturally understood as mod 1. In the case of topologicalinsulator, Q has time reversal symmetry, every states isKramers paired, so Tr v ( Q ) can be naturally understoodas mod 2. We will see in the following that it is Tr v ( Q )mod 2 [instead of Tr v ( Q ) itself for a fixed “chemical po-tential”] that has good properties. Also note that Q isnot trace class in general, so we cannot define Tr v ( Q ) asTr( Q ) mod 2.According to the above analysis, this expression is well-defined and the contributions mainly come from statesnear the vertex; it is a local expression. Interestingly,this local expression turns out to be independent of theflux-insertion point we choose; it only depends on thestate itself. Moreover, it is an integer and is topologicallyinvariant. Our main proposition is as follows: Main Proposition Tr v ( Q ) equals to the Kane-Mele in-variant.The derivation of our proposition is in Sec. VI and theappendix A. Before going on, we give three comments onour formula.1. There is another construction of Q , closely relatedto the one given by Eq. (10): Q = AA − AB − AC − BA BB − BC − CA − CB CC , (16)where AB means P AB , a block in the original ma-trix P . If we consider a circle with many sites onit, it still gives us a total phase − 1. In this case, Q − Q = 2 ACB ABCBCA CABCBA BAC , (17)where ACB means P AC P CB , etc. It is still concen-trated near the vertex (satisfies EDP), as long asthe partition is good . So we can follow the sameprocedure to define a new Tr v ( Q ).The Q defined here is not unitary equivalent to theone in Eq. (10)—they have different spectra in gen-eral. However, in Sec. VI we will show that Tr v ( Q )defined from them are equal (mod 2). We call the Q in Eq. (16) the topological one, denoted by Q t , be-cause its definition does not depend on the geomet-ric information such as “straight line segments”. For example, the one in Fig. 1(b) is good. However, if we rotate l towards l and deform it a little bit so they are parallel at infinity,then Q − Q does not satisfy EDP and convergence problem willoccur. The Q in Eq. (10) will be called the geometric one,denoted by Q g . It has the advantage of gauge in-variance and many quantities [like spec( Q )] definedfrom it are manifestly independent of the partition.2. For systems in the DIII class (TRS = − 1, PHS =1 where TRS is time-reversal symmetry and PHSis particle-hole symmetry), our formula can be sim-plified.Indeed, since the original system has PHS: K † ph (2 P − K ph = 2 P − , (18)where (2 P − 1) is the spectral-flattened Hamiltonianwith spectrum= {± } . K ph is an onsite action, andcommutes with the operation from P to Q , so thesame equation holds for Q . Therefore the spectraof Q is symmetric with respect to 1 / σ ( Q ) = 1 − σ ( Q ) . (19)Now, for a spectrum q such that q (cid:54) = , theKramers degeneracy and PHS provides us a four-fold { q, q, − q, − q } , which contributes 0 toTr v ( Q ) (mod 2). So ν = Tr v Q = No. { Kramers pairs at } (mod 2) . (20)3. The input P is the correlation matrix for an infi-nite system. If we start with a finite system, saya topological insulator on the sphere, then our for-mula always gives 0.Mathematically, this is because both Tr( Q ) andTr( Q ) = Tr( P ) are even due to time-reversal sym-metry. Physically, it is because when inserting aflux at some point, it is unavoidable to insert an-other flux at somewhere else for a closed geometry,then our formula counts the vertex states at bothpoints. To get the right invariant, we need to “iso-late” the physics at one vertex. IV. APPROXIMATION FROM FINITE SYSTEM Although the input of our formula is an infinite-dimensional operator P , our formula is a trace of vertexstates, which should only depends on the physics nearthe origin. Let us truncate the plane with a circle r [seeFig. 1(c)]. Denote P N and Q N to be the truncation of P and Q , where N is the number of sites inside the circle, N ∼ r . We expect that one can approximate the invari-ant with data near the origin, i.e. with matrix elementsof P N or equivalently Q N .However, a naive limit lim N →∞ Tr( Q N − Q N ) is wrong:it will give lim Tr( P N ) (mod 2) since Tr( Q N ) is even.This is because Q N (cid:54) = ( Q ) N . Physically, Q N and Q dohave similar “vertex states”. However, unlike Q , Q N also includes boundary contributions [see Eq. (22)], whichneed to be excluded.We claim that we can use the following algorithm toapproximate our invariant. • Construct a matrix V N by V N = − s x , y (cid:88) z ∈ ( ABC ) s xyz = − P x , z P z , y . (21) V N will be almost commuted with Q N and it willtell us whether a state is near the vertex or theboundary. • Find approximations Q (cid:48) N , V (cid:48) N for Q N , V N so thatthey indeed commute . • Simultaneously diagonalize Q (cid:48) N and V (cid:48) N to get pairsof eigenvalues ( q (cid:48) , v (cid:48) ). Sum over all the eigenvalues q (cid:48) such that v (cid:48) (cid:54) = 0.The summation will converge to Tr v ( Q ) as N → ∞ .In the following we explain the algorithm in detail.First of all, we have:( Q N − Q N ) x , y = − s x , y [2 (cid:88) z ∈ ABCs xyz = − P x , z P z , y + (cid:88) z ∈ D P x , z P z , y ] def = ( V N ) x , y + ( W N ) x , y . (22)Here, V N is supported near the center, while W N is sup-ported near the boundary. This means the deviation of Q N to a projection happens both near the vertex and theboundary.We can also work in the topological construction of Q .In this case, Q N − Q N = 2 (cid:104) ACB ABCBCA CABCBA BAC (cid:105) − (cid:104) ADA ADB ADCBDA BDB BDCCDA CDB CDC (cid:105) = V N + W N . (23)In both constructions, Q N , V N , W N should almost com-mute, and V N , W N are almost orthogonal, since they aremainly supported in different regions (“almost” meansrelevant expressions approach 0 as N → ∞ ). Proposition IV.1. (1) Q N , V N , W N defined above sat-isfies (cid:107) [ Q N , V N ] (cid:107) < (cid:15), (cid:107) [ Q N , W N ] (cid:107) < (cid:15), (cid:107) V N W N (cid:107) < (cid:15), (cid:107) W N V N (cid:107) < (cid:15), (24)where the norm (cid:107)·(cid:107) is the L norm (maximal singularvalue), (cid:15) ∼ p ( r ) e − r where p ( r ) is a polynomial of r . In practice, there are some arbitrariness to find Q (cid:48) N , V (cid:48) N . Whatwe do is a joint approximation diagonalization (JAD) and thenmake them commute according to some rules. For example, onemay make all v (cid:48) such that | v (cid:48) | > (cid:15) to zero. Another rule is indi-cated in Sec. V. (2) There exist Hermitian matrices Q (cid:48) N , V (cid:48) N , W (cid:48) N as ap-proximations of Q N , V N , W N in the sense that (cid:107) Q N − Q (cid:48) N (cid:107) < ρ , (cid:107) V N − V (cid:48) N (cid:107) < ρ, (cid:107) W N − W (cid:48) N (cid:107) < ρ, (25)such that[ Q (cid:48) N , V (cid:48) N ] = [ Q (cid:48) N , W (cid:48) N ] = V (cid:48) N W (cid:48) N = W (cid:48) N V (cid:48) N = 0 (26) Q (cid:48) N − Q (cid:48) N = V (cid:48) N + W (cid:48) N . (27)Here ρ can be chosen as F ( (cid:15) ) (cid:15) / (independent of N )where the function F ( x ) grows slower than any power of x . Proof. (1) Straight forward calculation.(2) It is easy to check || Q N || and || V N || are finite, in-dependent of N (one way to do this is to prove it forthe topological construction Q in Eq. (16) and use therelationship between two constructions as in PropertyVI.2). According to Lin’s theorem [37], ∃ Q (cid:48) N , V (cid:48) N suchthat (cid:107) Q N − Q (cid:48) N (cid:107) , (cid:107) V N − V (cid:48) N (cid:107) < δ and [ Q (cid:48) N , V (cid:48) N ] = 0.Moreover [38], we can choose δ = E (1 /(cid:15) ) (cid:15) / where thefunction E ( x ) grows slower than any power of x , inde-pendent of N .Define W (cid:48) N = Q (cid:48) N − Q (cid:48) N − V (cid:48) N , then W (cid:48) N , Q (cid:48) N , V (cid:48) N can besimultaneously diagonalized. Since W N = Q N − Q N − V N , (cid:107) Q N − Q (cid:48) N (cid:107) , (cid:107) V N − V (cid:48) N (cid:107) < δ , so we have (cid:107) W N − W (cid:48) N (cid:107) (cid:46) δ and (cid:107) V (cid:48) N W (cid:48) N (cid:107) = (cid:107) ( V (cid:48) N − V N + V N )( W (cid:48) N − W N + W N ) (cid:107) (cid:46) (cid:15) + δ ∼ δ. (28)This means for each pair of eigenvalues ( v (cid:48) N , w (cid:48) N ), at leastone of them should be smaller than √ δ . We manuallymake these eigenvalues to be 0, while fixing v (cid:48) N + w (cid:48) N .The new V (cid:48) N and W (cid:48) N would be strictly orthogonal, andstill commute with Q (cid:48) N , and still obeys Q (cid:48) N − Q (cid:48) N = V (cid:48) N + W (cid:48) N . Moreover, now (cid:107) V N − V (cid:48) N (cid:107) ∼ δ + √ δ ∼ √ δ def = ρ .Having Q (cid:48) N , V (cid:48) N , W (cid:48) N exactly commute, and V (cid:48) N , W (cid:48) N exactly orthogonal, we use them to distinguish vertexcontributions and boundary contributions. We simulta-neously diagonalize them and get triples ( q (cid:48) , v (cid:48) , w (cid:48) ). Dif-ferent contributions are then identified as follows (thereason for this identification is evident) [see Fig. 4]: • v (cid:48) (cid:54) = 0 , w (cid:48) = 0: vertex states • v (cid:48) = 0 , w (cid:48) (cid:54) = 0: boundary states • q (cid:48) = 0 or 1 , v (cid:48) = w (cid:48) = 0: bulk statesWe anticipate that the summation of q (cid:48) over vertex stateswill be a good approximation of Tr v ( Q ). Proposition IV.2 (finite size approximation) . After theabove procedure, lim N →∞ (cid:88) vertexstates q (cid:48) = Tr v ( Q ) . (29)The proof is in appendix A 2. (a) (b) (c) (d) FIG. 4. Numerical results. (a) The geometry for computing the topological invariants. The regions A, B, and C are representedby colours yellow, red, and navy blue respectively. (Region D is represented by cyan.) (b)–(c) Numerical results for theHamiltonian Eq. (30) with (b) m = 1 . m = 2 . 4. Shown in the plot are eigenvalues of V (cid:48) N vs. Q (cid:48) N . These resultsare generic; dots along the x -axis represent boundary states; dots near (0 , 0) and (1 , 0) are bulk states; dots on the parabolarepresent vertex states. (d) Numerical result from the Hamiltonian Eq. (31), which strongly breaks particle-hole symmetry. V. NUMERICAL RESULTS For our numerical results we use the Bernevig-Hughes-Zhang (BHZ) model [39] on a square lattice with Rashbacoupling and scalar/valley disorder: H = (cid:90) k ( H + H R ) d k + (cid:88) r H D ( r ) ,H ( k ) = v ( τ x σ z sin k x + τ y sin k y )+ ( m − t cos k x − t cos k y ) τ z ,H R ( k ) = r ( σ x sin k y − σ y sin k x ) ,H D ( r ) = V ( r , +) 1 + τ z V ( r , − ) 1 − τ z . (30)Here there are four degrees of freedom per site, with τ and σ acting on valley and spin space respectively.In Figs. 4(b) and 4(c), we show the computation of thetopological invariant of this model for v = t = 1, r = 1 / V ( r , ± ) sampled uniformly from theinterval [ − . , . m = 1 . 6, while Fig. 4(c) shows a trivial phase with m = 2 . q (cid:48) , v (cid:48) ) of the matrices Q (cid:48) N and V (cid:48) N respectively ( q (cid:48) is along the x -axis and v (cid:48) isalong y -axis). Recall that we have ( q (cid:48) ) − q (cid:48) = v (cid:48) + w (cid:48) and v (cid:48) w (cid:48) ≈ 0, hence points ( q (cid:48) , v (cid:48) ) either lives on theparabola y = x − x (if v (cid:48) (cid:54) = 0) or along the x -axis (if v (cid:48) =0). According to our analysis, points along the x -axisrepresent boundary states; points near (0 , 0) and (1 , q (cid:48) → − q (cid:48) ). Disorder only breaks this sym-metry very weakly. To break this mirror symmetry, we m total system size radius of ABC regionblue crosses 7 × × 13 4.5red line 25 × 25 8.9FIG. 5. The sum Tr v Q (approximating the invariant ν ) of afinite system for various m . The system sizes are shown in thetable. (For example, the green crosses shows data computedwith Fig. 4(a)’s geometry.) construct a spinful model with three valleys, H = 0 . λ + 0 . λ + (0 . λ + 0 . λ ) σ z sin k x + (1 . λ − . λ ) sin k y + (0 . λ − . λ − . λ + 1 . λ − . λ ) cos k x + (0 . λ − . λ − . λ + 1 . λ − . λ ) cos k y , (31)where λ , . . . , λ are the Gell-Mann matrices acting onvalley space. We retain the Rashba term with r = 0 . − . , . q (cid:48) , v (cid:48) ) is shown inFig. 4(d), with ν evaluating to 1.02 indicating a QSHphase.In Fig. 5, we plot the result of our formula Eq. (29)for model Eq. (30) as a function of m . (The data is com-puted for a single disorder realization.) For the compu-tation of Tr v Q , we distinguished the vertex states (fromthe bulk and edge states) by only considering points sat-isfying q (cid:48) < q (cid:48) > 1, or v (cid:48) < ( q (cid:48) ) − q (cid:48) . We see that thesystem is in the Quantum spin Hall (QSH) phase for thebulk of − (cid:46) m (cid:46) 2. As the Hamiltonian H + H R isgapless (Dirac-type) at m = 0, we expect a small sliverof metallic phase near m ≈ ξ ), theinvariant approaches 0 or 1 to the gapped phases. VI. PROPERTY AND PROOF In this section, we will investigate the property ofTr v ( Q ) and derive our main proposition step by step. Property VI.1 (gauge invariance) . Fixing the positionof the flux and working in the geometric definition, thenTr v ( Q ) does not depend on the actual partition of theplane. For example, the following partition and the or-dering of A, B, C will give the same Tr v ( Q ). A BC AB C = ⊗⊗ This is because different partitions correspond to dif-ferent gauge choices in the Pierls substitution. Indeed,fix a reference point ∗ , define u x = s ∗ , x /s (cid:48)∗ , x where s, s (cid:48) are the phases for two partitions. Since s ∗ , x s x , y s y , ∗ = s ∗ xy = s (cid:48)∗ , x s (cid:48) x , y s (cid:48) y , ∗ , we have: s (cid:48) x , y = u x s (cid:48) x , y ¯ u y (32)Thus Q (cid:48) = U QU † and they have the same spectra. Property VI.2. Tr v ( Q t ) defined from topological Q (forgood partitions) and Tr v ( Q g ) defined from geometric Q are equal (mod 2). Proof. We calculate Q t − Q g and find that( Q t − Q g ) x , y = − P x , y , x , y belong to different regionsand ( x , y ) intersects 2 boundaries0 , otherwise . (33)From geometry we can see if | x | > r , then the first con-dition is satisfied only if | y − x | > O ( r ) where r =max {| x | , | y |} . So Q t − Q g satisfies a EDP: | ( Q t − Q g ) x , y | 2) is the “vertex” term for two junctionsrespectively. As in Sec. III, V i ( i = 1 , 2) satisfies EDP forvertex i and S − S is trace class, so Tr v ( S ) is well-defined. A BC D1 20 ⊗ ⊗⊗ FIG. 6. Divide the plane into 4 regions A + B + C + D . Inserta 1/2-flux for each vertex (1 and 2). We will show that it is“equal” to insert a unit flux at the middle (point 0). As dist(1 , → ∞ , we have V V = V V → 0. In thelimiting case where V V = V V = 0 exactly, one cansimultaneously diagonalize them and at least one eigen-value for a common eigenstate should be 0. This meanseach “vertex state” of S (those states with S (cid:54) = 0 , 1) islocated at junction 1 or 2. Moreover, define Q as thematrix corresponding to a 1 / A + B + CD , Q corresponding to the in-sertion at point 2 with partition AB + C + D , then theeffect of S for a state near junction i will be close to theeffect of Q i , so one anticipates thatTr v ( S ) ≈ Tr v ( Q ) + Tr v ( Q ) . (38)In the case where dist(1 , 2) is large but not infinity,some vertex states of S may come from the coupling oftwo vertex states at different fluxes. However, it is notdifficult to convince that Eq.(38) still holds. The physicshere is similar to that for the two-states systems: due tothe weak but nonzero coupling (off diagonal elements),the eigenstates will be approximately | φ ±(cid:105) = 1 √ | (cid:105) ± | (cid:105) ) . (39)Here, we cannot say a vertex states of R is located at oneof the fluxes. However, the summation of eigenvalues for | φ ±(cid:105) is still equal to that for | (cid:105) and | (cid:105) . Property VI.3 (additivity) . Under technical assump-tions, as the distance between two fluxes goes to infinity,the vertex spectrum of S “decouples”:lim dist(1 , →∞ Tr v ( S ) − (Tr v ( Q ) + Tr v ( Q )) = 0 . (40)The proof is in appendix A 1. Property VI.4. Tr v ( S ) = 0 mod 2. This is an exactequation, no matter whether dist(1 , 2) is large or not.The idea is, if one looks from a large distance, insertingtwo 1/2-fluxes is approximately equivalent to insert a 1-flux, which is (singularly) gauge equivalent to 0-flux, sothat no vertex states appear in the spectrum at all. Proof. We works in the AB gauge, where the vector po-tential of a flux satisfies | A ( r ) | = flux2 πr , A ( r ) ⊥ r . (41)We still use S to denote the Hamiltonian in this gauge.Denote T to be the Hamiltonian for the case of inserting1-flux at the center of two half-fluxes. S and T shouldbe close outside the center. We will prove that S − T is trace class in the Appendix A 3. Then, similar to theproof of Property VI.2,Tr( S − T ) = lim N →∞ Tr( S N − T N ) = 0 , (42)since they have the same diagonal elements. S − T =( S − S ) + ( S − T ) is also trace class, soTr( S − T ) = Ind( S, T )= dim Ker( S − T − − dim Ker( S − T + 1)= 0 (mod 2) , (43) due to time reversal symmetry. Therefore,Tr v ( S ) = Tr( S − S ) = Tr( S − T ) − Tr( S − T ) = 0 (mod 2) . (44) Property VI.5. Tr v ( Q ) is an integer (mod 2) indepen-dent of the position of the flux. Proof. For every pair of points 1 and 2, we haveTr v ( Q ) − Tr v ( Q ) = [Tr v ( Q )+Tr v ( Q )] − [Tr v ( Q )+Tr v ( Q )] . (45)Due to Properties VI.3 and VI.4, it can be arbitrarilyclose to 0 (mod 2) as dist(1 , 3) and dist(2 , 3) goes to in-finity. So we must haveTr v ( Q ) = Tr v ( Q ) (mod 2) . (46)Plug this into Eq.(40) and use again property VI.4, weobtain that Tr v ( Q ) = Tr v ( Q ) ∈ Z .This property already shows that Tr v ( Q ) is an Z in-variant for topological insulators which only depends onthe state itself. The only natural identification is theKane-Mele invariant. Property VI.6. Tr v ( Q ) is equal to the Kane-Mele in-variant. Proof. Denote our invariant as ν . For two gapped time-reversal-symmetric systems A and B , we stack them(without hopping/interaction) and denote the new sys-tem A ⊕ B . Obviously ν ( A ⊕ B ) = ν ( A )+ ν ( B ). From theclassification of topological insulator [1–3] for AII class,the Kane-Mele invariant ν KM is the only invariant withthis property. It follows that ν = kν KM , (47)where k = 0 or k = 1.To prove k = 1, it is enough to verify the existence ofone system with ν = 1. To this end, consider a transla-tional invariant system in the DIII class with nontrivial Z invariant. In this case, according to Eq. (20), we have ν = No. { Kramers pairs at 12 } = 1 (mod 2) . (48)The last equation can be obtained by considering theband structure, due to translational invariance: AKramers pair at 1/2 correspond to the a band cross-ing. VII. CONCLUSION In this paper, we propose a formula Eq. (15) for the Z invariant for topological insulators in 2D, which remainsvalid with disorder. The intuition behind our formula isflux-insertion-induced spectral flow, which manifest itselfas the difference of numbers of electrons in the groundstates. The formula works by taking the single-body pro-jection matrix P (or ground state correlation function)as the input, performing a Peierls substitution (either ge-ometrical or topological), and then take the “trace oververtex states”. Our formula is a local expression, in thesense that the contribution mainly comes from quanti-ties near an arbitrarily but fixed point. The validity ofthis formula is proved indirectly, by showing its proper-ties (gauge invariance, addictivity, integrality, etc). Allproperties are physically explained and mathematicallyproved.Due to the local property of our formula, it can bewell approximated with partial knowledge of the projec-tion matrix. In this case, we construct “vertex matrix”and “boundary matrix” which almost commute. Usingan interesting parabola, the vertex contributions are sep-arated out. The validity of this algorithm is proved andverified numerically.Similar ideas may be used in the case of other symme-tries and other dimensions. For example, for class A in2D (Chern insulator), an infinitesimal flux insertion willreproduce Eq. (1). It would be interesting to work outother cases to see if one can get (and prove) a new for-mula. last but not least, the idea may also be extended tosome interacting cases. It would be interesting to exploresuch generalizations. VIII. ACKNOWLEDGEMENT This work is supported from NSF Grant No. DMR-1848336. ZL is grateful to the PQI fellowship at Univer-sity of Pittsburgh. Appendix A: More technicalities in the proof1. Proof of additivity In this appendix, we will prove PropertyVI.3 inSec. VI. The proof is a little bit technical, but the physicsidea is simple: vertex states of S comes from those of Q and Q . Lemma 1 (Almost orthogonal vectors) . For N unit vec-tors u n in a d -dimensional linear space s.t. | ( u i , u j ) | < σ for each i (cid:54) = j . If σ < √ d , then N < d − Proof. Let A = ( A i,j ) = (( u i , u j )) to be the Gram ma-trix of { u i } . Denote λ , · · · , λ N the eigenvalues of A .Since u i ∈ C d , rank( A ) ≤ d , at most d of them arenonzero. By Cauchy inequality, (cid:80) λ i ≥ ( (cid:80) λ i ) /d =(Tr A ) /d = N /d . On the other hand, (cid:80) λ i =Tr( A ) = (cid:80) | ( u i , u j ) | < N + N ( N − σ . Therefore, N d < N + N ( N − σ < N + N ( N − 1) 12 d . (A1)Solving this inequality, we find N < d − Lemma 2 (an estimation of the eigenvalue distribution) . Assume a Hermitian matrix A satisfying the exponen-tial decay property (EDP) | A i , j | < C p ( t ) e − C t where t = max {| i | , | j |} , p ( t ) is a monic (leading coefficient=1)polynomial. The number of eigenvalues outside ( − (cid:15), (cid:15) ) isbounded by O ( C ln C C α (cid:15) ), where α = 2 + deg p . Proof. For an unit eigenvector x : Ax = ax , | a | > (cid:15) , weseparate x into two parts x = y ⊕ z according to whetherthe label is inside or outside a circle: y i = 0 if | i | > r , z i = 0 if | i | < r . The radius r will be determined later(depend on (cid:15) ).We claim that (cid:107) z (cid:107) < δ def = C r / p ( r ) e − C r (cid:15) . Indeed, ay ⊕ az = Ax , || x || = 1. According to the Cauchy inequalitywe have (cid:107) (cid:15)z (cid:107) < (cid:107) az (cid:107) < (cid:88) | i | >r, j | A i , j | (cid:46) [ C r p ( r ) e − C r ] . (A2)Here, (cid:46) (means inequality up to constant) can be veri-fied by doing integral. Denote the number of eigenvalueslarger than (cid:15) to be N : Ax n = a n x n , n = 1 , · · · , N . With-out loss of generality, we can assume they are orthogonal,so ( x i , x j ) = 0 ⇒ | ( y i , y j ) | = | ( z i , z j ) | < δ . (A3)Now we have N unit vectors u n = y n √ − z n in dimension d ∝ r whose inner products with each other are lessthan σ def = δ − δ . We choose r large enough so that σ < √ d ∼ r . (A4)According to Lemma 1, N < d = O ( r ). We cansolve Eq. (A4) to estimate r (thus N ). Roughly, set σ ∼ ( C r / p ( r ) e − C r (cid:15) ) = r , let x = C r , we find e x = C C α (cid:15) x α , where α = 2+deg p . The exact solution canbe expressed using the Lambert W function [40]. Here weonly need the asymptotic expansion. Denote β = C C α (cid:15) ,take logarithm, we have the following iteration: (A5) x = ln β + α ln x = ln β + α ln(ln β + α ln x )= · · · = O (ln β ) . Thus, it is enough to choose r = O ( C ln C C α (cid:15) ), and N < d = O ( r ) = O ( C ln C C α (cid:15) ).As a corollary, it follows that (cid:88) | a | <(cid:15) | a | = (cid:88) | a | <(cid:15) (cid:90) (cid:15) θ ( | a |− x ) dx = (cid:90) (cid:15) (cid:88) | a | <(cid:15) θ ( | a |− x ) dx< (cid:90) (cid:15) C ln C C α x dx, (A6)0which converges to 0 as (cid:15) → 0. Therefore any EDP op-erator is trace class. Lemma 3. Assume A is Hermitian. If ∃ x (cid:54) = 0 s.t. (cid:107) ( A − λ ) x (cid:107) < (cid:107) (cid:15)x (cid:107) , then A has an eigenvalue a ∈ ( λ − (cid:15), λ + 2 (cid:15) ). Moreover, decompose x = x (cid:107) + x ⊥ withrespect to subspace ( λ − (cid:15), λ + 2 (cid:15) ), then (cid:13)(cid:13) x ⊥ (cid:13)(cid:13) < . If A is of finite size N × N , then an eigenvector x a of A witheigenvalue a ∈ ( λ − (cid:15), λ + 2 (cid:15) ) satisfies | ( x, y ) | > (cid:113) N . Proof. Denote B = A − λ , then (cid:107) Bx (cid:107) < (cid:107) (cid:15)x (cid:107) . Withoutloss of generality, assume (cid:107) x (cid:107) = 1. Let us diagonalize B ,so that B = diag { b , · · · , b n } . Then we have (A7) (cid:15) > (cid:88) b i | x i | = (cid:88) | b i |≥ (cid:15) b i | x i | + (cid:88) | b i | < (cid:15) b i | x i | > (cid:15) (cid:88) | b i |≥ (cid:15) | x i | + (cid:88) | b i | < (cid:15) b i | x i | . So we must have (cid:80) | b i |≥ (cid:15) | x i | < , thus (cid:13)(cid:13) x ⊥ (cid:13)(cid:13) < and ∃ i such that | b i | < (cid:15) . If N is finite, then from (cid:80) | b i | < (cid:15) | x i | < we know ∃ i such that | b i | < (cid:15), | x i | > (cid:113) N .Go back to the original proposition. We want to finda correspondence between vertex eigenvalues of S andthose of Q and Q . To make the notation simple, in thefollowing e − r means C p ( r ) e − C r and C , C can change.For each vertex state x of Q : Q x = qx , (cid:107) x (cid:107) = 1( q (cid:54) = 0 , x as x = y + z with respect to thedisk B (1 , r ). As in Eq.(A2), we still have (cid:107) z (cid:107) < e − r | q − q | .It is not difficult to show that (A8) (cid:107) ( S − q ) x (cid:107) = (cid:107) ( S − Q ) x (cid:107)≤ (cid:107) ( S − Q ) y (cid:107) + (cid:107) ( S − Q ) z (cid:107) (cid:46) (1 + 1 | q − q | ) e − r def = δ. According to Lemma 3, S has an eigenvalue in ( q − δ, q +2 δ ). The same argument also applies to Q . Thus, foreach vertex eigenvalue q of Q , Q , we get an eigenvalueof S in a neighbourhood of q .Denote U (cid:15) def = ( − (cid:15), (cid:15) ) ∪ (1 − (cid:15), (cid:15) ). For q / ∈ U (cid:15) , | q − q | > (cid:15)/ 2, so that δ < e − r (cid:15) . As r → ∞ , we willadjust (cid:15) accordingly so that δ is still small enough suchthat spectral gaps outside U (cid:15) are always greater than δ .Then we can describe the spectra structure of Q , Q inFig. 7. The shaded windows have width ∼ δ and are ofthree types. For types 1 and 2, we already get the cor-respondence. For type 3, we claim the dimension of thesubspace X corresponding to eigenvalues (of S ) in suchwindow is at least 4. Indeed, denote x , x the eigenstatesof Q , x , x the the eigenstates of Q , then | ( x i , x j ) | < δ .Let x i = u i + v i where u i ∈ X and v i ⊥ X , then | ( y i , y j ) |≤ | ( x i , x j ) | + | ( z i , z j ) | < δ + 14 . (A9) -ϵ ϵ 1+ϵ1-ϵ Q Q S FIG. 7. Spectra structure of Q i and S . We only considereigenvalues outside ( − (cid:15), (cid:15) ) ∪ (1 − (cid:15), (cid:15) ). The dots • areeigenvalues. Each dots represents a Kramers pair due to timereversal symmetry. The shaded windows are of width δ ∼ e − r (cid:15) and are (from left to right) of three types. Similar to Eq. (A1) (here N = 4), we get d ≥ S with s ∈ U (cid:15) is generatedin this way. Indeed, assume Sx = sx , (cid:107) x (cid:107) = 1 ( s (cid:54) = 0 , x = 1 s − s ( V x + V x ) (A10)is a decomposition of x . At least one of (cid:107) V i x (cid:107) shouldbe larger than | s − s | / 2, say V x . Note that V i x and( S − s ) V i x are mainly supported near vertex i , and( S − s ) V x + ( S − s ) V x = ( S − s )( S − S ) x = 0 , (A11)and both terms must be small: (cid:107) ( S − s ) V i x (cid:107) < e − r . (A12)Therefore, (A13) (cid:107) ( Q − s ) V x (cid:107) ≤ (cid:107) ( S − s ) V i x (cid:107) + (cid:107) ( S − Q ) V x (cid:107) < e − r ≤ e − r | s − s | (cid:107) V x (cid:107) def = δ (cid:48) || V x || . According to Lemma 3, Q has an eigenvalue in ( s − δ (cid:48) , s + 2 δ (cid:48) ). Therefore, s must be near (within ∼ δ ) awindow, and an argument similar to (A9) shows that x cannot be a new eigenstate.Due to this correspondence, the last term in the de-composition (A14) | Tr v ( S ) − (Tr v ( Q ) + Tr v ( Q )) |≤ | (cid:88) s ∈ U (cid:15) s | + | (cid:88) q ∈ U (cid:15) q | + | (cid:88) q ∈ U (cid:15) q | + | (cid:88) s/ ∈ U (cid:15) s − (cid:88) q / ∈ U (cid:15) q − (cid:88) q / ∈ U (cid:15) q | , is bounded by δ ln (cid:15) and goes to 0 as r → ∞ . Thesecond and third term can be bounded due to Eq. (A6),since Q i − Q i obeys EDP with the same C and C .For the first term, S − S also obeys EDP (with respectto the central point 0), however with a new constant C (cid:48) ∼ C e r . This is because (see Eq. (37)) the EDP of V i is with respect to vertex i , so the decay property of1 S − S with respect to vertex 0 need to be estimated by e − dist( x , < e r e − dist( x , . Luckily, similar to Eq. (A6),we still have | (cid:88) s ∈ U (cid:15) s | < (cid:90) (cid:15) C ln C e r C α x dx → , (A15)as r → ∞ as long as (cid:15) = o ( r ). Thus we have provedthat lim r →∞ Tr v ( S ) − (Tr v ( Q ) + Tr v ( Q )) = 0 . (A16)Note that we have assumed that the requirement (cid:15) = o ( r ) is compatible with δ = e − r (cid:15) < gap. This technicalassumption is reasonable. Indeed, according to Lemma2, the spectral gap at (cid:15) is roughly ( dd(cid:15) ln (cid:15) ) − ∼ (cid:15) ln (cid:15) in average. In order for δ < (cid:15) ln (cid:15) , it is enough to set (cid:15) > Ω( e − C (cid:48) r ), which is exponentially smaller than r forlarge r . Even if we consider the fluctuation of the spectralgaps and even if the level statistics is Poissonian (so thatno level repulsion), the probability for this to be true is1 from the following estimation: (A17)Pr(at least one gap < p ( r ) e − r (cid:15) ) < (cid:88) x>(cid:15) p ( r ) e − r /(cid:15)x/ ln x< p ( r ) e − r ln (cid:15)(cid:15) × ln (cid:15) → 2. Proof of the finite size approximation In this section, we prove Proposition IV.2. The tech-nique used will be similar to the above section. We needto compare vertex eigenvalues of Q and Q (cid:48) N . Recall that (cid:107) Q (cid:48) N − Q N (cid:107) ≤ ρ , (cid:107) V (cid:48) N − W N (cid:107) < ρ, (cid:107) V (cid:48) N − W N (cid:107) < ρ where ρ ∼ e − Cr .Temporarily fix (cid:15) , and only consider eigenvalues out-side U (cid:15) def = ( − (cid:15), (cid:15) ) ∪ (1 − (cid:15), (cid:15) ). For any α / ∈ U , | α − α | > (cid:15)/ Q − q ) x = 0, q / ∈ U (cid:15) , we separate it as x = y + z with respect to circle r/ 2, again (cid:107) z (cid:107) < p ( r ) e − Cr (cid:15) def = δ . Inthe following δ means “everything that goes like p ( r ) e − Cr (cid:15) with perhaps different C ”. Thus (cid:107) ( Q (cid:48) N − q ) x (cid:107) = (cid:107) ( Q (cid:48) N − Q N + Q N − Q ) x (cid:107)≤ (cid:107) ( Q (cid:48) N − Q N ) x (cid:107) + (cid:107) ( Q N − Q ) y (cid:107) + (cid:107) ( Q N − Q ) z (cid:107) (cid:46) δ, (A18)so Q (cid:48) N has an eigenvalue q (cid:48) ∈ ( q − δ, q + 2 δ ) with eigen-state x (cid:48) satisfies | ( x, x (cid:48) ) | > (cid:113) N (Lemma 3). This implies x (cid:48) must contain a vertex eigenstate. Indeed, if not, we have V (cid:48) N x (cid:48) = 0 so x (cid:48) = q (cid:48) − q (cid:48) W (cid:48) N x (cid:48) which is concentratednear boundary r , thus (A19) | ( x, x (cid:48) ) | = 1 | q (cid:48) − q (cid:48) | ( x, W (cid:48) N x (cid:48) ) < (cid:15) [ | ( y, W (cid:48) N x (cid:48) ) | + | ( z, W (cid:48) N x (cid:48) ) | ] (cid:46) δ, a contradiction as r → ∞ .On the other hand, if ( Q (cid:48) N − a ) x = 0( a (cid:54) = 0 , x is a vertex state: W (cid:48) N x = 0, then x = V (cid:48) N xa − a . We have( Q − a ) x = 1 a − a ( Q − Q (cid:48) N ) V (cid:48) N x = 1 a − a [( Q − Q N ) V N x + ( Q N − Q (cid:48) N ) V (cid:48) N x ] (cid:46) δ. (A20)So, Q has an eigenvalue in ( a − δ, a + 2 δ ).Now we choose r according to the same technical as-sumption above, so that there is a correspondence outsideregion U (cid:15) for ∀ (cid:15) . Then, similarly we have | (cid:88) vertex q (cid:48) − Tr v ( Q ) |≤ | (cid:88) q ∈ U (cid:15) q | + | (cid:88) q (cid:48) ∈ U (cid:15) q (cid:48) | + | (cid:88) q / ∈ U (cid:15) q − (cid:88) q (cid:48) / ∈ U (cid:15) q (cid:48) | . (A21)The last term is bounded by δ ln (cid:15) which goes to 0 as (cid:15) → 0. The first term also goes to 0 since Q − Q istrace class. The second term is bounded by No. { q (cid:48) } (cid:15) .Obviously No. { q (cid:48) } < dim Q N ∼ r , so this term alsoconverges to 0 since (cid:15) = o ( r ). 3. Proof that S − T is trace class In this section, we prove the claim used in PropertyVI.4.The first step is to figure out the decay behaviour ofthe matrix elements of S − T . According to the Peierlssubstitution [36], S ij = P ij e i (cid:82) ji A · d r , T ij = P ij e i (cid:82) ji A (cid:48) · d r , (A22)where A and A (cid:48) are the vector potential for the twoflux configurations. See Fig. 8(a) (all angles here aredirected), we have: A · d r ∝ θ, A (cid:48) · d r ∝ θ + θ , ( S − T ) x , y ∼ P x , y e i ( θ + θ ) ( e i (2 θ − θ − θ ) − . (A23)From geometry, 2 θ − θ − θ = ( α − α ) − ( α − α ).2 21 21 3 45 6 θ θθααα ααα α α r (a) (b) rr y x d d β FIG. 8. Relevant geometries in the proof. Let us calculate ( α − α ). See Fig. 8(b), we have: α − α = arctan d sin βr + d cos β − arctan d sin βr − d cos β = arctan d sin βr + d cos β − d sin βr − d cos β d sin βr + d cos β d sin βr − d cos β = − arctan d sin 2 βr − d = − d sin 2 βr + O ( 1 r ) . (A24)Due to the energy gap, | A x , y | (cid:46) e − C | x − y | . Assuming r ≥ r , we claim that e − C | x − y | | e i (2 θ − θ − θ ) − | < e − C | x − y | O ( 1 r ) . (A25)Indeed, if e − C | x − y | < r + r ) , there is nothing neededto prove. If not, then | x − y | < C ln( r + r ) < C ln r (asymptotically). In this case, from geometry, we know | β i − β j | = | θ | (cid:46) | x − y | r . Then 2 θ − θ − θ = ( α − α ) − ( α − α ) = sin 2 β r − sin 2 β r + O ( r ) will be of order O ( | x − y | r )as can be seen from Taylor expansion. Then it is easy tosee that the claim also holds.The result is (in a more symmetric fashion, ignore con-stants) as follows: the operator A = S − T satisfies thefollowing decay property: | A x , y | < | x | + | y | ) e −| x − y | . (A26) Now, we prove this kind of operator must be traceclass. Let us denote the n th singular value (decreasingorder) to be s n . According to the Courant min-max prin-ciple [34], we have s n = min Y n − max u ⊥ Y n − ( Au, Au )( u, u ) , (A27)where Y n − means a subspace of dimension n − 1. Thus,for any given n-dimensional subspace Y n − , we have s n ≤ max u ⊥ Y n − ( Au, Au )( u, u ) = max u ⊥ Y n − || u || =1 || Au || . (A28)Let us choose the subspace Y n − to be spanned by the n components nearest to the center (so that the label ofthe components are approximately in the disk of radius r ∼ √ n ). Denote the columns of A to be v x ( v x = Ae x , x ∈ Z is the label). With Eq. (A26) it is easy to show(note that here e −| x − y | means e − C | x − y | for a different C ) | ( v x , v y ) | = | ( A ) x , y | (cid:46) e −| x − y | | x | | y | . (A29)Thus, (A30) (cid:107) Au (cid:107) = (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:88) | x | >r u x v x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) = (cid:88) | x | , | y | >r ¯ u x u y ( v x , v y )= ( (cid:88) | x − y |≥ l | x | , | y | >r + (cid:88) | x − y |