Absence of equilibrium edge currents in theoretical models of topological insulators

Our purpose is to investigate the equilibrium currents in TIs produced by the entire band structure, including both the edge state Dirac cone and the bulk bands, by means of minimal lattice models that can take into account these features. For this purpose we start from the low energy effective Hamiltonian of 2D and 3D TIs described by a linear Dirac model

$\displaystyle H=\sum_{\ell=1}^{D}A\,k_{\ell}\Gamma^{\ell}+M\,\Gamma^{0},$

(1)

where the $\Gamma^{\ell}$ matrices satisfy the Clifford algebra $\left\{\Gamma^{\ell},\Gamma^{m}\right\}=2\delta_{\ell m}$, whose precise forms depend on the symmetry classesSchnyder et al. (2008); Ryu et al. (2010); Chiu et al. (2016); Chen and Schnyder (2019). The kinetic terms $k_{\ell}$ signifies the usual linear band crossing, and $M$ is the mass term that controls the topology. We proceed to regularize the Hamiltonian to cover the entire Brillouin zone (BZ) by

	$\displaystyle A\,k_{\ell}\rightarrow A\sin k_{\ell}$
	$\displaystyle M\rightarrow M+B\sum_{\ell=1}^{D}k_{\ell}^{2}\rightarrow M+4B-2B\sum_{\ell=1}^{D}\cos k_{\ell},$		(2)

where the $k^{2}$ term is important to regularize the model on a latticeAraújo et al. (2019). The basis of the Hamiltonian depends on the symmetry classes, which may be electrons in different orbitals with or without spinsSchnyder et al. (2008); Ryu et al. (2010); Chiu et al. (2016). We then Fourier transform the second quantized Hamiltonian into real space to construct the 2D square lattice or 3D cubic lattice Hamiltonian. For instance, if the basis are electron creation and annihilation operators $c_{iI\sigma}^{{\dagger}}$ and $c_{iI\sigma}$ with momentum $k$, orbital $I$, and spin $\sigma$, then the Fourier transform leads to

	$\displaystyle\sum_{k}\cos k_{\ell}c_{kI\sigma}^{{\dagger}}c_{kJ\sigma}\rightarrow\frac{1}{2}\sum_{i\sigma}\left\{c_{iI\sigma}^{{\dagger}}c_{i+\ell J\sigma}+c_{i+\ell I\sigma}^{{\dagger}}c_{iJ\sigma}\right\},$
	$\displaystyle\sum_{k}i\sin k_{\ell}c_{kI\sigma}^{{\dagger}}c_{kJ\sigma}\rightarrow\frac{1}{2}\sum_{i\sigma}\left\{c_{iI\sigma}^{{\dagger}}c_{i+\ell J\sigma}-c_{i+\ell I\sigma}^{{\dagger}}c_{iJ\sigma}\right\},$

and so follows the lattice Hamiltonian defined on sites $i$.

The 2D class A models are the simplest ones to demonstrate the absence of edge current in the regularized lattice models introduced in Sec. II A. We use the prototype Chern insulator to demonstrate our statements in this class. The model is described by the spinless basis $\psi=\left(c_{{\bf k}s},\;c_{{\bf k}p}\right)^{T}$ and HamiltonianQi and Zhang (2011); Bernevig and Hughes (2013)

	$\displaystyle H({\bf k})$	$\displaystyle=$	$\displaystyle A\sin k_{x}\sigma^{x}+A\sin k_{y}\sigma^{y}$		(4)
		$\displaystyle+$	$\displaystyle\left(M+4B-2B\cos k_{x}-2B\cos k_{y}\right)\sigma^{z},$

which leads to the following lattice model

	$\displaystyle H=\sum_{i}t\left\{-ic_{is}^{{\dagger}}c_{i+ap}+ic_{i+as}^{{\dagger}}c_{ip}+h.c.\right\}$
	$\displaystyle+\sum_{i}t\left\{-c_{is}^{{\dagger}}c_{i+bp}+c_{i+bs}^{{\dagger}}c_{ip}+h.c.\right\}$
	$\displaystyle+\sum_{i\delta}t^{\prime}\left\{-c_{is}^{{\dagger}}c_{i+\delta s}+c_{ip}^{{\dagger}}c_{i+\delta p}+h.c.\right\}$
	$\displaystyle+\sum_{i}\left(M+4t^{\prime}\right)\left\{c_{is}^{{\dagger}}c_{is}-c_{ip}^{{\dagger}}c_{ip}\right\}-\sum_{iI}\mu\,c_{iI}^{{\dagger}}c_{iI},$		(5)

where $t=A/2$, $t^{\prime}=B$, $I=\left\{s,p\right\}$ are the orbitals, and $\delta=\left\{a,b\right\}$ are the lattice constants along planar directions. Although the chemical potential $\mu$ is frequently ignored because it does not alter the topology of the system, we will elaborate its importance for the existence of the edge current. We impose periodic boundary condition (PBC) in the ${\hat{\bf x}}$ direction and open boundary condition (OBC) in the ${\hat{\bf y}}$ direction, such that the edge currents are localized at the two edges $y=1$ and $y=N_{y}$. The band structure at parameters $t=t^{\prime}=-M=1$ with a finite width $N_{y}=8$ is shown in Fig. 1 (a). For each eigenstate $|k_{x},n_{y}\rangle$ we calculate the weight of its wave function close to the $y=1$ boundary minus that close to the $y=N_{y}$ boundary

