Topological invariants from zero modes in systems of topological insulators (TI) and superconductors (TS)

By a regularization method, a topological charge density formula expressed by the form of Dirac delta functions was derived. From this analytic approach, it will be easy to tell how topological charges are related to zero-field points (or zero modes) and depend on the differential properties of the associated vector field in the proximity of isolated zero modes. With this result, topological integrals (winding numbers) can be turned into simple sums contributed by isolated zero modes. Through this method, we investigated the physical features of zero modes for several cases including the integral form for the Teo-Kane model in dimensional reductions, the case of Majorana zero modes bound to a vortex line in a topological superconductor, the Bernevig-Hughes-Zhang (BHZ) and its extended (EBHZ) models, an extended Haldane model, and the additivity of topological charges. In these systems, the physical implications of zero modes which reflect their singular nature were also discussed. Without advanced knowledge, we provide a handy formula for calculating topological numbers accompanied by intuitive understandings of origins of topological charges which are identical to the zeros of the vector field.


Introduction
Physicists can distinguish the existence of various ordered states of matter by symmetry-breaking, including superconductors, Bose-Einstein condensation, ferromagnets, charge density waves, and many other systems. These states can be explained by Landau's theory of phase transitions of finite order which is related to a local order parameter. However, in past decades, many new phases [1][2][3][4][5] of matter discovered in condensed-matter physics, like the fractional quantum Hall states and Majorana fermion states, do not have a local order parameter. These give rise to newly found phases of matter with nonlocal order parameters, known as topological phases. The phrase, topological phase, directs to the existence of a bulk invariant to distinguish trivial and nontrivial phases of matter even though they have the same symmetry. These topological phases cannot be identified by broken symmetries but by the topologies of the Hamiltonians of systems. In recent years, it was realized by an enhanced band theory called topological band theory which takes account of the concepts of the Berry connection and curvature in the Brillouin zone (BZ) [6]. The Chern number is defined as the integral of Berry curvature [7]. Moreover, some discrete symmetries, like time reversal and charge conjugate, exist in some certain topological insulating and superconducting phases. A topological number can classify these phases and distinguish them from the trivial ones.
The integral of the Berry curvature due to its Berry connection originated from the Bloch function of the Hamiltonian over the BZ is related to the Hall conductance of the insulator. It is identified as the Chern number of the filled band. The Hall conductance of the Chern insulator can be identified as the coefficient of the Chern-Simons action for an applied field [7]. If the Berry vector potential is globally well defined over the BZ, the Hall conductance as well as Chern number will be null. Nonzero values of Chern numbers are consequences of the nontrivial structure of the Berry connection and there exist singularities of Berry curvature at some points in the BZ. Therefore, it is necessary to calculate Chern numbers to distinguish if a TI is trivial or nontrivial.
Recently, tremendous interests have been raised by the studies of topological states of matter because of their potential applications [8,9]. Consequently, classifying physical systems in terms of their topological numbers becomes a key practical task. To the best of our knowledge, the topological number for the TI and TS system with Dirac-like Hamiltonian means the winding number of its field vector within a certain region. However, sometimes it is not easy to obtain the topological numbers by directly calculating so complicated integrals, especially over the BZ.
However, recent studies show that the zero-energy modes play crucial roles in determining the topological numbers in the pure TI and TS systems [7]. Meanwhile, the topological defects such as magnetic vortex lines and dislocations in crystals can provide a π flux to certain wave-number electronic states, and host the zero modes. The interplay of topological defects and gapless modes was therefore investigated, which led to the emerging of the periodic table for defect classifications [10][11][12].
There are other ways to link the topological number to zero modes, such as zero eigenvalues of electronic local in-gap Green's functions in the presence of impurities [10], and Majorana zero modes hosted by a vortex line in a topological superconductor [13]. In addition, the topological line defects like dislocations in 3D TI can serve as the probes of weak TI states [14,15].
Here we took a different route to obtain the topological winding number by a topological charge density formula through a special regularization method to regularize a vector field which exists zero-field points. From the analytic form of the formula, it was easy to tell that topological charges are directly related to zero modes for Dirac-like Hamiltonians and get intuitive understandings of the topological number.
In this paper, we would introduce an approach based on the regularization method to calculate topological numbers. The approach was taken to regularize a vector field to obtain the topological charge density formula in section 2. The formula was applied to various cases including the integral forms for the Teo-Kane model in dimensional reduction, the case of Majorana zero modes bound to a vortex line in a TS, the BHZ and EBHZ models, an extended Haldane model, etc in sections 3-7. Finally, section 8 is devoted to summary and conclusion. Details of the revisit of the residue theorem by the regularization method are represented in appendix A. Specially, caution is required in determining signs of the Jacobian at two of zero field points for the EBHZ model, and the readers are referred to appendix B.

