Topological winding properties of spin edge states in Kane-Mele graphene model

# Topological winding properties of spin edge states in Kane-Mele graphene model

## Abstract

We study the spin edge states in the quantum spin-Hall (QSH) effect on a single-atomic layer graphene ribbon system with both intrinsic and Rashba spin-orbit couplings. The Harper equation for solving the energies of the spin edge states is derived. The results show that in the QSH phase, there are always two pairs of gapless spin-filtered edge states in the bulk energy gap, corresponding to two pairs of zero points of the Bloch function on the complex-energy Riemann surface (RS). The topological aspect of the QSH phase can be distinguished by the difference of the winding numbers of the spin edge states with different polarized directions cross the holes of the RS, which is equivalent to the topological invariance proposed by Kane and Mele [Phys. Rev. Lett. 95, 146802 (2005)].

###### pacs:
74.43.-f, 72.25.Hg, 85.75.-d
1

Since the discovery of the integer quantum Hall effect (IQHE) about three decades ago (1), topological invariants, which classify the electronic states, have been well accepted as a powerful tool for understanding the quantum many-body phases which have bulk energy gaps (2). Historically, the pioneer work was performed by Thouless, Kohmoto, Nightingale, and den Nijs (TKNN) (3), who recognized that the IQHE can be understood in terms of topological invariants known as Chern numbers (4), which are integrals of the -space Berry curvatures of the bulk states over the magnetic Brillouin zone. While IQHE finds its elegant connection through the adiabatic curvature with bulk topological invariants, Halperin (5) first stressed that the existence of the sample edges, which produces the current-carrying localized edge states in the Landau energy gap, is essential in the Laughlin’s gauge invariance argument (6). Hatsugai further developed a topological theory of the edge states (7), in which topological invariants are the winding numbers of the edge states on the complex-energy Riemann surface (RS).

Very recently, another striking topological quantum phenomenon, i.e., the quantum spin Hall effect (QSHE), was identified (8); (9) after long-distance efforts on metallic and conventionally insulating spin Hall effects (10); (11); (12), and soon after has been attracting extensive current interest (13); (14); (15); (16); (17); (18); (19); (20); (21); (22); (23); (24); (25); (26) due to its basic physics and its potential application in dissipationless spintronics. Unlike charge IQHE, whose presence fundamentally rely on the breaking of the time-reversal (-) symmetry via external magnetic field (1) or intrinsic magnetic gauge flux (27); (28), QSHE does not violate the -symmetry, which implies the absence of the non-zero Chern invariants in QSHE insulators. Then, it turns out that a -valued topological invariant could be associated with QSHE (9). This topology is, as one selective choice, characterized by whether the number of Kramers doublet localized at the edges in a strip geometry is even (nonzero) or odd. If even, the insulating phase is an ordinary Bloch insulator; otherwise, the insulating phase is a QSHE insulator.

To study QSHE and topological order, Kane and Mele introduced the model of graphene (8), which consists of two copies of Haldane’s model (27), one for spin-up electrons along some axis and one for spin-down electrons. -symmetry can be maintained if the intrinsic IQHE magnetic fields are opposite for the two spin components. To make the model more physical in that realistic mixing of the two spin components should emerges, a Rashba spin-orbit coupling term is further included (9), which now becomes the well-known Kane-Mele (KM) model. The KM model is the simplest possible model in that it has four spin-split bands, which is the minimum number required for the nontrivial phase to exist (26). For this reason, the KM model has received the most attention in the other studies (16); (17); (22); (29), in which the main focus is on the boundary phase twist and disorder effect on QSHE.

