Spin generation via bulk spin current in three-dimensional topological insulators

To date, spin generation in three-dimensional topological insulators is primarily modelled as a single-surface phenomenon, attributed to the momentum-spin locking on each individual surface. In this article, we propose a mechanism of spin generation where the role of the insulating yet topologically non-trivial bulk becomes explicit: an external electric field creates a transverse pure spin current through the bulk of a three-dimensional topological insulator, which transports spins between the top and bottom surfaces. Under sufficiently high surface disorder, the spin relaxation time can be extended via the Dyakonov–Perel mechanism. Consequently, both the spin generation efficiency and surface conductivity are largely enhanced. Numerical simulation confirms that this spin generation mechanism originates from the unique topological connection of the top and bottom surfaces and is absent in other two-dimensional systems such as graphene, even though they possess a similar Dirac cone-type dispersion.

T opological insulators (TIs) have attracted world-wide attention because of their intriguing fundamental physics and exciting application opportunities in spintronics 1 . Three-dimensional (3D) TIs 2,3 are of particular technological importance as the unique spin generation can be realized in single crystals rather than in complex heterogeneous structures 4 . TIs are considered as efficient spin generators 5 , yet the spin generation is generally regarded as a pure surface phenomenon. Namely, the electronic momentum and spin are locked at the TI surface, and a net charge current leads to a net spin polarization at the surface, whose magnitude is directly proportional to the charge current 6 . In this view, all physics occur independently at the top and bottom surfaces of a TI and the role of the bulk is passive, which simply separates the top and bottom surfaces. The surface conductivity is understood through density of states and scattering rate, just like in other two-dimensional (2D) systems such as graphene and 2D electron gas. The conductivity behaviour governs the spin generation on the surface of a 3D TI, and spin accumulation is merely a side product of conductivity. Although this interpretation of spin generation in TIs is most mathematically straightforward, it is far from satisfactory in the sense that the most amazing feature of a TI-surface-bulk correspondence does not explicitly enter this physical picture.
On the other hand, there is an alternative viewpoint of spin generation. The external electric field induces a transverse pure spin current through the bulk, which acts as a bridge for transporting spins between top and bottom surfaces. Opposite spins are thus accumulated on the two surfaces, which lead to charge current in the same direction of the electric field because of the opposite chirality of the momentum-spin textures on the top and bottom surfaces (Fig. 1a). An empirical formula for the bulk spin current can be written down as where j s is the spin current density, E is the electric field and s s is the spin Hall conductivity tensor 7 . A system that is electrically insulating but can carry a pure spin current is termed a spin Hall insulator 8 . The bulk of a 3D TI has been demonstrated to be a spin Hall insulator because of its Z 2 topological order 9 . Analogous to Hall effect, the transverse spin Hall current leads to surface spin accumulation in a slab geometry. Yet unlike electric charge, spin is usually a nonconserved quantity in a spin Hall insulator. The ultimate spin accumulation induced on the surface closely depends on the spin relaxation mechanism. In the low disorder limit mt=' ) 1 with m being the Fermi level, t being the momentum relaxation time, it has been demonstrated that the spin relaxation time t s on the surface of a 3D TI is identical to the momentum relaxation time t because of the momentum-spin locking, and the traditional Dyakonov-Perel spin relaxation is absent 10 . Charge-spin dynamics in the high disorder limit mt=' t1, however, has rarely been discussed in the literature so far.
The exact behaviour of these spin transport coefficients under high disorder is of crucial importance to the application of 3D TI-based spintronic devices, because unlike the bulk, surface is extremely vulnerable to various kinds of defects, especially when placed in ambient environment. Even for a material that is generally considered 'inert', the top most layer of atoms could still suffer from high concentration of impurities 11,12 . A clear physical model of charge and spin transport in this case is highly desired for the design of novel 3D TI-based spintronic devices.
In this article, we demonstrate that the unique topological connection of the surface bands in a 3D TI demands a pure spin current through the insulating bulk. Spins are thus generated and accumulated on the two surfaces. Sufficiently high surface disorder can suppress spin relaxation and result in an increase of the spin relaxation time t s in a manner similar to the traditional Dyakonov-Perel mechanism. Consequently, both electro-spin susceptibility k yx (surface spin density s y divided by the electric field E x ) and the electrical conductivity s xx should increase with the increase of disorder, as illustrated in Fig. 1b.