Topological number of N-dimensional vector field
In the study of topological properties for vector fields, it is often required to calculate the flux through a closed boundary. Sometimes, it can be very difficult to integrate topological flux directly. In this section, we adoptd the particular regularization method mentioned above to deal with topological vector field problems. Since this method was inspired by observing the resemblance between the residue theorem and calculation of topological numbers, it would be illuminating to revisit the residue theorem before investigating the calculation of topological numbers. The details were left in appendix A.
We would like to derive the key formula to convert integrals of topological fluxes into those of topological charges with delta function forms. Let X be an N-dimensional smooth manifold, and its local coordinate is x μ (μ=1, K, N). Assume a smooth mapping  f X R : N , leading to an N-dimensional smooth manifold vector field, f a =f a (x), with a=1, 2, K, N, where R N is the real N-dimensional Euclidean space. The normalized unit vector of f(x) is expressed asf and its a-th component For simplicity, the partial derivative off a with respect to the local coordinate x μ is denoted by ¶ mf a . The topological flux density characterized by field f(x) can be expressed by^( ) f x as follows: are N-dimensional Levi-Civita symbols whose value equals 1 if their indices are an even permutation of 1, 2, K, N and equals −1 if they are an odd permutation; otherwise, it equals zero 3 . Let M be a submanifold in X and M i be the neighborhood of the i-th zero-field point, then the generalized winding number can be given by [16] This metric-free form is independent of the choice of local coordinates, which reflects the nature of topological invariance.
The divergence of the topological flux density (or topological charge density), is the integrand in equation (2). It will vanish in the region wheref ʼs are well defined (or regular) [16]. However, it is crucial to observe the fact that the direction of the vector field f is not well defined at the locations f being a zero field. As we will show later, the singular behavior of these zero-field points will be connected to the rise of topological charges. To capture this singular feature of the vector field near the zero-field points, we add an infinitesimal quantity ε 2 to the denominator off a as a regulator and the regularized vector field 2 , and its derivative can be written as e e d e ¶ = ¶ Within this regularization scheme, the new version of equation (3) is expressed as We then substitute equation (4) into equation (5) and use the antisymmetry of ò μν... and the symmetry of ∂ μ ∂ ν when two indices are exchanged to simplify the expression and obtain: The Jacobian of field components f 1 , f 2 , K, f N with respect to coordinates x 1 , in equation (6) can be replaced by  (7) is null if ¹ f 0 and goes to infinity if f=0. In addition, its integration over the whole N-dimensional vector space f is easily found as A(S N−1 ). Therefore, as e  0, the function is reduced to a generalized function: where δ(f)=δ( f 1 )δ( f 2 ) ... δ( f N ). Similarly, the limit of equation (7) is also a generalized function: where the topological charge density formula was obtained and takes the form of delta functions by the regularization method. From equation (3) to equation (9), we derived a relation:  Equation (10) was also found in a different approach [16]. However, the derivation here is more direct and simpler. Details are provided in the Supplementary Material available online at stacks.iop.org/JPCO/2/ 085014/mmedia. As we mentioned before, ∇ · J is always zero except for points at the locations f(x) = 0, where topological charges are situated. The calculation of a topological number is related to the integral of the topological field around a closed boundary and sometimes it is difficult to calculate directly. According to the divergence theorem, the flux of the topological vector field through a closed hypersurface is equivalent to the topological charge enclosed by the closed surface. In calculations of topological numbers, we can convert the associated integrals to the sums of contributions from isolated zero-field points just like the residue theorem where the contour integrals are changed into the sums of residues of isolated poles. In the following sections, this formula will be applied to several cases to show that it is simple in calculations of topological numbers for various topological systems and the physical implications of the zero modes seem clearer.