$\displaystyle\tilde{n}_{k_{x},n_{y}}=\left(\sum_{1\leq y\leq N_{y}/2}-\sum_{N_{y}/2+1\leq y\leq N_{y}}\right)|\psi_{k_{x},n_{y}}(y)|^{2}.$

(6)

The color in the band structure in Fig. 1 (a) indicates a positive (red) or negative (green) $\tilde{n}_{k_{x},n_{y}}$, i.e., whether the $|k_{x},n_{y}\rangle$ is more localized in the $y=1$ or the $y=N_{y}$ boundary. The result clearly indicates that the branch of the Dirac cone with positive group velocity is localized at the $y=1$ edge (red), whereas the branch of negative group velocity is localized at the $y=N_{y}$ edge (green). In addition, the valence bands at $k_{x}>0$ that have negative group velocities are also more localized at the $y=1$ edge, indicating that they contribute to a current against that of the edge states, and likewise the valence bands at $k_{x}<0$ that have positive group velocity are more localized at $y=N_{y}$ edge and also contribute to a current against that of the edge states.

We proceed to quantify the contribution from the edge states and that from the valence bands by calculating the expectation value of the local current operator, constructed from the equation of motion of the density operator $n_{i}=\sum_{I}c_{iI}^{{\dagger}}c_{iI}$ written in the form of continuity equations

$\displaystyle\dot{n}_{i}=\frac{i}{\hbar}\left[H,n_{i}\right]=-{\bm{\nabla}}\cdot{\bf J}_{i}^{0}=-\frac{1}{a}\sum_{\delta}\left(J_{i,i+\delta}^{0}+J_{i,i-\delta}^{0}\right),\;\;\;$

(7)

which defines the charge current operator $J_{i,i+\delta}^{0}$ running from site $i$ to $i+\delta$, and that run from $i$ to $i-\delta$. We define the currents flowing in the positive bonds as $J_{x}^{0}\equiv J_{i,i+a}^{0}$ and $J_{y}^{0}\equiv J_{i,i+b}^{0}$, which have the expressions

	$\displaystyle J_{x}^{0}$	$\displaystyle=$	$\displaystyle\frac{a}{\hbar}t\left\{c_{is}^{{\dagger}}c_{i+ap}+c_{i+ap}^{{\dagger}}c_{is}+c_{ip}^{{\dagger}}c_{i+as}+c_{i+as}^{{\dagger}}c_{ip}\right\}$
		$\displaystyle+$	$\displaystyle\frac{a}{\hbar}it^{\prime}\left\{-c_{is}^{{\dagger}}c_{i+as}+c_{i+as}^{{\dagger}}c_{is}+c_{ip}^{{\dagger}}c_{i+ap}-c_{i+ap}^{{\dagger}}c_{ip}\right\}$
	$\displaystyle J_{y}^{0}$	$\displaystyle=$	$\displaystyle\frac{b}{\hbar}it\left\{-c_{is}^{{\dagger}}c_{i+bp}+c_{i+bp}^{{\dagger}}c_{is}+c_{ip}^{{\dagger}}c_{i+bs}-c_{i+bs}^{{\dagger}}c_{ip}\right\}$
		$\displaystyle+$	$\displaystyle\frac{b}{\hbar}it^{\prime}\left\{-c_{is}^{{\dagger}}c_{i+bs}+c_{i+bs}^{{\dagger}}c_{is}+c_{ip}^{{\dagger}}c_{i+bp}-c_{i+bp}^{{\dagger}}c_{ip}\right\}.$

To separate the contribution from the Dirac cone and that from the bulk bands, in calculating the expectation value we introduce an energy cutoff $E_{cut}$ within which we sum the eigenstates

$\displaystyle\langle J_{x}^{0}\rangle_{cut}=\sum_{n}\langle n|J_{x}^{0}|n\rangle f(E_{n})\theta(E_{cut}-|E_{n}|),$

(9)

where $E_{n}$ and $|n\rangle$ are eigenenergy and eigenstate of the lattice model, and $f(E_{n})=\left(e^{E_{n}/k_{B}T}+1\right)^{-1}$ is the Fermi distribution, $\theta(E_{cut}-|E_{n}|)$ is the step function that selects the energy window within which the states are included, and we choose $k_{B}T=0.03$ that corresponds to the room temperature throughout the article, assuming the energy unit $t=$ eV.

The results shown in Fig. 1 (a) indicates that choosing an energy cutoff smaller than the bulk gap $E_{cut}<|M|=1$, which corresponds to counting the contribution from the Dirac cone alone, indeed yields a finite current localized at the two edges. As increasing $E_{cut}$, the spatial profile of the current starts to alter, and eventually the current vanishes everywhere at a large $E_{cut}=8$ that includes contributions from all the valence bands, proving that the two contributions cancel out exactly. This conclusion holds regardless the temperature, distance between the two edges $N_{y}$, and band parameters $\left\{t,t^{\prime},M\right\}$ within our regulated lattice model, provided the system remains in the topologically nontrivial phase $M<0$. Thus unless one has certain experimental probe that can discern the edge current contributed only within an energy window near the chemical potential, one cannot resolve the contribution from the edge states alone and the total current should be zero. Remarkably, the edge current profile at a specific choice of energy cutoff $E_{cut}$ (each line in Fig. 1 (b)) remains the same at any temperature, even up to a unrealistically high temperature $k_{B}T=1\sim 10000$K. The difference at different temperatures remains less than $10^{-10}$ in our numerical calculation.