Results
Spectral function and electrical conductivity. To begin with, we consider a realistic four-band tight-binding model 13 built on a slab of a tetragonal lattice, as shown in Fig. 2a. The slab is infinite in xy directions and has a total number of N ¼ 10 layers in the z direction (c-axis). With four states on each site, the bulk Hamiltonian in the 3D k-space is where a i (i ¼ x, y, z), b are the Dirac matrices, a and a z are the lattice constants in the xy and z directions, D is the mass term and A, A z , B, B z are nearest-neighbour hopping amplitudes. In the slab configuration, inverse Fourier transform is performed in z-direction to comply with the finite thickness. We use a typical 3D TI Bi 2 Se 3 as our prototype and adopt parameters as obtained in ref. 2 to best fit the band structure of Bi 2 Se 3 . The resultant band structure of the surface is shown in Fig. 2b, which clearly has a Dirac cone near the G point. Owing to the z-inversion symmetry of the slab, all bands are doubly degenerate.
To account for surface disorder, atoms in the top and bottom layers of the slab are subject to a typical kind of impurityvacancies. Each site at the surface has a probability of c to be occupied by a vacancy where the on-site energy is brought to infinity so as to forbid electrons from this site. As c may not be small, the first Born approximation does not apply. Here we adopted the coherent potential approximation (CPA) method for binary alloys 14,15 in computing the Green's function G(k, o) and self-energy S(o). A typical spectral function À (1/p)ImG(k, o) obtained by CPA at impurity concentration c ¼ 0.001 is plotted in Fig. 2c. The evolution of the spectral function with increasing impurity concentration is consistent with results obtained in ref. 13. Subsequently, transport coefficients were calculated via the standard linear response theory. More details of the numerical simulation can be found in the Methods section. Figure 2c,d shows the electrical conductivity calculated from the CPA Green's function via the Kubo-Greenwood formalism. With an impurity concentration c ranging from 5 Â 10 À 4 to 0.5, the Fermi level dependence of conductivity gets weaker and the magnitude of conductivity reaches a minimum at around c ¼ 0.006. Further increasing the impurity concentration leads to an increase of conductivity at a given Fermi level position. Such anomalous increase of conductivity with impurity concentration is difficult to understand based on a single surface model 16,17 , which suggests the essential role of the bulk of a 3D TI in surface conduction. In the following, we reveal that the anomalous increase of conductivity is a signature of a different type of spin dynamics and manifests a spin generation mechanism in 3D TIs.
Spin generation via bulk spin current. To start discussions on spin dynamics, we notice that spin is not a predefined quantity in Hamiltonian (equation 2). Although common 3D TIs such as Bi 2 Se 3 are known to have chiral spin texture on the surface states, the spin polarization is not 100% (ref. 18). Nevertheless, one can always talk about a pseudo-spin, which is defined to exactly match the energy eigenstates and has all essential features of the real spin 19 . Here we take the definition It can be verified that such definition satisfies all symmetry requirements of the real spin and shows a chiral spin texture near the G point, as shown in Fig. 3a. Note that a unique spin polarization can be specified for all points in the k-space except for those time-reversal invariant momenta (TRIM), where the Kramers theorem asserts the degeneracy of the two opposite spin polarizations.
With the above definition of spins, the bulk spin current density j s zy denoting the transport of y-spins in the z-direction can be defined as Each atom on these layers has a probability of c to be occupied by a vacancy and (1 À c) to be intact. The on-site energy of a vacancy is brought to infinity to forbid electrons from entering this site. (b) The energy dispersion of the surface branch of a clean system. A Dirac cone exists around the G point. Parameters used for simulation: ARTICLE where P þ y and P À y are spin projection operators in the þ y and À y directions, v z is the velocity operator in the z direction and O 3 is the volume of the slab serving as a normalization factor. With the above expression, the spin Hall conductivity s s zyx defined by can be calculated via the standard linear response theory. Owing to the even symmetry of spin current j s under time reversal J , the resultant expression for the spin Hall conductivity is different from that for the electrical conductivity (Kubo-Greenwood formula), but contains a term which involves all states below the Fermi level. This term has been thoroughly reviewed in ref. 20 for the calculation of electrical conductivity in a J -symmetry broken system and also discussed in a recently published article 21 for the calculation of spin Hall conductivity. The emergence of this term in our system indicates the non-dissipative nature of the spin current, which has already been demonstrated possible for a wide class of traditional semiconductors 22,23 . We leave the details of derivation to the Methods section and plot the calculated spin Hall conductivity s s zyx in Fig. 3b. It is clear that the magnitude of the spin Hall conductivity is independent of both the Fermi level position (must be within the bulk bandgap) and the surface impurity concentration.
Although the existence of a bulk spin current in 3D TI has been predicted analytically through topological argument 9 , a visualization of the spin transfer mechanism is not yet available so far. Neither has its relevance to the transport behaviour of the gapless surface states been studied ever. In the following, we present an intuitive picture of the spin transfer in a 3D TI slab and uncover its close relationship with the surfacebulk correspondence of a 3D TI.
We notice that despite the chiral spin texture over the entire Brillouin zone (BZ), only states with small magnitude of momentum are truly localized on the surface. Figure 3c,d shows the evolution of electronic wave functions as the wave vector k approaches the BZ boundary from the G point. It is seen that, beyond a certain point, electronic wave functions become extended through the entire bulk and the surface band has essentially merged into bulk bands. States beyond this merging point should be classified as bulk states although they lie on the same branch of energy sub-band as true surface states.
Imagine applying a weak electric field to this system in the þ x direction and examine the G À X line in the extended BZ view. Because of the inversion symmetry in the z direction, all bands are doubly degenerate. We notice, however, in order for the spin texture to be continuous, every top surface branch must be connected to the adjacent bottom surface branch and vice versa, as shown in Fig. 3e. This alternating structure exists across all TRIM points in our system, and is distinctively different from a normal band, which smoothly connects to itself at the BZ boundary. Consider an electron on the bottom surface with its spin polarized in þ y direction. Under the driving of the electric field, this electronic state drifts to À x direction in k space and merges into the bulk valance band. Upon further drifting, this electron finally enters the top surface with its spin in þ y direction unchanged. Simultaneously, an electron with spin polarized in the À y direction will drift from the top surface to bottom surface. The drift motion across the X point is similar to the Klein tunnelling of Dirac Fermions in the sense that, in order for a certain spin to be continuous, the electron must tunnel to another band rather than return to its original band. Overall, each of these processes corresponds to a unit spin-pair exchange between the bottom surface and the top surface. Thus, a longitudinal electric field induces a transverse pure spin current through the bulk, which plays the role of a spin injector for the two surfaces, as described in Fig. 1a. During this process, it is essential that the Fermi level lies within the gap, because there exists another pair of merging points near the conduction band edge. If the Fermi level is above these points as well, there would be an opposite process that leads to the cancelation of net spin current. This is of course consistent because the system in this case is not a TI any more.
Spin relaxation on the surface. Unlike charge, spin is not a conserved quantity in our system. The spins injected onto the surface suffer from immediate relaxation. The scenario is slightly different from both the Hall effect and the 2D spin Hall effect. In 3D, the spin relaxation is actually necessary for the system to reach a steady state, as detailed in Supplementary Note 1. The ultimate spin density accumulated on the surface is determined by the spin relaxation time t s . In the following, we provide an intuitive physical picture of the spin relaxation dynamics under the eigenbasis defined by H 0 . More details of this picture can be found in Supplementary Note 2. This is essentially an interaction picture that splits the Hamiltonian into a free part H 0 and an interaction part U. Because of the momentum-spin locking, each electron senses an effective magnetic field B eff according to its wave vector k. When an electron with its spin aligned with B eff suffers from a momentum change :Dk because of the scattering of an impurity potential, its spin may no longer align with the new B eff . If scatterings are rare, that is, the time it takes for momentum to change by a unit amount is long, the adiabatic perturbation theory predicts that the new spin must evolve to the new energy eigenstate, that is, rotate to the direction of the new B eff . If scatterings are frequent, however, the spin does not have time to follow B eff and will precess about the instantaneous B eff , as shown in Fig. 4a. Frequent scatterings constantly change the precession axis and the spin ends up doing a random walk on a unit sphere, as shown in Fig. 4b,c. The more frequent momentum scattering is, the less effective the random walk is, and the spin will preserve its original direction for a longer time. Therefore, the spin relaxation time t s inversely depends on the momentum relaxation time t, just like in the traditional Dyakonov-Perel mechanism 24 . One point to note is that as disorder is only present on the surface, states outside the dashed box of Fig. 3e are unaffected by scattering and the previously discussed spin transfer mechanism remains valid even under strong surface disorder.
If the above physical picture is correct, the ultimate spin density accumulated on the surface should increase with the increase of disorder. We calculated the spin relaxation time t s and electro-spin susceptibility k yx via standard linear response theory. The results are shown in Fig. 4d,e, which perfectly agree with the expectation. Combined with the fact that velocity operator is proportional to spin on the surface, it is not difficult to understand the anomalous increase of conductivity as well.