Dimensional reductions of the integral form for the Teo-Kane model
In this section, we will demonstrate the application of equation (10) to the issue raised by Teo and Kane [18] which concerned about Majorana fermionic excitations in the 3D Bogoliubov-de Gennes (BdG) Hamiltonians with particle-hole symmetry in the presence of physical defects. For convenience, it is termed as the Teo-Kane model here. The Hamiltonian of the system ( ) k r , is a function of the wavevector k and the position vector r, where k is defined in a 3D BZ (a torus T 3 equivalent to S 3 in the sense of strong topology) and r is on a 2-sphere S 2 . An adiabatic deformation  l ¢( ) k r , , parametrized by λ connects ( ) k r , at λ=0 to a trivial Hamiltonian independent of k and r at λ=1 and for more details the readers are referred to [18,19]. The Hamiltonian of the system can be expressed by where γ a and Γ a are 8×8 Dirac matrices, f a is the a-th component of the vector field f and a=1, 2, 3. These Dirac matrices satisfy {Γ a , Γ b }={γ a , γ b }=2δ ab and {Γ a , γ b }=0, and should be arranged to make the Hamiltonian possess particle-hole symmetry. The Hamiltonian  needs a six-component unit vector whered is a unit vector on -S N 2 1 . Then, the topology of the Hamiltonian  is characterized by its integervalued N-th Chern character which is calculated as the volume on S 2 N swept out by l ( ) g k r , , [18]. As we have mentioned in section 2, the N-th Chern number can be obtained from the generic formula equation (2). By replacing the local coordinate x and vector fieldf with q andĝ, respectively, the Chern number of the Majorana fermion (MF) system, W MF , can be rewritten as , , , , , , , , , , , , and θ i ʼs are S N−1 angular variables, and the subscript MF describes the system of Majorana fermions. In order to simplify equation (13), the partial derivatives ofĝ a ʼs with respect to q μ ʼs are listed in table 1. It is apparent to tell from table 1 that the nonzero terms in the integrand should take the form of ,where all the indices in the second parenthesis range from 2 to N+1, and those in the last parenthesis from N+2 to + N 2 1. The corresponding Levi-Civita symbols are consequently reduced to the form of , and the rest for (θ 1 , θ 2 , K, n m n m m n stands for combinations. Therefore, with expressionsĝ in equations (12), (13) can be rewritten as Table 1. Partial derivatives ofĝ a with respect to q μ . Here l n l n = - where all the indices, such as α i , β i , ..., range from 1 to N, and ∂ α stands for ∂/∂ k α in the first square bracket, and ∂ N =∂/∂ν, ∂ β =∂/∂ θ β , as β=1, K, N−1 in the last one. After integrating out the variables λ and ν, the integral is reduced to By replacingf a byk a and x μ by k μ in equation (10), the integration over k-variables can be performed as follows: Finally, we obtain the winding number as When N is equal to three, equation (14) reduces to ò ò By this method, the result is the same as described in [18]. For a nontrivial topological phase, the value of equation (15) should be nonzero. Since the zero-field points of thef or f field exist in real space, Thus there will be physical defects corresponding to zero-field points of the f field.

Majorana zero modes in a topological superconductor
Takahiro Fukui has studied a line defect in a 3D TS proposed by Teo and Kane which exists Majorana zero modes. The Hamiltonian of the system takes the form of equation (11) and the field vector considered for the system [13] is given by where D ¥ = D > ( ) 0 0 is required, q is vorticity of a generic vortex and the continuous function Z(z) is a mass term controlling topological phase transition. Here we assume that the mass term asymptotically is The line defect is a vortex line where the coordinates of zero modes satisfy Δ(r)=0. For simplicity, we adopt the same assumption in [13] and assume Δ(0)=0 so that the core of the vortex is located on the z-axis. To derive the topological character of the Hamiltonian , one can use equation (15) directly or the method developed by Weinberg [20]. It has been shown that the topological index of the Hamiltonian  for the system can be expressed as  can be referred to [13,21]. It is straightforward to calculate the partial derivatives of the topological field vector and derive 3 We calculate the Jacobian of f with respect to cylindrical coordinates directly and arrive at , and q is the degree of mapping. In addition, δ(f)=δ(Δ(r))δ(Θ(θ))δ(Z(z))/Δ(r), where Θ(θ)=qθ. To make the mapping between Θ and θ as a one-to-one mapping, and the integration range of θ in polar coordinates can be divided into | | q branches. Therefore, It is straightforward to rewrite d d , where we assume that z 0 is a single zero of Z. Putting all of them together and noting that Δ r 0 near r=0, where the symbol ( ) x sgn represents the sign function of x. In the presence of many zeros, the index of  will be the sum of equation (18) over all individual zero modes. If the signs of m + and m − are opposite, there will be an odd number of zeros of Z(z). In general, they are topologically equivalent to the system with a single zero because the signs of Z z for any adjacent zeros are opposite and will be canceled in pairs. Similarly, if m + and m − share the same sign, there will be an even number of zero modes in the system and the ind  will be zero. Notice in passing that Zʼs asymptotical behavior can determine the sign of Z z at a single zero point z 0 . Therefore, in general, the index for the system with Majorana zero modes can be rewritten as The result in [13] is a special case of equation (19) and independent of the detailed structure of the profile function. It is evident that a defect is a point defect situated somewhere along the z-axis domain wall, where the sign of mass profile Z(z) changes.

The BHZ and EBHZ models
Some kinds of topological systems, such as the 2D HgTc/CdTe quantum wells, respect time reversal symmetry and realize a  2 TI [6,22]. One can also add the next nearest neighbor (NNN) term into the BHZ model, which becomes the EBHZ model. This model has been studied in [10] and [12] for the investigation of the space group classification of the topological band-insulators. Its reduced Hamiltonian determining the  2 topological charges can be expressed by a 2×2 Hermitian matrix as s s s s s s sin cos 4 1 cos cos , 21 x y x y x y 1 2 3 where in equations (20) and (21) each Hamiltonian inside the first and second curly bracket represents the Hamiltonian of the BHZ model and the NNN term, respectively. The topological vector field for the EBHZ system is h= cos sin , sin sin cos , 2 2 cos cos , , and M depend on materials. The first Chern number of the system can be expressed as where òB Z M denotes the integration over the k-space and the mass coordinate, andĥ is the unit vector of h. The Jacobian in equation (22) is directly calculated as J = - cos cos sin sin cos cos . Therefore, we can rewrite Plugging equation (23) into (22) and integrating over the BZ, we finally get h v sgn v i is the sign of the Jacobian at the i-th zero-field point. They are −1, 1, 1, −1 for i=1, 2, 3, 4, respectively. The determination for signs of the Jacobian at v 1 and v 4 requires more effort, and readers are referred to appendix B for details.
We note in passing that the formalism from equation (22) to equation (24) is still valid for the BHZ model, where the vector field h=(h 1 , cos cos x y . For this model, each zero-field point, as well as its sign of Jacobian, is the same as that of the EBHZ model with z=0. Apparently, the value of M will determine the number of the zero-field points included in the integration range and therefore the value of the Chern number.
In the case of z<1/2, In particular, the phase with = - is termed as the valley phase which possesses the trivial  2 invariant and is protected by C 4 rotational symmetry [10].
The phase diagram of the BHZ model is the same as that of the EBHZ model when z=0. In this regard, our result is consistent with that of the two-band model [22]. The system can transit between different topology phases by tuning M. There are topologically nontrivial  2 quantum-spin Hall phases with W BHZ =±1 as M varies from 0 to B 8 . Outside this range, a trivial insulator phase exists. When M crosses one of endpoints of the range, a band inversion occurs and the trivial phase becomes a nontrivial one.

An extended Haldane model (EHM)
In the late 1980s, Haldane imitated the integer quantum Hall effect expected in the Landau-level problem while preserving the translational symmetry of the lattice [23]. He proposed a topological insulator on a hexagonal lattice consisting of two inter-penetrating triangular sub-lattices. The Hamiltonian contains various hopping terms In addition to a couple of hopping terms in the Haldane model, the EHM is investigated by considering one more hopping term [22]. Up to a constant term, the Hamiltonian of the system can be written as , where h i is the i-th component of the vector field h over the first BZ with a parameter M defined below. Here nearest-neighbor electron hopping with the hopping integral t 1 , next-nearest-neighbor hopping integral t 2 with a phase term f, and next-next-nearest-neighbor hopping integral t 3 are taken into account, and therefore the three components of the h field can be expressed as: where M plays the role of an on-site inversion symmetry-breaking, and a 1 and a 2 are simply the Bravais lattice vectors. When one vertex of a hexagonal lattice is chosen as the origin and the y-axis is assumed to be aligned with the longest diagonal from the origin of the lattice, then = (  The Chern number shares the same form of equation (22) with the expression of J( | ) h v in equation (26). After integrating over the k-space, we obtain the Chern number for the EHM The sign of J( | )| h v v i at the zero-field points can be found in table 2 and the system's phase diagram of the EHM is shown in figure B1. From the property of delta functions in equation (27), it is readily to see that the first zero- . Similarly, points v 2 ,v 3 and v 4 share the same value of M and the three curves associated with them coincide with each other on curve (ii), where M−msinf=0. Therefore, the multiplicity of curve (ii) is 3. Likewise, v 5 , v 6 and v 7 are on curve (iii) and v 8 on curve (iv). The value of Chern number depends on the zero-field points within the integration range which is determined by the upper limit M. For convenience, we define J ( ( | )| ) h v sgn v i as the intersection number of the curve associated with the points v i and a straight line with a constant f value. Therefore, the intersection number pertaining to curve . However, the intersection number pertaining to curve (ii) is 3 rather than 1 since the sum of h v sgn v 4 is equal to 3. Likewise, the intersection numbers of curves (iii) and (iv), should be −3 and 1, respectively. To find the Chern number of an interested point in figure B1, we draw a line with a constant f joining ¢ = -¥ M and this point. Next, we can get the Chern number by adding up the intersection numbers of the curves with the line. For instance, for the Chern number in the region between curves (i) and (ii) and −π<f<0, the above-mentioned line with one endpoint in this region will intersect with curve (i). The integration range only includes the first zero-field point v 1 and the Chern number of the region is the intersection number of curve (i), which equals to negative one. Similarly, in the region between curves (ii) and (iii) and −π<f<0, the line intersects with both curves (i) and (ii). Furthermore, the integration range includes points v 1 , v 2 , v 3 , v 4 , and the Chern number of this region is the sum of the intersection numbers of curves (i) and (ii), which is two. Compared with the Haldane model, the additional hopping term creates new nontrivial topological phases with high Chern number W EHM =±2, which means the appearance of additional edge modes [22]. By extracting the information of zero modes, it is straightforward to extend the calculation for systems with more hopping terms These systems have more degrees of freedom in the parameter space and were led to the occurrence of high-Chern-number phases.

Additivity of topological charges
Berry charges (Topological charges) due to the Berry vector potentials in the Hamiltonians of some semiconductors are closely related to their Hall conductivities [24]. In this section, we would like to demonstrate the additivity of topological charges exemplified by a bilayer graphene system and heavy holes in quantum wells. These systems are not TI. Their Hall conductivities are not quantized but are still related to (or proportional to) the topological charges of the systems.
Before discussing the additivity of topological charges, we first deal with a single charge system. The electronic Hamiltonian of the combined Rashba-Dresselhaus (RD) system [24] with the out-of-plane polarization induced by the electric field is given by where α and β are the respective spin-orbit coupling strengths, and (k x , k y ) is the electron wavevector. Here the Pauli matrices σ 1 , σ 2 , and σ 3 stand for the spin operators of electrons. The momentum dependent effective magnetic field of the system can be written as For a b ¹ | | | |, it is apparent that f has an only zero-field point at k=(k x , k y )=(0, 0), and Jacobian is We then obtain Integrating the product of equations (29) and (30) over k-space, the topological charge of the system, W RD , becomes In the case of a b = | | | |, the topological charge density in equation (9) is a zero generalized function because the Jacobian vanishes. The result is related to the spin-Hall conductivity of the system [24,25].
Now we turn to demonstrate the additivity of topological charges with following cases. First, the bilayer graphene system [24] is modeled as two coupled monolayer graphene sheets, with each layer having two inequivalent lattice sites. If the Bernal stacking configuration [26] is adopted, the Hamiltonian of the bilayer graphene (BG) system can be written as x y x y BG 2 2 2 1 2 in the low energy limit. The Hamiltonian operates on pseudospin wave functions, where the upper and lower components describe electronic amplitudes on the B-site in the top layer and the A-site in the bottom layer respectively. Only two possible ways of hopping via the dimer state are taken into account. The field vector of this system [13] can be written as . 33 x y x y 2 2 For simplicity, we add a positive number η to equation (33) to split the multiplicities, hence it becomes where η goes to zero, as f η approaches f. It is apparent that f η has two zero-field points at 2 so that its sign is always positive. Therefore, we have Finally, the topological charge of the system, W BG , is obtained as From the calculation, it is obvious that the topological charges are additive. As an alternative, the positive number η can also be added to the second term of the field vector and the resulting W BG will still be two. Hence the total topological charge is independent of this separation scheme of zero-field points. The second example is the system of p-doped quantum wells in III-V semiconductors, the energy gap between the light-and heavy-hole bands diverges with the reduction of well width. In a sufficiently narrow quantum well with low doping densities and low temperature, only the heavy-hole band is occupied The effective Hamiltonian for heavy-holes in quantum wells (QW) [24] is m  k  k k  k  k k  2  3  3  ,  3 5   y  x y  x  x y  QW   2 2  3  2  1  3  2  2 where k x and k y are the components of wave vectors, and λ is the spin-orbit coupling strength [24]. The corresponding field vector can be written as Both components of the vector f are zero at (k x , k y )=(0, 0). In order to split the multiplicities, we add a positive number η to the field vector for simplicity: When η is null, f η is the same as f. Then it is easy to see there are three zeros, , its sign is always positive. Therefore, we have x y x y x y where we break the product of a delta function and the Jacobian into the sum of three delta functions by the regularization method. Then the topological charge of the system is , 36 x y QW 0 which comes from the contribution of three individual topological charges and is related to the spin-Hall conductivity of the system. As we mentioned in the previous example, adding the positive number η to the second term of the field instead of the first term will not change the result of equation (36) because the topological charge is independent of the separation scheme of zero-field points.

Summary and conclusion
Calculation of topological numbers is crucial in identifying nontrivial phases of topological materials for potential applications. In nontrivial topological phases, there are topological obstructions depicted by zero-field points (charges), where the direction of an associated unit vector field is not well defined and singular. We have developed a simple and systematical method for calculating topological numbers and investigated the physical implications of zero modes by calculating topological numbers for systems of TS and TI. We used the regularization method to deal with the singularity of zero-field points in the BZ of TI. Within the framework of generalized functions, we obtained a topological charge density formula which can turn topological integral problems into simple sums contributed by isolated zero modes, similar to the residue theorem. Armed with this topological charge density formula, we explored the zero-mode features for various TS and TI systems. Firstly, through this formula, the dimension of the topological integral in the Teo-Kane model was lowered. This result of dimensional reduction indicates the fact that zero modes should be situated on the defects in real low dimensional space (S 2 ) for nontrivial topological phases. Meanwhile, a generalized version of dimensional reduction formula (from S 2 N to S N−1 ) was also obtained for possible applications. Secondly, Majorana zero modes in TS were identified as zero-field points in topological charge density formula. The zeros are hosted by the flux line and the domain walls of the mass function profile Z(z). The topological number of this system was obtained through the information of the identified Majorana zero modes and independent of the detailed structure of the profile function. Through our calculation, it is straightforward to see that systems with an odd number of zeros are topologically equivalent to one with a single zero and systems with an even number of zeros are topologically trivial. For the EBHZ model, one more term (the NNN term) was added to the BHZ model and an additional topological phase with a high Chern number was found which possesses the trivial  2 invariant and is protected by C 4 rotational symmetry. When the NNN term was neglected, the BHZ results were revisited. There are topologically nontrivial  2 quantum-spin Hall phases as the integration range of the mass term includes an odd number of zeros otherwise the system is trivial. After that, we investigated the EHM with more hopping terms which increases the complexity of the band structure and the possibility of zero-mode creation. From the topological charge density formula, one can tell that the EHM will possess a high Chern number for the presence of a large number of zero modes. Finally, we demonstrated the additivity of topological charges by exemplifying three various systems, including the combined Rashba-Dresselhaus system, bilayer graphene system, and heavy holes in quantum wells. The BG system (or heavy holes in QW ) possesses a single high topological charge which results from the additivity of separate unit charges. Moreover, the topological charge is independent of the separation scheme of zero-field points. The properties of Hall conductance are related to the topological charges of the systems, even though they are not quantized.
In conclusion, we have investigated the physical implications of zero modes of TS and TI systems as summarized in the last paragraph. Through the regularization method we proposed, we applied the charge density formula to the calculation of topological numbers for several various systems. This method facilitated the calculation of the topological numbers compared to the conventional methods [18,22]. Apparently, this regularization method can greatly simplify the mathematical process and avoid complicated integrations. Moreover, it offers a physical picture of zero modes playing in various topological systems. From the Teo-kane model and the case of Majorana zero modes bound to a flux line, localization of zero modes may indicate that they sitting on the defects or domain wall of the system. In a general sense, the edge can be treated like the domain wall between the system and the vacuum. This result is related to edge-bulk correspondence. Besides, it may be applied to other complicated topological systems.
Appendix A. Revisit of the residue theorem by the regularization method As mentioned in section 2, it is illuminating to revisit the residue theorem before investigating the calculation of topological numbers for the resemblance between the residue theorem and calculation of topological numbers. In this appendix, we would like to use the regularization method to revisit the results of the residue theorem. First, we demonstrate the usage of the regularization method by considering two complex functions with different topological properties. According to the well-known residue theorem, if the integrand of a line integral is analytic within and on a closed contour C except for a finite number of poles, its integral will be 2πi times the sum of the residues of the integrand at its poles within C. However, it is abstract to understand the sum of residues. As an intuitive alternative, we would like to convert the contour integral into the surface integral with Dirac delta functions originating from the regularization of the Pólya vector field associated with the integrand [21]. Here Dirac delta functions are nonzero only at the poles of the original integrand. The derivation of their integrals is usually much easier than calculating residues directly. When dealing with the topological properties of a vector field, it is often required to calculate an integral and poles in the integrand corresponding to zero-field points in the vector field.
The regularization method is described as follows. Let us express the complex value z as x+iy where x and y are real variables. A complex function F(z) can be constructed as u(x, y)+iv(x, y) where u(x, y) and v(x, y) are real functions. The vector field of the complex mapping F(z) is defined as F=[u, v] T . The Pólya vector field [21] associated with F(z) is defined as the complex conjugate of the mapping F(z) and is denoted by * = - is a complex function without poles, we can apply Green's theorem to the complex integral of F(z) over a closed contour C and obtain x y , and R is the enclosed region bounded by C. In general, F(z) has isolated finite poles within C and Green's theorem cannot be applied to the contour integral directly. In order to make equation (A.1) applicable, we will first replace F(z) by a smooth function F ε (z) characterized by a small positive number ε called the regulator. Then we can apply equation (A.1) to F ε (z) directly. After that, we take the limit as ε tends to zero on both sides of the equation, . Then the order of the limit and integral can be exchanged for generalized functions [27], and = . We then obtain Next, we make use of equation (A.2) to demonstrate the regularization method by two examples. For instance, when F(z)=1/z=(x−iy)/(x 2 +y 2 ), there is a singularity at z=0. To identify the nature of this singularity, we introduce an infinitesimal quantity ε 2 to its denominator, so that F ε =(x−iy)/(x 2 +y 2 +ε 2 ), which is equivalent to when the contour C includes the pole z=0 and it is null otherwise. The value of ∮ z z 1 d C is independent of the details of contour C as long as the same pole is encircled by the contour revealing the topological property of the integral. This is the same result as the one from the residue theorem. With regularization, we calculate the integral of delta functions instead of the contour integral on the left of equation (A.2).
In another case, we considered the function G(z)=i/z and its contour integral. It is apparent that there is a pole at z=0 and e = + + + e ( ) It is straightforward to generalize this regularization method to other complex functions with isolated poles.
The functions in the previous two examples, F(z) and G(z), represent curl-free and divergence-free vector fields on the complex plane, respectively. We have shown the regularization method is a useful and effective tool to derive the results of the residue theorem by replacing the residue calculation at the poles with twodimensional delta function integrals. The same approach could be adopted to regularize a vector field in the derivation of a topological charge density formula. It can transform the integration of topological numbers into the weighted sum over isolated zero modes. Similar to the residue theorem in complex variables, it simplifies the calculation procedure.
Appendix B. Signs of the Jacobian at points Γ and M for the EBHZ model Conventionally, the four zero-field points (0, 0), (0, π), (π, 0), and (π, π) are termed as points Γ, X, Y, and M, respectively. For the determination of the signs of Jacobian at points Γ and M of the EBHZ model, we present two different schemes to continuously deform its Hamiltonian by an insertion of a tuning parameter into the topological vector field. From the expression of the Jacobian given in section 5, we have J . Therefore, the signs of the Jacobian at points Γ and M (or v 1 and v 4 ) are uncertain and they will be remedied by introducing a tuning parameter into the EBHZ Hamiltonian, whose physical properties include signs of the Jacobian will be treated as the physical properties of the system in some parameter limit.

B.1. η-Hamiltonian approach
In the first approach, the EBHZ Hamiltonian is deformed into the η-Hamiltonian by inserting a tuning parameter η in the vector field components given by h v , their sign values of the Jacobian are −1, 1, 1, and −1 respectively. Now we focus on the sign of Jacobian at point Γ. We might simply take the limit sign of Jacobian at point Γ as that of the EBHZ model. It will be h -( ) sgn 1 2 in the limit of h  -1 , which is negative. However, for consistency, we need to investigate the limit from interval II . From equations (B.1) and (B.3), it is apparent that points Γ, X, Y, and M are solutions. In addition to these four points, there are two sets of zerofield points of interest in interval II of η, and each set contains two points. They can be found by recognizing the fact h =  k cos 1 and thus negative in interval II of η. The Jacobian at point Γ, η 2 −1, is positive in interval II. As η goes to 1 + , both of the first set of points will merge with point Γ. Thus the sum of Jacobian signs from these three points is −1, which equals that from another side of the limit, h  -1 . The consistency of the two-sided limit for the sign of Jacobian at point Γ is reached. The

B.2. λ-Hamiltonian approach
Another scheme-consistency issue may be raised for the sign determination of points Γ and M. This part of the appendix will address this issue by employing another scheme for deforming the EBHZ Hamiltonian into the λ-Hamiltonian.
In this scheme, a small parameter λ is introduced into h 1 , i.e. h , can be treated as a continuous deformation of H λ=0 . Namely, the vector field of the EBHZ model, h, can be treated as the limiting case of h λ as λ goes to zero. The zeros of h of the EBHZ model will also be these of h λ in the limit, l  0.
Now we turn to find the signs of the Jacobian for the λ-Hamiltonian. The null values of h 1 λ and l h 2 give algebraic relations of k x and k y of the zero-field points, expressed by: At first, we seek for the determination for the sign of Jacobian at point Γ in the EBHZ model, and focus on roots (or one root) of equations (B.5) and (B.6) near k y =0 and k x =0 for a small value of λ. Through the perturbative method, we found that there are two roots of equation (B.7) in the vicinity of k y =0. Their leading order expansion in the limit of l  0 are l k 2 . Therefore, the value of k x 1 is about π−λ/2, which is close to π. Thus, ( ) k k , x y 1 1 is close to point Y as λ approaches zero. It will determine the sign of Jacobian at point Y for the EBHZ model, which is already known as positive one.
On the other hand, substituting the latter solution k y 2 into equations (B.5) and (B.6) gets l k sin . Thus with substituted values of k x 2 and k y 2 , the Jacobian at ( ) k k , x y 2 2 is found as l -3 2 3 as l  0. Thus in both limits, l  + 0 and -0 , each limit sign value at ( ) k k , x y 2 2 exists and is −1, which gives the sign of Jacobian at point Γ in the EBHZ model. The same fashion will also work for that at the point M.
To conclude, both schemes have their own consistent two-sided limit in the determining Jacobian signs at points Γ and M, and reach the same conclusion for sign values as well.