It is also immediately clear that if the Dirac point is shifted away from the zero energy by adjusting chemical potential $\mu$, then a finite edge current occurs, since some states below the Dirac cone are omitted if the chemical potential lies below the Dirac point, and some more states are included if the chemical potential lies above the Dirac cone, hence the cancellation from the valence bands is not exact. This argument also implies that the edge current must change sign as chemical potential sweeps through the Dirac point $\mu=0$, as confirmed by the total current close to the $y=1$ edge $\langle J_{x}^{0}\rangle_{h}=\sum_{1\leq y\leq N_{y}/2}\langle J_{x}^{0}\rangle$ shown in Fig. 1 (c).

Moreover, Fig. 1 (c) also indicates that the edge current as a function of the chemical potential $\mu$ is stepwise quantized within the bulk gap $|\mu|<|M|=1$ (the steps are smeared by temperature). This behavior is due to the fact that the Dirac cone is geometrically quantized by the small size of the strip, as indicated by the grey circles inside the dashed box in Fig. 1 (a). As a result, every time $\mu$ sweeps through one of these quantized edge states, the edge current increases by one unit. The number of plateaux $N_{p}$ and the edge current quantum $\Delta J_{x}^{0}$ can be understood in the following manner. The width of the dashed box in Fig. 1 (a) is the insulating gap divided by the group velocity of the Dirac cone $k_{p}=M/\hbar v_{F}$, and the spacing between states due to geometric quantization is $\Delta k=2\pi/L_{x}=2\pi/N_{x}a$, so the number of circles that can fit into the dashed box in Fig. 1 (a) is

$\displaystyle N_{p}=\frac{k_{p}}{\Delta k}=\frac{MN_{x}a}{2\pi\hbar v_{F}}.$

(10)

On the other hand, because each quantized edge state spreads out the whole edge of length $L_{x}=N_{x}a$ and decays into the bulk with length $\xi$, one may approximate the 2D density of the edge state by $n_{2D}=1/N_{x}a\xi$. Consequently, each quantized edge state contributes to a current density $\Delta j_{2D}^{0}=en_{2D}v_{F}$. The edge current quantum is this value multiplied by the cross section $\xi$ of the current flow

$\displaystyle\Delta J_{x}^{0}=\Delta j_{2D}^{0}\xi=\frac{ev_{F}}{N_{x}a}.$

(11)

Equations (10) and (11) well explain the number and height of the plateaux in Fig. 1 (c). Finally, we remark that the finite edge current only occurs in the topologically nontrivial phase $M<0$ in this model. In the topologically trivial phase $M>0$ there is no edge state, and neither do the valence bands contribute to a current. Finally, an experimental setup for measuring this quantized edge current is proposed in Fig. 1 (d), which will be discussed in the next section for the more realistic quantum dot geometry.

The results in Sec. II.2 prompts us to explore the effects of weak impurities and edge confining potentials, since it helps to understand the situation when the chemical potential is altered locally. For this purpose, we consider random impurities of density $n_{imp}=N_{imp}/N_{x}N_{y}$ by adding the term

$\displaystyle H_{imp}=\sum_{i\in imp,I}U_{imp}c_{iI}^{{\dagger}}c_{iI},$

(12)

into our Hamiltonian in Eq. (5), where $i\in imp$ denotes impurity sites. In Fig. 2 we present the simulation in an open island of Chern insulator quantum dot (OBC imposed in both ${\hat{\bf x}}$ and ${\hat{\bf y}}$ directions), which confirms that a local current is induced around the impurities. The local current forms a vortex circulating the impurity, and the direction of the vorticity depends on the sign of $U_{imp}$. The vortices interfere if the impurities are getting too close to each other. Comparing small $U_{imp}=\pm 0.3$ and large $U_{imp}=\pm 0.6$, one sees that the magnitude of the local current is enhanced at large impurity potentials, at least within this weak impurity strength we investigate. For the $U_{imp}=\pm 0.6$ case in Fig. 2 (b) we also include a confining potential at the edge sites $i\in edge$

$\displaystyle H_{con}=\sum_{i\in edge,I}V_{con}c_{iI}^{{\dagger}}c_{iI},$

(13)

with $V_{con}=0.6$ chosen the same as the impurity potential. One sees that the confining potential indeed induces an edge current circulating around the edge, since it locally shifts the chemical potential, and the magnitude of the current (size of the arrows in Fig. 2 (b)) is generally larger than that induced by point-like impurities. Interestingly, Fig. 2 (b) also reveals that the current on the boundary sites (e.g., $y=1$) and the sites next to the boundary (e.g., $y=2$) actually have opposite directions of flow, meaning that the edge current induced by the edge confining potential is a laminar flow whose direction of flow depends on the distance to the edge. Note that this confining potential induced edge current bears a striking similarity with the edge current in the quantum Hall effect of 2D electron gas, which is also generated by edge confining potentialHalperin (1982); Smrcka (1984); MacDonald and Středa (1984); Streda et al. (1987); Büttiker (1988); Chklovskii et al. (1992), although there are no Landau levels in our problem.