In this paper, we give a topological study of the spin edge states and its relation with QSHE by using the KM model. This study closely parallels with Hatsugai’s topological theory of edge-state IQHE in that we are seeking the winding numbers of the spin edge states on the complex-energy RS. We show that the spin-edge-state energy loops cross the holes of the RS, generating winding numbers and for spin-up and spin-down electrons, respectively. The quantized charge (spin) Hall conductance could be expressed as a summation (difference) of the spin-up and spin-down winding numbers. Thus we propose an edge-state topological invariant to distinguish a quantum spin-Hall (QSH) insulator from an ordinary insulator. If is zero, the insulating phase is an ordinary insulator; otherwise, the insulating phase is a QSH insulator. We stress that this classification between topological and ordinary insulating phases survives but is strongly modified by the weak mixing of the two spin components, which reflects the fact that although the exact quantization of the spin Hall conductance is destroyed by the spin-nonconserved perturbation, the QSH phase is still topologically distinct from the ordinary insulating phase.

Now we consider the tight-binding KM model of graphene (8); (9), which generalizes Haldane’s model (27) to include spin with -invariant spin-orbit interactions,

 H =t∑⟨ij⟩αc†iαcjα+iλSO∑⟨⟨ij⟩⟩αβνijc†iαszαβcjβ (1) +iλR∑⟨ij⟩αc†iα(s×^dij)zcjβ+λυ∑iξic†iαciα.

The first term is the nearest neighbor hopping term on the graphene (honeycomb) lattice, while the second term is the mirror symmetric spin-orbit interaction which involves spin-dependent second neighbor hopping. ==, where and are next-nearest neighbors, and and are unit vectors along the two bonds that connect to . is a Pauli matrix describing the electron’s spin. The third term, which will arise due to a perpendicular electric field or interaction with a substrate, is a nearest-neighbor Rashba coupling term, which explicitly violates the mirror symmetry. The fourth term is a staggered sublattice potential (=), which violates the symmetry under twofold rotations in the plane.

For the bulk graphene, this Hamiltonian can be written in the momentum space. For each the Bloch wave function is a four-component eigenvector of the Bloch Hamiltonian matrix ,

 H(k)=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝Z−λυY−iX0ia−Y+iX−Z+λυ−ia∗+00ia+−Z−λυY−iX−ia∗−0Y+iXZ+λυ⎞⎟ ⎟ ⎟ ⎟ ⎟⎠ (2)

with =, =, =, and =. Here =, =, and = represent the vectors from sublattice site (red circles in Fig. 1) to its three nearest sublattice sites (green circles in Fig. 1), respectively, and , , and represent the vectors between the nearest same sublattice sites. No doubt, the Hamiltonian (2) can also be written in terms of the SO(5) Clifford algebra, which is a very helpful technique to deduce -valued topological invariants (30).

To explore the edge topological invariant characterizing the QSH phase, we now turn to study the graphene ribbon with zigzag edges (see Fig. 1), which is periodic in the direction while it has two edges in the direction. In the following, we replace index with to denote the lattice sites, where label the unit cells and label the sublattice and in this cell. Since the graphene ribbon is periodic in the direction, we can use a momentum representation of the electron operator

 cmns,α=1√Lx∑keikX(mns)γns,α(k), (3)

where represents the coordinate of the site in the unit cell , and is the momentum along the direction. For simplification, let us first consider the spin-conserved case, i.e., the Rashba term in Hamiltonian (1) is set to be zero. Also, we let the staggered sublattice potential vanishes and only consider the first two terms in the Hamiltonian (1). Now let us consider the one-particle state =. Inserting it into the Schrödinger equation =, one can get the following eigenvalue equations for - and -sublattice sites:

 (ϵ+p2)ψnA↑ =p1ψnB↑+tψ(n−1)B↑ +p3[ψ(n+1)A↑+ψ(n−1)A↑], (4) (ϵ−p2)ψnA↓ =p1ΨnB↓+tΨ(n−1)B↓ −p3[Ψ(n+1)A↓+Ψ(n−1)A↓], (5) (ϵ−p2)ψnB↑ =p1ψnA↑+tψ(n+1)A↑ −p3[ψ(n+1)B↑+ψ(n−1)B↑], (6) (ϵ+p2)ψnB↓ =p1ψnA↓+tψ(n−1)A↓ +p3[ψ(n+1)B↓+ψ(n−1)B↓], (7)