Discussion
It is worth noting that, in contrast to a common belief, the behaviour of conductivity of a Dirac system under high disorder is not simply governed by the dispersion relationship. The anomalous increase of conductivity is closely related to the spin generation and relaxation mechanism. To illustrate this point, we calculated the conductivity of a single layer of atoms within the same model (Fig. 5a). By setting the mass term D ¼ 0, the band structure of this 2D system has almost identical shape as the previous surface state, as shown in Fig. 5b, but the spin generation mechanism discussed above is obviously absent. With the increase of the impurity concentration, the conductivity of such system monotonically decreases towards zero, even in the high disorder range. Similar behaviour has also been shown in graphene 25,26 , where no anomalous increase of conductivity was found.
The spin dynamics in a 3D TI under strong surface disorder makes delicate connection with the scenario of a 2D TI. In fact, it is straightforward to apply the previously discussed spin transfer mechanism to a 2D TI, and obtain the universal quantized spin Hall conductivity s s 2D ¼e= p' ð Þ. Different from the 3D case, in the edge channel of a 2D TI, the spin (understood as pseudo-spin when necessary) is a conserved quantity, which does not relax. Consequently, the electro-spin susceptibility is infinite, which means no external field is needed to support the edge spin accumulation and charge current. The quantized and finite channel conductivity e 2 /h is actually a contact effect, while the channel itself is dissipationless 27 . On the surface of a 3D TI, however, spin is not conserved due to an additional angular degree of freedom of the wave vector k. The electro-spin susceptibility is thus finite and transport is dissipative. Strong surface disorder greatly suppresses spin relaxation and brings the system closer to the situation of a 2D TI, leading to a more efficient spin generation. In this view, strong surface disorder can be beneficial for spintronic devices, in contrast with the common belief. Technically, it is obviously of more convenience to induce high disorder on a surface than to make it pure and pristine.
Experimentally, the anomalous increase of conductivity in 3D TIs has already been hinted by results from several groups, yet researchers do not generally regard it as an intrinsic property of the TI surface. Field effect measurements in TIs have often shown a high minimum conductivity even when the Fermi level is tuned to the charge neutral point [28][29][30][31][32] , which is much larger than expected by normal transport theory assuming low disorder 16,17 . Although they are often attributed to bulk conduction 29 or surface electron puddle formation 30,33 , a closer look at these models reveals several problems, which are discussed in details in Supplementary Note 3. Theory presented in this paper, however, provides a simple and natural way to understand these observations. Recently, it has been directly observed in exfoliated BiSbTeSe 2 nanoflakes that after argon ion milling treatment to create more surface defects, the sample becomes more conductive 34 , although this effect was not understood.
Our simulation results also suggest that under strong surface disorder, with the magnitude of transport coefficients increased, their Fermi level sensitivity has dropped, which is also beneficial for the design of 3D TI-based spintronic devices. This is because in ambient environment, the Fermi level on the surface is subject to unintentional change because of contamination and degradation 35,36 . Sensitive Fermi level dependence renders the device less stable and robust in air.