An experimentally relevant system that belongs to 2D class A is the thin film of magnetic TIs of few quintuple layers that manifests quantum anomalous Hall effect (QAHE)Chang et al. (2013); Liu et al. (2016); He et al. (2018). Although the equilibrium edge current will not show up in the nonequilibrium transport measurements, the magnetic field it produces should in principle be measurable. Assuming a Chern insulator quantum dot of dimension $N_{x}a\sim 100$nm can be realized by these systemsFerreira and Loss (2013), the edge current in the quantum dot is quantized as a function of gate voltage $\mu$ similar to that shown in Fig. 1 (c). A typical Fermi velocity $v_{F}\sim 10^{5}$m/s gives the edge current quantum $\Delta J_{x}^{0}\sim 10^{-7}$A. Ampere’s law then gives an increase of magnetic field at the center $\Delta B=\mu_{0}\Delta J_{x}^{0}/2\pi N_{x}a\sim 10^{-7}$T, and consequently an increase of magnetic flux $\Phi_{B}\sim\pi(N_{x}a)^{2}B\sim 10^{-13}$Wb$\sim 100\Phi_{B}^{0}$ per edge current quantum, where $\Phi_{B}^{0}$ is the magnetic flux quantum. Thus the magnetic flux through the quantum dot as a function of gate voltage is also quantized, which should be measurable by magnetic flux detectors such as the superconducting quantum interference device (SQUID), as indicated schematically in Fig. 1 (d).

We now examine a Chern insulator quantum dot that is partially gated, and show that the edge current created in the gated area can extend into the ungated area over length scale $\xi$. To elaborate this statement, in Fig. 3 we show the local current pattern in a $N_{x}\times N_{y}=32\times 8$ quantum dot with OBC imposed in both ${\hat{\bf x}}$ and ${\hat{\bf y}}$ directions. A quarter of the area is gated to have $\mu=0.1$, whereas the ungated area is kept at $\mu=0$. The result indicates that the boundary between the gated and ungated areas serves as an edge, and an edge current is created in both sides of the edge even though only one side has a finite chemical potential $\mu$. The circulating pattern of the current decays into the ungated area over length $\xi$, as concluded from comparing the $M=-1$ (small $\xi\approx a$) and $M=-0.2$ (large $\xi\approx 5a$) cases that show different decay lengths. The current flowing along ${\hat{\bf x}}$ direction at the $y=1$ edge $\langle J_{x}^{0}(y=1)\rangle$ also confirms this decaying behavior, as shown in Fig. 3. Thus we anticipate that this partially gated TI quantum dot can serve as a nonlocal edge current generator, with the range of the edge current controllable by the insulating gap $M$. Finally, we remark that in all cases, the continuity equation in Eq. (7) is satisfied $\dot{n}_{i}=0$ at any site $i$, provided that the current flowing along positive $J_{i,i+\delta}^{0}$ and negative bonds $J_{i,i-\delta}^{0}$ are both taken into account.

For 2D class AII, we use the prototype Bernevig-Hughes-Zhang (BHZ) model as an exampleBernevig et al. (2006), which has basis $\psi=\left(s\uparrow,p\uparrow,s\downarrow,p\downarrow\right)^{T}$ and the Dirac matrices

$\displaystyle\Gamma^{\ell}=\left\{\sigma^{z}\otimes s^{x},I\otimes s^{y},I\otimes s^{z},\sigma^{x}\otimes s^{x},\sigma^{y}\otimes s^{x}\right\}.$

(14)

where $\sigma^{b}$ and $s^{b}$ are Pauli matrices in the spin and orbital spaces, respectively. The continuum model reads

	$\displaystyle H({\bf k})$	$\displaystyle=$	$\displaystyle A\sin k_{x}\Gamma^{1}+A\sin k_{y}\Gamma^{2}$		(15)
		$\displaystyle+$	$\displaystyle\left(M-4B+2B\cos k_{x}+2B\cos k_{y}\right)\Gamma^{3}\;,$

which regularizes to give the square lattice model

	$\displaystyle H=\sum_{i\sigma}t\left\{-i\sigma c_{is\sigma}^{{\dagger}}c_{i+ap\sigma}-i\sigma c_{ip\sigma}^{{\dagger}}c_{i+as\sigma}+h.c.\right\}$
	$\displaystyle+\sum_{i\sigma}t\left\{-c_{is\sigma}^{{\dagger}}c_{i+bp\sigma}+c_{ip\sigma}^{{\dagger}}c_{i+bs\sigma}+h.c.\right\}$
	$\displaystyle+\sum_{i\sigma}\left(M+4t^{\prime}-\mu\right)c_{is\sigma}^{{\dagger}}c_{is\sigma}+\sum_{i}\left(-M-4t^{\prime}-\mu\right)c_{ip\sigma}^{{\dagger}}c_{ip\sigma}$
	$\displaystyle-\sum_{i\sigma\delta}t^{\prime}\left\{c_{is\sigma}^{{\dagger}}c_{i+\delta s\sigma}-c_{ip\sigma}^{{\dagger}}c_{i+\delta p\sigma}+h.c.\right\},$		(16)

where $t=A/2$ and $t^{\prime}=-B$. Here $I=\left\{s,p\right\}$ is the orbital index, $\delta=\left\{a,b\right\}$ denotes the lattice constant along the two planar directions, $\sigma=\left\{\uparrow,\downarrow\right\}=\left\{+,-\right\}$ is the spin index, $i=\left\{x,y\right\}$ denotes the planar position. We use the parameters $-t=t^{\prime}=-M=1$ and temperature $k_{B}T=0.03$. Once again the chemical potential $\mu$ will be crucial for the existence of edge spin current.