where =, =, and . Eliminating the ()-sublattice sites, we obtain the difference equation for ()-sublattice sites

 f1ψn=f2[ψn+2+ψn−2]+f3[ψn+1+ψn−1], (8)

where =, =, = and () was replaced by . Eq. (8) is the so-called Harper equation (31). Note that since the spin is conserved, the spin-up and spin-down electrons satisfy the same Harper equation (8).

Now we represent Eq. (8) in the transfer matrix form. For this purpose we introduce a new wave function , which is a linear transformation of the original wave function ,

 φn=ψn±t∓ψn−1+ψn−2, (9)

where = with = and =. The new wave function can be written in the following transfer matrix form

 (φnφn−1)=~M(ϵ)(φn−1φn−2),

with

 ~M(ϵ)=(±t±−110). (10)

Although is a transformation of the original wave function , here we also can take and as the wave functions of two edges, which are linked by a reduced transfer matrix as follows:

 (φN+1φN)=M(ϵ)(φ1φ0), (11)

where

 M(ϵ)=[~M(ϵ)]N=(M11M12M21M22). (12)

All kinds of solutions of Eq. (12) are obtained by different choices of and . Now we study the energy spectrum of the graphene ribbon with special attention to the edge states. The general boundary condition is

 φN+1=φ0=0. (13)

With Eqs. (11), (12) and the boundary condition (13), one can easily get that the solutions satisfy

 M21(ϵ)=0, (14)

and

 φN=M11(ϵ)φ1. (15)

If we use a usual normalized wave function, the state is localized at the edges as

 Unknown environment '% (16)

Because the analytical derivation is very difficult, we now turn to start a numerical calculation from Eq. (1) with and setting to be zero. Figure 2(a) shows the energy bands for the graphene with zigzag edges with the intrinsic spin-orbit coupling strength =, same as that in Ref. (8). From Fig. 2(a) one can see that there are two branches of spin-degenerate dispersed energy bands with two pairs of spin states lying in the energy gap. These two pairs of spin states cross at the -invariant point and are localized at the edges of the graphene ribbon. To show the spin localization features of these gapless spin states, we plot in Figs. 2(b) and 2(c) the spin (-component) distribution along the graphene strip for the two pairs of edge states at 0.99 and 1.01, respectively. At the left side of the crossing point [Fig. 2(b)], the spin-up (-down) edge state with lower energy is localized at the down (up) edge of the graphene ribbon, while the spin-up (-down) edge state with upper energy is localized at the up (down) edge of the graphene ribbon. On the other hand, at the right side of the crossing point [Fig. 2(c)], the spin-up (-down) edge state with lower energy is localized at the up (down) edge of the graphene ribbon, while the spin-up (-down) edge state with upper energy is localized at the down (up) edge of the graphene ribbon. Thus, if we trace the spin-up flow of the edge states when varying the momentum in one period, it can be found that the electrons with different spin propagate in opposite directions and that the electronic state of graphene is a QSH state. The spin Hall conductivity has been shown as = in units of by different methods.

Before introducing the winding number of the spin edge states on the complex energy RS, let us try to interpret the above two numerical characteristics. The first one is that there are always two energy-degenerate edge states appearing localized at opposite edges. To explain this issue, we suppose the spin edge states to be exponentially localized on the boundary with the following ansatz (33); (34)

 ψn=λnψ, (17)

where is a complex number. Inserting Eq. (17) into the Harper equation (8), one can easily get the complex number satisfying the following equation

 Missing or unrecognized delimiter for \left (18)

That means

 λ+λ−1=±t±. (19)

It is trivial to see from Eqs. (18) and (19) that if is a solution, then is also to be. That means if there is an edge state localized near one boundary, at the same time there is another edge state localized near opposite boundary. Moreover, note that the conservation of the spin , i.e., the spin-up and -down electrons satisfy the same Harper equation (8), one easily find that if is a solution describing an edge state with spin-up localized near one boundary, then is another solution describing an edge state with spin-down localized near opposite boundary. That is the intrinsic reason for the second feature that there are always two edge states with opposite spins localized near opposite edges.