Methods
CPA and self-energy. The numerical calculation results were obtained via CPA and the standard linear response theory. In the following, we use convention : ¼ 1.
The Hamiltonian of a clean lattice in the main text can be written in a block diagonal form where each element in this matrix is a 4 Â 4 matrix. The N rows and columns Overall, the Hamiltonian is 4N Â 4N. Diagonalizing this matrix gives the band structure of a clean lattice without impurities.
We consider impurities of an on-site scalar potential U on the top and bottom layers, which takes the form of U ¼ diag u; u; u; u; 0; 0; . . . ; 0; 0; u; u; u; u ½ ð 8Þ Each site has a probability of c subject to potential U and probability (1 À c) subject to potential 0, which makes the entire system essentially a binary alloy. Such configuration assumes non-physical correlation between the appearance of an impurity on the top and bottom layers. Yet from practical consideration, for a sufficiently thick slab, the crosstalk between the top and bottom layers should vanish, which justifies the binary alloy model of the above.  (1 À c) to be intact. The on-site energy of a vacancy is brought to infinity. Only nearest-neighbour hopping is allowed. (b) The energy dispersion of this single layer system. Simulation parameters are identical to before except for the mass term D is set to zero in order to have a Dirac cone near the G point. The band structure looks almost identical to the surface bands of the 3D TI slab discussed before. (c) The DC conductivity s xx (0) of this single layer versus Fermi level at different impurity concentrations. Despite the similarity in band structures, this single layer does not have anomalous increase of conductivity at high disorder, in contrast with our previous system.
In CPA, the configurationally averaged impurity potential is denoted by a k-independent self-energy S(o). In the binary alloy case, S(o) is determined by the iterative equation where G(o) is the on-site Green's function Because of the symmetry of this problem, S(o) is actually a scalar on the top and bottom layers only. The real and imaginary parts of the self-energy are plotted in Supplementary Fig. 1.
Conductivity. The velocity operator v takes the form v ¼ @H @k ð11Þ According to the Kubo formula, the direct current (DC) conductivity at T ¼ 0 is The conductivity of the top surface is half of this value Definition of spin operators and the chiral spin texture. Spin is not a predefined quantity in the original Hamiltonian. Unless we are interested in a direct coupling to an external magnetic field, the four state basis vectors can be thought to have arbitrary spin polarizations, which is essentially a 'pseudo-spin'. Even in real TIs such as Bi 2 Se 3 , the pseudo-spin does not exactly match the real spin. Nevertheless, the definition of a 'pseudo-spin' S must satisfy a couple of restrictions: (i) S is an Hermitian operator, S w i ¼S i . (ii) The components of S satisfy the anti-commutation rules, S i S j ¼d ij þ iE ijk S k . (iii) S is a pseudo-vector. It transforms like a vector under in-plane (xy) rotation but does not flip sign under space inversion, bS i b ¼ S i . (4) To comply with the chiral surface spin texture, we require S be polarized along y-direction when k y ¼ 0. Thus, [S y , H(k x , 0)] ¼ 0.
Note that we do not include :/2 in our definition, so this pseudo-spin has dimension 1 instead of angular momentum. It can be easily verified that the following expressions are a good representation of 'spin' in our system.
The spin in an arbitrary direction y within the xy plane is thus If there exists an angle y for an arbitrary point in k-space (k x , k y ) such that then S(y) and H 0 (k x , k y ) have common eigenstates and the energy eigenstates can be assigned a unique spin polarization. It is not difficult to obtain When the magnitude of k is small, we recover to the well-known chiral spin texture tan y % À k x k y ð21Þ whereas when k gets large, the spin texture deforms to adapt to the tetragonal symmetry of the BZ, as shown in the main text.
Surface spin density and electro-spin susceptibility. The spin density accumulated on the top surface is At T ¼ 0 the electro-spin susceptibility is calculated similar to the electrical conductivity Note that the definition of s top y contains a factor of 1/O and thus the above expression is actually independent of the box size O.
Bulk spin current and spin Hall conductivity. The z-position operator is z ¼ a z Á diag À 1; À 1; À 1; À 1; À 2; À 2; À 2; À 2; À 3; À 3; À 3; À 3; . . . ½ Thus, the z-velocity operator As we are interested in the flux across a certain intermediate layer, we restrict the velocity operator to be only between two adjacent layers in the middle To discuss the spin current, we need to define a spin projection operator P þ y ¼ I S y ¼ þ 1 S y ¼ þ 1 ð29Þ P À y ¼ I S y ¼ À 1 S y ¼ À 1 ð32Þ The spin current operator is then j s zy ¼ P þ y v m z P þ y À P À y v m z P À y = a z O ð Þ ð35Þ -----------------------------Spin relaxation time. Imagine applying a pulse electric field to our system Empirically, the induced spin has the asymptotic form s t ð Þ ¼ s 0 e À t ts y t ð Þ ð43Þ where the function y(t) ¼ 0 when to0 and y(t)-1 as t-þ N. Consider the Fourier transformŝ For sufficiently small o, the details of the rising part of s(t) characterized by y(t) becomes unimportant, thus we replace y(t) with the step function Y(t) and obtainŝ The Fourier transform of the electric field is just a constant Thus, the electro-spin susceptibility k o ð Þ / 1 io À 1 ts ð49Þ and the spin relaxation time can be extracted as It is necessary to point out that the spin relaxation process can be thought as an eigen mode with a complex frequency o* on the lower half plane. o* is a pole of the response function w(o*) ¼ N. The spin relaxation time is determined by the imaginary part of the pole closest to the real axis The equation 51 is based on the low-frequency expansion of w(o), which may not give the exact pole position, and ref. 10 has shown that for a perfect Dirac-cone dispersion, the spin relaxation time found by the exact pole is twice as the value found by low frequency expansion. Nevertheless, apart from an order 1 factor, low-frequency expansion should give a reasonable estimate of the true spin relaxation time.