The band structure at $N_{y}=8$ is shown in Fig. 4. The edge states consist of counter propagating spins polarized along ${\hat{\bf z}}$, as can be seen from the expectation value of $\sigma^{z}$ of $|k_{x},n_{y}\rangle$ near the $y=1$ edge

$\displaystyle\tilde{m}_{k_{x},n_{y}}^{z}=\sum_{1\leq y\leq N_{y}/2}\sigma_{k_{x},n_{y}}^{z}.$

(17)

The color in the band structure of Fig. 4 (b) indicates the spin up (blue) and down (green) of $\tilde{m}_{k_{x},n_{y}}^{z}$. The branch of Dirac cone of positive group velocity has spin up and that of negative group velocity has spin down, as expected (note that the Dirac cone of black color is that localized at the other edge $y=N_{y}$). Besides, the topmost valence band at $k_{x}>0$ with negative group velocity has more spin up component, whereas that at $k_{x}<0$ with positive group velocity have more spin down component. This indicates at least some valence bands contribute to a spin current against that contributed from the edge states. We further quantify the two contributions by considering the local spin current operators constructed from the equation of motion of the spin density $m_{i}^{a}=\sum_{I}c_{iI\alpha}^{{\dagger}}\sigma_{\alpha\beta}^{a}c_{iI\beta}$

	$\displaystyle\dot{m}_{i}^{a}=\frac{i}{\hbar}\left[H,m_{i}^{a}\right]=-{\bm{\nabla}}\cdot{\bf J}_{i}^{a}=-\frac{1}{a}\sum_{\delta}\left(J_{i,i+\delta}^{a}+J_{i,i-\delta}^{a}\right).$
			(18)

Because the edge states are spin polarized along ${\hat{\bf z}}$, we investigate the local spin current operators $J_{x}^{z}\equiv J_{i,i+a}^{z}$ and $J_{y}^{z}\equiv J_{i,i+b}^{z}$ that are given by

	$\displaystyle J_{x}^{z}=\frac{a}{\hbar}\sum_{\sigma}t\left\{c_{is\sigma}^{{\dagger}}c_{i+ap\sigma}+c_{ip\sigma}^{{\dagger}}c_{i+as\sigma}+h.c.\right\}$
	$\displaystyle+\frac{a}{\hbar}\sum_{\sigma}t^{\prime}\left\{-i\sigma c_{is\sigma}^{{\dagger}}c_{i+as\sigma}+i\sigma c_{ip\sigma}^{{\dagger}}c_{i+ap\sigma}+h.c.\right\},$
	$\displaystyle J_{y}^{z}=\frac{b}{\hbar}\sum_{\sigma}t\left\{-i\sigma c_{is\sigma}^{{\dagger}}c_{i+bp\sigma}+i\sigma c_{ip\sigma}^{{\dagger}}c_{i+bs\sigma}+h.c.\right\}$
	$\displaystyle+\frac{b}{\hbar}\sum_{\sigma}t^{\prime}\left\{-i\sigma c_{is\sigma}^{{\dagger}}c_{i+bs\sigma}+i\sigma c_{ip\sigma}^{{\dagger}}c_{i+bp\sigma}+h.c.\right\}.$		(19)

The spatial profile of the edge spin current evaluated with an energy cutoff $E_{cut}$, as in Eq. (9), is essentially the same as that in Fig. 1 (b) for the Chern insulators, up to an overall prefactor. This is not surprising, since the BHZ model is basically two copies (spin up and down) of the Chern insulator, and we also see that the contribution from the valence bands cancels out exactly that from the edge states, rendering no edge spin current if the Dirac cone resides at the chemical potential $\mu=0$. The quantized edge spin current $\langle J_{x}^{z}\rangle_{h}=\sum_{1\leq y\leq N_{y}/2}\langle J_{x}^{z}\rangle$ as a function of chemical potential $\mu$ is shown in Fig. 4 (b), which has the same number of plateaux and edge spin current quantum as that in Eqs. (10) and (11), except one replaces the electron charge by the Bohr magneton $e\rightarrow\mu_{B}$ in these equations.

The equilibrium edge spin current $J_{x}^{z}$ does not show up in nonequilibrium transport or spin transport measurements, but the electric dipolar fieldHirsch (1990, 1999); Meier and Loss (2003); Schütz et al. (2003); Sun et al. (2004); Chen (2014) it produces should in principle be measurable. This can be understood by considering that a magnetic moment ${\bf m}$ moving with velocity ${\bf v}$ produces an electric dipole moment in the rest frame

$\displaystyle{\bf p}=\frac{\gamma}{c^{2}}{\bf v}\times{\bf m},$

(20)

where $\gamma=(1-v^{2}/c^{2})^{-1/2}$, and subsequently a dipolar field at position ${\bf r}$ from the moving magnetic moment

$\displaystyle{\bf E_{m}(r)}=\frac{3{\hat{\bf r}}({\bf p}\cdot{\hat{\bf r}})-{\bf p}}{4\pi\epsilon_{0}r^{3}}.$

(21)

Approximating the edge along ${\hat{\bf x}}$ as a long straight wire, the edge spin current quantum $\Delta J_{x}^{z}=\mu_{B}v_{F}/N_{x}a$ corresponds to an electric dipolar fieldMeier and Loss (2003)