Now we show that similar to charge IQHE, QSHE of the graphene ribbon is associated with winding number of the spin edge states on a complex energy RS. Within the topological edge theory (7), this number is, as well as the Chern number, a well-defined topological quantity. In order to well understand the winding numbers of the spin edge states in QSHE, let us first consider that of the edge states in charge IQHE. First we ignore the open condition and consider the bulk Bloch function at sites with -coordinate of . For Bloch function,  and  compose an eigenvector of  with the eigenvalue ,

 M(ϵ)⎛⎜⎝φ(b)1φ(b)0⎞⎟⎠=ρ(ϵ)⎛⎜⎝φ(b)1φ(b)0⎞⎟⎠. (20)

In order to discuss the wave function of the edge state, we extend the energy to a complex energy. In the following, we use a complex variable  instead of real energy. From Eq. (20) we get

 ρ(z)=12[Δ(z)−√Δ2(z)−4] (21)

and

 φN=−M11(z)+M22(z)−ω−M11(z)+M22(z)+ωM21(z)φ1, (22)

where =Tr and =. Clearly,

 detM(ϵ)=1 (23)

since =. From Eq. (22), one can find that the analytic structure of the wave function is determined by the algebraic function =. The RS of = on the complex energy plane can be built by the conglutination between different analytic brunches. The close complex energy plane can be obtained from the open complex energy plane through spheral pole mapping. Now let us discuss the analytic structure of on the open complex plane. If the system has  energy bands, i.e.,

 ϵ∈[λ1,λ2],...,[λ2j−1,λ2j],...,[λ2q−1,λ2q], (24)

where  denote energies of the band edges and , , then  can be factorized as

 ω=√Δ2(z)−4= ⎷2q∏j=1(z−λj). (25)

In the present graphene system, there are two energy bands, so =. The two single-valued analytic branches are defined on the plane with two secants down or up the bank of which we choose corresponding complex angles, which is shown in Fig. 3.

In order to ensure (=,,=) [Here, is the energy of the edge state in the gap ], we divide the two single-value analytic branches based on the parity of . Let us consider the case that the energy lies in the band [see Fig. 3(a)]. The left side of this energy band is the th gap, and the right side is the th gap. When the energy  moves from the th band to the th gap (the th gap) along an arbitrary path (), only the singularities and  have contributions to the variance of the principal value of argument of . On the up bank of the secant, we distinguish two branches  and as the following:

For is even, if we set = and =, which corresponds to () when moves along (), then the branch is defined as a complex plane with secants. If we set = and =, which corresponds to () when moves along (), then the branch is defined as a complex plane with secants. For is odd, the definitions of and are reverse to those for is even. So, if  lies in the th gap from the lowest one on the real axis,

 α(−1)jω≥0, z (real) on Rα (α=+,−), (26)

and at energies of the edge states we have

 ω(μj)=α(−1)j∣∣M11(μj)−M22(μj)∣∣, (27)