$\displaystyle\Delta{\bf E}_{m}({\bf r})=\frac{\mu_{0}\Delta J_{x}^{z}}{2\pi r^{2}}(0,\cos 2\phi,-\sin 2\phi),$

(22)

where $\sin\phi=y/r$, $\cos\phi=z/r$, and $r=\sqrt{y^{2}+z^{2}}$. Assuming a quantum dot of size $N_{x}\sim 10^{2}$ and a typical Fermi velocity $v_{F}\sim 10^{5}$m/s, the electric dipolar field at a distance $r\sim$ nm away from the edge is $|\Delta{\bf E}_{m}|\sim 1$V/m. Thus two points that are $\left\{{\bf r}_{1},{\bf r}_{2}\right\}\sim$ nm away from the edge and are ${\bf r}_{1}-{\bf r}_{2}\sim$ nm apart would experience a voltage drop $\Delta V\sim$ nV, which should be measurable by local high impedance voltmeters, assuming the gate voltage does not affect the measurement.

The local spin current promoted by impurities and edge confining potential follows that discussed in Sec. II.3. The nonlocal edge spin current in partially gated BHZ quantum dot is also given by that of the Chern insulator in Fig. 3, and the extension of the edge spin current into the ungated area is again characterized by the decay length $\xi$. Note that the BHZ model at parameters $M\sim-0.01$eV and $\hbar v_{F}/a\sim$ eV is relevant to the HgTe quantum wellKönig et al. (2007, 2008) or recently proposed III-V semiconductor quantum wellCandido et al. (2018) at few nanometers thickness, which yields a decay length of the order of $\xi\sim 100a\sim 100$nm. Thus if a partially gated HgTe quantum well of area $\sim 100$nm$\times 100$nm can be fabricated, it can serve as a dissipationless nonlocal spin current generator of $\sim 100$nm range, which is fairly significant. In contrast, the recently proposed monolayer WTe₂ has $M\sim-0.1$eV and $\hbar v_{F}/a\sim 0.1$eVQian et al. (2014); Zheng et al. (2016); Tang et al. (2017), so $\xi\sim a$ is of the order of lattice constant, which may not be the ideal material for the nonlocal spin current generator.

We next consider 3D TIs in class AII, which are relevant to materials such as Bi₂Se₃ and Bi₂Te₃. These materials as made have a single edge state Dirac cone with the Dirac point located $\sim 0.1$eV away from the chemical potentialChen et al. (2009); Analytis et al. (2010); Zhang et al. (2010); Pan et al. (2011), but the position of Dirac point can be efficiently engineered by dopingHsieh et al. (2009); Zhang et al. (2011); Kondou et al. (2016). We will consider the theoretical model for the low energy sector described by the $\Gamma$-matricesZhang et al. (2009); Liu et al. (2010)

$\displaystyle\Gamma^{\ell}=\left\{\sigma^{x}\otimes\tau^{x},\sigma^{y}\otimes\tau^{x},\sigma^{z}\otimes\tau^{x},I_{\sigma}\otimes\tau^{y},I_{\sigma}\otimes\tau^{z}\right\},$

where the spinor is written as $\psi_{\bf k}=\left(c_{{\bf k}P1_{-}^{+}\uparrow},c_{{\bf k}P2_{+}^{-}\uparrow},c_{{\bf k}P1_{-}^{+}\downarrow},c_{{\bf k}P2_{+}^{-}\downarrow}\right)^{T}\equiv\left(c_{{\bf k}s\uparrow},c_{{\bf k}p\uparrow},c_{{\bf k}s\downarrow},c_{{\bf k}p\downarrow}\right)^{T}$, where $s$ and $p$ are abbreviations for the $P1_{-}^{+}$ and $P2_{+}^{-}$ orbitals in real materials. The low energy Hamiltonian given by ${\bf k\cdot p}$ theory is

	$\displaystyle\hat{H}$	$\displaystyle=$	$\displaystyle\left(M+M_{1}k_{z}^{2}+M_{2}k_{x}^{2}+M_{2}k_{y}^{2}\right)\Gamma^{5}$		(24)
		$\displaystyle+$	$\displaystyle B_{0}\Gamma^{4}k_{z}+A_{0}\left(\Gamma^{1}k_{y}-\Gamma^{2}k_{x}\right),$

where we keep only lowest order terms. The regularization on a cubic lattice yields

	$\displaystyle H=-\sum_{iI\sigma}\mu c_{iI\sigma}^{{\dagger}}c_{iI\sigma}+\sum_{i\in TI,\sigma}\tilde{M}\left\{c_{is\sigma}^{{\dagger}}c_{is\sigma}-c_{ip\sigma}^{{\dagger}}c_{ip\sigma}\right\}$
	$\displaystyle+\sum_{i\in TI,I}t_{\parallel}\left\{c_{iI\uparrow}^{{\dagger}}c_{i+a\overline{I}\downarrow}-c_{i+aI\uparrow}^{{\dagger}}c_{i\overline{I}\downarrow}+h.c.\right\}$
	$\displaystyle+\sum_{i\in TI,I}t_{\parallel}\left\{-ic_{iI\uparrow}^{{\dagger}}c_{i+b\overline{I}\downarrow}+ic_{i+bI\uparrow}^{{\dagger}}c_{i\overline{I}\downarrow}+h.c.\right\}$
	$\displaystyle+\sum_{i\in TI,\sigma}t_{\perp}\left\{-c_{is\sigma}^{{\dagger}}c_{i+cp\sigma}+c_{i+cs\sigma}^{{\dagger}}c_{ip\sigma}+h.c.\right\}$
	$\displaystyle-\sum_{i\in TI,\sigma}M_{1}\left\{c_{is\sigma}^{{\dagger}}c_{i+cs\sigma}-c_{ip\sigma}^{{\dagger}}c_{i+cp\sigma}+h.c.\right\}$
	$\displaystyle-\sum_{i\in TI,\delta,\sigma}M_{2}\left\{c_{is\sigma}^{{\dagger}}c_{i+\delta s\sigma}-c_{ip\sigma}^{{\dagger}}c_{i+\delta p\sigma}+h.c.\right\},$		(25)

where $\tilde{M}=M+2M_{1}+4M_{2}$, $t_{\parallel}=A_{0}/2$, $t_{\perp}=B_{0}/2$, $I=\left\{s,p\right\}$ and $\overline{I}=\left\{p,s\right\}$ are the orbital indices, $\delta=\left\{a,b,c\right\}$ denotes the lattice constants, $\sigma=\left\{\uparrow,\downarrow\right\}$ is the spin index, and $\mu$ is the chemical potential crucial for the existence of the surface spin current. We use the band parameters $t_{\parallel}=t_{\perp}=M_{1}=M_{2}=-M=1$ such that conclusions can be drawn by numerics done in a accessible lattice size of the order of $N_{x}\times N_{y}\times N_{z}\sim 10\times 10\times 10$, but we emphasize that the conclusions are robust against changing the parameters within the same order of magnitude. Choosing OBC along ${\hat{\bf z}}$ and PBC along ${\hat{\bf x}}$ and ${\hat{\bf y}}$, the surface state of momentum ${\bf k}_{\parallel}=(k_{x},k_{y})$ is spin polarized along ${\hat{\bf z}}\times{\hat{\bf k}}_{\parallel}$, which speculates a spin current $J_{x}^{y}$ flowing along ${\hat{\bf x}}$ and polarized along ${\hat{\bf y}}$, and a $J_{y}^{x}$ flowing along ${\hat{\bf y}}$ and polarized along ${\hat{\bf x}}$ of equal magnitude. In Fig. 5 (a), we use the band structure along $k_{x}$ at $k_{y}=0$ to elaborate the spin-momentum locking, where the sign of the spin polarization along ${\hat{\bf y}}={\hat{\bf z}}\times{\hat{\bf x}}$ close to the surface at $z=1$

$\displaystyle\tilde{m}_{k_{x},k_{y}=0,n_{z}}^{y}=\sum_{1\leq z\leq N_{z}/2}\sigma_{k_{x},k_{y}=0,n_{z}}^{y}.$

(26)

is indicated by blue (positive) and green (negative). Once again the spin polarization of the Dirac cone is evident, and it is also clear that some valence bands have the same spin polarization as the Dirac cone but opposite group velocities, hence producing a spin current against that produced by the surface states.

To quantify the effect of valence bands, we turn to the spin current operators constructed from Eq. (18). The spin currents flowing along the positive bonds $J_{x}^{y}\equiv J_{i,i+a}^{y}$ and $J_{y}^{x}\equiv J_{i,i+b}^{x}$ are given by

	$\displaystyle J_{x}^{y}=\frac{ia}{\hbar}\left\{it_{\parallel}\sum_{I}\left[c_{iI\uparrow}^{{\dagger}}c_{i+a\overline{I}\uparrow}+c_{iI\downarrow}^{{\dagger}}c_{i+a\overline{I}\downarrow}\right]-h.c.\right\}$
	$\displaystyle+\frac{ia}{\hbar}\left\{iM_{2}\sum_{\sigma}\left[\sigma c_{is\sigma}^{{\dagger}}c_{i+as\overline{\sigma}}-\sigma c_{ip\sigma}^{{\dagger}}c_{i+ap\overline{\sigma}}\right]-h.c.\right\}$
	$\displaystyle J_{y}^{x}=\frac{ib}{\hbar}\left\{it_{\parallel}\sum_{I}\left[-c_{iI\uparrow}^{{\dagger}}c_{i+b\overline{I}\uparrow}-c_{iI\downarrow}^{{\dagger}}c_{i+b\overline{I}\downarrow}\right]-h.c.\right\}$
	$\displaystyle+\frac{ib}{\hbar}\left\{M_{2}\sum_{\sigma}\left[-c_{is\sigma}^{{\dagger}}c_{i+bs\overline{\sigma}}+c_{ip\sigma}^{{\dagger}}c_{i+bp\overline{\sigma}}\right]-h.c.\right\}.\;\;\;\;$		(27)