where , =. In addition, one can easily obtain

 Δ(ϵ){≤−2 for j odd≥2 for j even, (28)

where the energy  (on ) is in the th gap.

When the branches and  on the open complex energy plane are mapped to the close complex energy plane through spheral-pole-projection, one can get two single-value analytic spherical surfaces. The RS is obtained by gluing the two spherical surfaces at these branch cuts. Make sure the banks face the banks of other sphere. (see Figs. 3(b)). Note that there are two real axes after gluing. In our model, the genus of the RS is = for spin-up (or spin-down) electrons, which is the number of energy gaps. The wave function is defined on the = RS . The branch of the Bloch function is specified as , which we have discussed above. With Eqs. (28), (27) and (22), and using the fact when  for and  for , one can obtain that when the zero point is on the upper sheet of RS (), the edge state is localized at the down edge; when the zero point is on the lower sheet of RS (), the edge state is localized at the up edge.

In Fig. 4, on the RS of our system the energy gap corresponds to the loop around the hole of the and the energy bands correspond to the closed paths vertical to the energy gap loop on the . The Bloch function is defined on this surface. For the fixed there is always = zero point at the down-edge-state energy . Since there are two real axes on the , correspondingly, there is = zero point at the up-edge-state energy .

The above considerations are for the fixed . Now, let us consider a family of the RSs parametrized by changing in one of its periods. can be modified by this change. However, all the RSs with different are topologically equivalent if there are stable energy gaps in the 2D energy spectrum.

On the genus = RS, the first homotopy group is generated by = generators,  and , =. The intersection number of the curves (including directions) (7) is given by (see Fig. 5)

 I(αi,βj)=δij. (29)

Any curves on the RS are spanned homotopically by  and . When   moves  times around the th hole with some integer , one has

 C(μj)≈βpj, (30)

which means

 I(αi,C(μj))=tδij. (31)

When the Fermi energy  of the 2D system lies in the th energy gap, the charge Hall conductance is given by the winding number of the edge state, which is given by the number of intersections () between the canonical loop  on the RS and the trace of , i.e.,

 σj,edgexy=−e2hI(C(μj)). (32)

Similarly, the above expression can be obtained by the Byers-Yang (32) and Laughlin-Halperin (6); (5) gauge arguments. The system with the periodic boundary condition in the direction and open condition in the direction can be considered as a cylinder. By Laughlin gauge invariance argument, the vector potential  has to have the form =( is an integer). When the flux  threading the cylinder is adiabatically turn on from = to =  with = a flux quantum, = accords with the gauge argument, therefore, () maps the system back to itself. Basing on the single-electron assumption, When the Fermi energy lies in the th energy gap, there are states (electrons) transferring from the down edge (=) to the up edge (=) in net. The energy change during the adiabatic process is =, where  is a voltage in the direction. This gives the charge Hall current  as follows

 Ix=cΔEΔΦ=σxyVy.

Then we get an expression for  as Eq. (32).

The above analysis is for the winding number of the edge states in IQHE and it can be easily generalized to the winding number of the spin edge states in QSHE when the spin degree of freedom is considered. In the latter case, the winding number of the spin edge state is given by the number of intersections (, , for spin-up and -down, respectively) between the canonical loop  on the RS and the trace of . From Fig. 4 one can observe that  moves one time across the hole (=), which means   and =. Considering the winding direction (see Fig. 5), one obtain = for spin-up electrons while = for spin-down ones. Thus, the charge Hall conductance is given by the summation over two spin channels,

 σ(c)edgexy=(I↑+I↓)e2h≡Ice2h=0, (33)

while the spin Hall conductance is given by the difference between them,

 σ(s)edgexy=(I↑−I↓)ℏ2e2h≡Ise4π=2e4π. (34)

When the finite on-site energy is considered while the Rashba coupling is kept as zero, the two degenerate energy bands will split. For the bulk graphene (2), one can easily obtain the four dispersion energy bands,

 ϵk=±√X2+Y2+(Z±λv)2, (35)

between which there is an energy gap with magnitude =. We plot in Fig. 6 the energy dispersion with different and . The condition provides a finite bulk energy gap in the graphene. However, when the Fermi energy lies in the gap, the graphene’s phases are different with different topology: for , the system is an ordinary insulator and for , the system is a topological insulator with QSH phase. Kane and Mele (9) have proposed a index describing the QSH phase, which is determined by counting the number of pairs complex zeros of =Pf with the wave functions corresponding to the bulk Hamiltonian (2) and time reversal operator. To understand these different phases, Kane and Mele also studied the edge states of the graphene ribbon with zigzag edges. However, the winding properties of the edge states in the KM model has not been previously considered, which just is the special focus of our present study.

Considering nonzero , now let us investigate the topological winding numbers with zigzag edges. In this case, the Harper equation for site turns to have the same form as Eq. (8), provided that the in the coefficients and is replaced by for spin-up and -down channels, respectively. Similarly, one can introduce a new wave function, as a linear transformation of the original wave functions (9), satisfying the matrix form with the transform matrix (10). The in the coefficients , , and is now also replaced by for different spins. In evidence, the edge-state energies become non-degenerate for the finite on-site energy . Figure 7 shows the energy spectrum of the “zigzag” graphene ribbon in a QSH phase with = and =, which satisfy . From Fig. 7, one can clearly see that the gapless edge states are localized near the system boundary with different spins. The corresponding complex energy RS of the Bloch wave functions is topologically equivalent to Fig. 3. And no doubt this QSH phase can be described by the difference of the winding numbers of the edge states with different spins, i.e., ==, in units of . On the other hand, one can obtain the charge Hall conductivity is zero since =. As a comparison, we draw in Fig. 8 the energy spectrum of an ordinary insulator phase with = and =. One can find that in this case (see Fig. 9), there are no gapless edge states connecting the two energy bands. The traces of the edge states in the corresponding complex energy RS are sunk in the bulk bands, i.e., ==, which means that the system in this case is an ordinary insulator without QSHE. We thus conclude that while the previously studied topological winding number index (=+) (7); (35); (36) has been verified to describe the IQH insulator, our present studied topological winding number index (=) can be used to characterize the QSH insulator.

Finally let us turn to consider the most complex case of being nonconserved, i.e., the Rashba coupling , as well as the on-site energy , is finite in the KM graphene model (1). After a straightforward derivation, we obtain the following couple eigenvalue equations for :

 (ϵ±p2±)ΨnA± =p1ΨnB±+t1Ψ(n−1)B± ±p3[Ψ(n+1)A±+Ψ(n−1)A±] +ip∓ΨnB∓−iλRΨ(n−1)B∓, (36)
 (ϵ∓p2±)ΨnB± =p1ΨnA±+t1Ψ(n+1)A± ∓p3[Ψ(n+1)B±+Ψ(n−1)B±] −ip±ΨnA∓+iλRΨ(n+1)A∓, (37)

where =, =, and “” in label the spin channels. Although the Harper equation and the corresponding transform matrix are very difficult to derive in this case, we can also distinguish the QSH phase from the ordinary insulator phase by the topological winding index introduced in the present paper directly from the energy spectrum. As an example, we reinvestigate the QSH phase and the ordinary insulator phase with the Rashba coupling being a finite value =. The other parameters are chosen be the same as thosed used in Fig. 7. The calculated energy spectrum (not shown here) turns to be essentially the same as that depicted in Fig. 7. As a result, the corresponding RS of the Bloch wave functions is topologically equivalent to Fig. 3. This means that in spite of finite Rashba term, the system is still in the QSH phase with the topological index =. Therefore, although in this case the spins are nonconserved (see Table I) and the exact quantization of the spin Hall conductance is destroyed by the spin-nonconserved perturbation, the QSH phase of the graphene system is still topologically distinct from the ordinary insulating phase once provided that the Rashba term does not change the system’s topological properties. This conclusion keeps consistent with Ref. (9). If the Rashba spin-orbit coupling turns much strong, the concept of the “spin edge states” becomes faint and the winding number of the spin edge states becomes meaningless.

Before ending the present paper, we would like to stress that the topological index describing the QSH phase can also be interpreted in terms of Laughlin and Halperin’s arguments (6); (5). In fact, in the QSH system with topological integer =, there are spin-up electrons transferred from one edge to the other when a unit magnetic flux is adiabatically through a cylindrical system by branches of gapless edge states. And at the same time there are spin-down electrons to be transferred with the opposite direction (see Fig. 10).

In summary, by theoretically studying within the KM Hamiltonian a graphene strip with zigzag edges, we have presented an alternative topological index to characterize the QSH phase. The topological index describing QSH phases has been defined as the difference of the winding numbers of the spin-resolved edge states crossing the holes of the complex-energy RS. Based on this topological index, we have discussed different phases by modulating different parameters in KM model, which agree well with the previous studies in terms of the conventional topological invariance.

###### Acknowledgements.
This work was supported by NSFC under Grants Nos. 10604010, 10534030, and 60776063, and the National Basic Research Program of China (973 Program) under Grant No. 2009CB929103.

### References

1. K. v. Klitzing, G. Dorda, and M. Peper, Phys. Rev. Lett. 45, 494 (1980).
2. X. G. Wen, Phys. Rev. B 40, 7387 (1989).
3. D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
4. D. J. Thouless, Topological Quantum Numbers in Nonrelatisvistic Physics (World Scientific, Singapore, 1998).
5. B. I. Halperin, Phys. Rev. B 25, 2185 (1982).
6. R. B. Laughlin, Phys. Rev. B 23, 5632 (1981).
7. Y. Hatsugai, Phys. Rev. Lett. 71, 3697 (1993); Phys. Rev. B 48, 11851 (1993).
8. C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
9. C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
10. S. Murakami, N. Nagaosa, and S.C. Zhang, Science 301, 1348 (2003).
11. J. Sinova, D. Culcer, Q. Niu, N. A. Sinitsyn, T. Jungwirth, and A. H. MacDonald, Phys. Rev. Lett. 92, 126603 (2004).
12. S. Murakami, N. Nagaosa, and S.-C. Zhang, Phys. Rev. Lett. 93, 156804 (2004).
13. M. Onoda and N. Nagaosa, Phys. Rev. Lett. 95, 106601 (2005).
14. B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
15. X.-L. Qi, Y.-S. Wu, and S.-C. Zhang, Phys. Rev. B 74, 085308 (2006).
16. D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane, Phys. Rev. Lett. 96, 036808 (2006).
17. T. Fukui and Y. Hatsugai, Phys. Rev. B 75, 121403(R) (2007).
18. L. Fu and C. L. Kane, Phys. Rev. B 74, 195312 (2006).
19. L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
20. L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
21. S. Murakami, Phys. Rev. Lett. 97, 236805 (2006).
22. M. Onoda, Y. Avishai, and N. Nagaosa, Phys. Rev. Lett. 98, 076802 (2007).
23. C. Wu, B. A. Bernevig, and S.-C. Zhang, Phys. Rev. Lett. 96, 106401 (2006).
24. C. Xu and J. E. Moore, Phys. Rev. B 73, 045322 (2006).
25. B. A. Bernevig, T. L. Hughes, S.-C. Zhang, Science 314, 1757 (2006).
26. J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (R) (2007).
27. F. D. M. Haldane, Phys. Rev. Lett. 61, 2015 (1988).
28. K. Ohgushi, S. Murakami, and N. Nagaosa, Phys. Rev. B 62, R6065 (2000).
29. A. M. Essin and J. E. Moore, Phys. Rev. B 76, 165307 (2007).
30. S. Murakami, N. Nagaosa, and S. C. Zhang, Science 301, 1348 (2003); Phys. Rev. B 69, 235206 (2004).
31. P. G. Harper, Proc. Phys. Soc. London Sect. A 68, 874 (1955).
32. N. Byers and C. N. Yang, Phys. Rev. Lett. 7, 46 (1961).
33. M. König, H. Buhmann, L. W. Molekamp, T. L. Hughes, C.-X. Liu, X.-L. Qi, and S.-C. Zhang, arXiv/0801.0901.
34. M. Creutz and I. Horvath, Phys. Rev. D 50, 2297 (1994); M. Creutz, Rev. Mod. Phys., 73, 119 (2001).
35. Z. Wang and P. Zhang, Phys. Rev. B 77, 125119 (2008).
36. N. Hao, P. Zhang, Z. Wang, W. Zhang, and Y. Wang, Phys. Rev. B 78, 075438 (2008).
101060