where $I=\left\{s,p\right\}$ and $\overline{I}=\left\{p,s\right\}$. Using Eq. (9), we obtain a finite spin current as $E_{cut}<|M|$ that includes only the surface states, and eventually vanishes at large $E_{cut}$ that includes all the bands. The spin current is again finite when the chemical potential $\mu$ is nonzero, and changes sign when $\mu$ sweeps through the Dirac point, as shown in Fig. 5 (b). Thus we also expect that several factors in realistic 3D TIs, such as impuritiesWray et al. (2011), surface band bendingBahramy et al. (2012), and local electron and hole puddlesBeidenkopf et al. (2011), helps to promote the surface spin current, since they shift Dirac point away from the chemical potential globally or locally. As an example, the 3D vortex of the spin currents formed around a single impurity is shown in Fig. 5 (c), where one sees that the spin current $J^{a}$ that is polarized along ${\hat{\bf a}}$ wraps around the ${\hat{\bf a}}$ axis. We should also emphasize that despite the absence of the equilibrium edge spin current at $\mu=0$, the edge states can still perform nonequilibrium transport and spintronic effects, such as the conductance due to a bias voltageZhou et al. (2017) and the current-induced spin accumulationChen (2020); Zegarra et al. at $\mu=0$. This is because these nonequlibrium spintronic effects only involve the edge states near the chemical potential due to the derivative of the Fermi function, and hence the valence bands are irrelevant.

The quantization of edge spin current $\langle J_{x}^{y}\rangle_{h}=\sum_{1\leq z\leq N_{z}/2}\langle J_{x}^{y}\rangle=-\langle J_{y}^{x}\rangle_{h}$ as a function of $\mu$ is again evident, as shown in Fig. 5 (b). The number of plateaux $N_{p}$ and the surface spin current quantum $\Delta J_{x}^{y}$ are irregular functions of $\mu$ (the steps are not evenly spaced), and changes with the shape of the quantum dot, as can be seen by comparing the case of $N_{x}\times N_{y}\times N_{z}=12\times 12\times 8$ with that of $12\times 8\times 8$. This is presumably due to the irregularity in the number of states that can fit into the Dirac cone, which is best described in a cylindrical coordinate, but the lattice is cubic. Nevertheless, the trend of more plateaux and smaller surface spin current quantum at larger quantum dot, similar to the 2D cases described by Eqs. (10) and (11), is qualitatively correct.

The surface spin current also produces an electric dipolar field according to Eqs. (20) and (21). Consider a cubic quantum dot of 3D TI whose top surface area $L^{2}$ is confined in the region $-L/2\leq x\leq L/2$ and $-L/2\leq y\leq L/2$. A uniform $\langle J_{x}^{y}\rangle$ (e.g., produced by doping) flowing in the top surface can be regarded as a current $I_{m}=\mu_{B}\langle J_{x}^{y}\rangle/a=m\rho v$ of magnetic moment ${\bf m}$ of planar density $\rho$ polarized in ${\hat{\bf y}}$ and moving along ${\hat{\bf x}}$. A point located at $z$ above the center of the quantum dot, as shown schematically in Fig. 5 (d), experiences an electric field along ${\hat{\bf z}}$ direction

	$\displaystyle E^{z}(z)$	$\displaystyle=$	$\displaystyle\frac{\mu_{0}I_{m}}{4\pi}\left\{\frac{8\sqrt{2}L^{2}}{\sqrt{L^{2}+2z^{2}}(L^{2}+4z^{2})}\right.$		(28)
			$\displaystyle+\frac{4}{z}{\rm arccot}\left(\frac{2\sqrt{2}z\sqrt{L^{2}+2z^{2}}}{L^{2}}\right)$
			$\displaystyle+\left.-\frac{4}{z}{\rm arctan}\left(\frac{L^{2}}{2\sqrt{2}z\sqrt{L^{2}+2z^{2}}}\right)\right\}.$

In the limit $z\ll L$, i.e., very close to the top surface of a large quantum dot, one has $E^{z}\approx\mu_{0}I_{m}/L$. Given that $\langle J_{x}^{y}\rangle$ in Fig. 5 (b) is of the order of $0.01t/\hbar$, the corresponding magnetization current is $I_{m}\sim 0.1$Cm/s². For a quantum dot of dimension $L\sim 100a\sim 10^{-7}$m, the magnetization current produces an electric field $E^{z}\approx\mu_{0}I_{m}/L\sim 1$V/m near the top surface. Thus for two points $\left\{z_{1},z_{2}\right\}$ that are nm away from the surface and are nm apart, as that in in Fig. 5 (d), the voltage difference between them is $\Delta V\sim$nV, which is readily measurable. Interestingly, Eq. (28) also implies that the surface spin current of an infinitely large TI $L\rightarrow\infty$ does not produce a electric dipolar field $E^{z}\rightarrow 0$. So the aforementioned voltage difference can only be observed in TIs of finite size.

The band structures resulted from regularizing the low energy Dirac Hamiltonian to the entire BZ, as shown in Fig. 1 (a), 4 (a), and 5 (a) certainly cannot capture the band structures in real TIs. Thus our conclusion of the exact cancellation of the edge current by valence bands is valid only within these regularized lattice models. In real TIs, we suspect that the valence bands can contribute to an edge current flowing in either direction depending on details of the material. To settle this issue, a feasible approach would be that one starts from the realistic tight-binding models constructed from certain ab initio calculations for specific materialsZhang et al. (2009); Liu et al. (2010), extract the current operators from Eqs. (7) and (18), then calculate the equilibrium edge currents from Eq. (9). How the edge current depends on the chemical potential $\mu$ can be easily verified, and the effect of inhomogeneity can be simulated by adding appropriate impurities or edge confining potentials into these realistic lattice models. Although it is not our purpose to perform this kind of calculation in the present work, the important point drawn from our simplified approach is that one cannot judge the amount and direction of the edge currents solely from the edge state Dirac cone, and the valence bands generally contribute to all equilibrium properties of TIs, even for those at the edge or surface.