LETTER pubs.acs.org/NanoLett
Helical States of Topological Insulator Bi2Se3 Yonghong Zhao,*,†,‡ Yibin Hu,† Lei Liu,§ Yu Zhu,§ and Hong Guo*,† †
Centre for the Physics of Materials and Department of Physics, McGill University, Montreal, Quebec H3A 2T8, Canada College of Physics and Electronic Engineering, Sichuan Normal University, Chengdu 610068, China § Nanoacademic Technologies Inc. 7005 Boulevard Taschereau, Brossard, Quebec J4Z 1A7, Canada ‡
ABSTRACT: We report density functional theory analysis of the electronic and quantum transport properties of Bi2Se3 topological insulator, focusing on the helical surface states at the Fermi level EF. The calculated Dirac point and the tilt angle of the electron spin in the helical states are compared quantitatively with the experimental data. The calculated conductance near EF shows a V-shaped spectrum, consistent with STM measurements. The spins in the helical states at EF not only tilts out of the two-dimensional plane, they also oscillate with a 3-fold symmetry going around the two-dimensional Brillouin zone. The helical states penetrate into the material bulk, where the first quintuple layer contributes 70% of the helical wave functions. KEYWORDS: Topological insulator, helical states, SOI, NEGF-DFT
T
he notion of topological insulator (TI) has attracted tremendous recent attention in communities of physics,1,2 chemistry,3 materials science,4 and electrical engineering.5 TI has an energy gap in its bulk band structure and metallic helical states on its surface.6,7 In the helical states, the direction of electron spin is locked to perpendicular to the electron momentum k: electrons moving to positive k have their spins pointing to one direction while those moving to the negative k have spins pointing to the opposite direction, and the helical states are protected by the time-reversal symmetry.6,7 The peculiar electronic structure of TI is caused by a strong spinorbit interaction (SOI) which, indeed, respects the time-reversal invariance. It is believed that Bi2Se3 and Bi2Te3 are 3D TIs.810 Nanowires and nanoribbons of these materials can now be fabricated by chemical means.3,4 Measurements by spin and angle resolved photoemission spectroscopy (ARPES) on these materials revealed several properties9 predicted by theoretical TI models11 such as the existence of a Dirac cone and the locking of spinmomentum. Electrical transport measurements on patterned nanoribbons of Bi2Se3 have also been carried out4 although, to the best of our knowledge, quantum transport that is only mediated by the helical surface states has yet to be achieved for 3D TIs. Scanning tunneling microscopy (STM)12,13 did show that the tunnel conductance to the Bi2Se3 surface has a V-shaped spectrum which is consistent with the Dirac band dispersion.13 On the theoretical side, density functional theory (DFT) has been applied to calculate electronic band structure of Bi-based TIs in bulk6,14 or in the slab geometry.1517 In particular, helical states in the slab geometry, the behavior of spin operators, Sharvin conductance, and the surface density of states (DOS)15,17 were calculated by planewave DFT methods. Interesting properties of 3D TI occur in the energy range of the material band gap that includes the Dirac point. ARPES data r 2011 American Chemical Society
for Bi2Se3 and Bi2Te3 showed that the Dirac point of these materials lies below the Fermi level810 EF. Since equilibrium quantum transport properties are contributed by electrons at EF, helical states at EF warrant very careful study in addition to those at or near the Dirac point.17 Motivated by the rapidly evolving science of 3D TI6,7 and the perspective of TI-based spintronic devices,4,17 in this work we carry out an atomistic analysis of the helical states of Bi2Se3. We qualitatively and quantitatively compare our results to the corresponding experimental data on several very important properties including the position of the Dirac point in the surface band structure, the out-of-plane tilting angle of the electron spin in the helical states at EF, and the surface conductance versus electron energy. Our calculation predicts that the out-of-plane tilting angle oscillates in the surface Brillouin zone (BZ) with a 3-fold rotational symmetry, and the helical states penetrate into the bulk which we characterize quantitatively. Figure 1a plots the crystal structure of a slab Bi2Se3 in our analysis. The slab contains six quintuple layers (QLs) and terminates by a selenium atomic layer on both surfaces of the slab. A vacuum of 20 Å thickness is added to the slab to form the supercell which is periodically extended to form a Bi2Se3 film. The atomic structure is fully relaxed by the full-potential linearly augmented plane waves method as implemented in the WIEN2k electronic package.18,19 The atomic structure is considered fully optimized if the force acting on each atom is smaller than 1 mRy/ bohr. The optimized lattice vectors of the bulk structure are found to be a = 4.107 Å and c = 27.89 Å. These compare well with the experimental values of 4.138 and 28.64 Å, respectively.15,21 Received: February 18, 2011 Revised: March 31, 2011 Published: April 15, 2011 2088
dx.doi.org/10.1021/nl200584f | Nano Lett. 2011, 11, 2088–2091
Nano Letters
Figure 1. (a) The supercell of the slab structure of Bi2Se3 with six quintuple layers. A vacuum (the empty region) of 20 Å thickness is added to the slab to form the supercell in the calculation. (b) Calculated band structure of the Bi2Se3 slab shown in (a): the red solid line is with SOI; the blue dotted line is without SOI. The horizontal solid line marks the Fermi level EF of the system, and the dashed line closely passes the Dirac point.
Using the relaxed atomic structure, we calculate the electronic and quantum transport properties including the surface helical states, by carrying out DFT within the Keldysh nonequilibrium Green’s functions (NEGF) using the transport package nanodcal.22,20 Here, a linear combination of atomic orbital basis (LCAO) is used to expand physical quantities, the standard norm-conserving pseudopotentials23 are used to define the atomic cores, the spinorbit coupling is handled at the atomic level,24,25 and the noncollinear spin capability of nanodcal22 is ideal for analyzing the helical spin states. To self-consistently calculate SOI of Bi2Se3 accurately, a triple-ζ polarization LCAO basis (TZDP) is necessary for all the atoms. A k-point mesh of 20 20 1 is employed for k-sampling, and the energy cutoff for the real space grid is taken as 600 Ry. Self-consistent calculations are considered converged when every element of the density matrix is converged to less than 105 a.u. The local spin density approximation (LSDA)26 is used for the exchange-correlation potential in all the calculations.27 The calculated band structure of the Bi2Se3 slab (film) is shown in Figure 1b with (red solid line) and without SOI (blue dotted line). A Dirac cone is clearly formed around the Γ point of the surface Brillouin zone due to a strong SOI of Bi2Se3, which is one of the most important characteristics of 3D TI. Without SOI, the Dirac cone disappears as shown by the blue dotted line. For the slab, the calculated Fermi level EF is ∼72 meV above the Dirac point. Experimentally, the most recent data28 on Bi2Se3 film grown on sapphire gives EF to be 135 meV above the Dirac point. The calculated EF is qualitatively consistent with the experimental result, the quantitative difference is possibly due to residual impurity doping effects in experimental systems and the approximation of the LSDA exchange-correlation functional in theoretical calculations.27
LETTER
Figure 2. Properties of the surface helical spin states of the Bi2Se3 slab at the Fermi energy EF. (a, b) Calculated spin-momentum locking of the helical states on the top (a) and bottom (b) surfaces of the slab. The hexagon is the first Brillouin zone (BZ) of the surface. Γ,K,M are the high symmetry points in the first BZ. The tail of each arrow indicates a Bloch state, the head of arrow indicates its spin direction. (c) The tilting angle Φ versus polar angle θ for helical states on top (red cross) and bottom (blue circle) surfaces. (d) Conductance G versus electron energy E for the Bi2Se3 slab around the Fermi energy (shifted to EF = 0). Unit of the conductance is G0 = e2/h where e is the electron charge and h Planck’s constant. Inset of (d) defines angles Φ,Ψ which are used to characterize the spin directions (S) going around the surface BZ (indicated by the polar angle θ).
Due to SOI, the electronic structure of the Bi2Se3 slab must be analyzed with general spin—as opposed to collinar spin, for which the wave functions are two-component spinors, namely, |ψnkæ = (|ψvnkæ, |ψVnkæ), where n is the band index and k the Bloch wave vector index. We calculate ψnk using the NEGF-DFT method22,20 and construct the k-space density matrix 0 1 jψvnk æÆψvnk j, jψvnk æÆψVnk j A ð1Þ ^nk ¼ jψnk æÆψnk j ¼ @ V F jψnk æÆψvnk j, jψVnk æÆψVnk j From nk, charge contribution from each atomic orbital can be deduced 0 1 1 @ Fvvnk Sk þ Sk Fvvnk , FvVnk Sk þ Sk FvVnk A ð2Þ Qnk ¼ Tr FVvnk Sk þ Sk FVvnk , FVVnk Sk þ Sk FVVnk 2 where S is the overlap matrix of atomic orbitals and Tr stands for trace. The charge contribution can be further decomposed with the help of Pauli matrix σB Qnk ¼ Q I þ Px σx þ Py σ y þ Pz σz
ð3Þ
where Q is the charge of the Bloch state labeled by (n,k), I the unit matrix, and P = (Px,Py,Pz) is the spin polarization vector. The azimuth angle Ψ and tilting angle Φ of the spin polarization P of the state (n,k) can be calculated to show the spin-momentum locking of the helical states. Finally, from the real space density 2089
dx.doi.org/10.1021/nl200584f |Nano Lett. 2011, 11, 2088–2091
Nano Letters
LETTER
matrix we obtain real space charge distribution which can be decomposed by Pauli matrices, then we obtain the spin polarization vector in the real space grids 0 1 Fvvnk ðrÞ FvVnk ðrÞ A Fnk jræ ¼ @ Vv ð4Þ Fnk ðrÞ ¼ Ærj^ Fnk ðrÞ FVVnk ðrÞ
There are two identical surfaces (top and bottom) on the Bi2Se3 slab; hence two degenerate states are found at each k-point. We calculated Bloch states ψnk at EF (see Figure 1b) for many values of k = (kx,ky), by projecting out the spin polarization vector P as discussed above, the spin-momentum locking is directly investigated. The energy EF cuts a circle in the two-dimensional Brillouin zone surrounding the Γ point, giving the allowed k = (kx,ky) values on the circle for the Bloch states. For every allowed k value on the circle, we found that the spin polarization vector P of that Bloch state is always along the tangent of the circle to a very high degree. By calculating many Bloch states this way, the direction of P is plotted in panels a and b of Figure 2 for the top and bottom surfaces, respectively. The bottom surface states are clockwise in k while the top surface states are counterclockwise. It can be seen clearly that the spinpolarization of the states has opposite direction at þk and k on the 2d Brillouin zone, namely, P(k) = P(k). These states at the Fermi level EF are the helical surface states which are the hallmark of topological insulators. Interestingly and quantitatively, at EF, the spin polarization vector P is found to not completely lie in the 2D BZ plan. In the inset of Figure 2d we use polar angle θ which runs around the circle to distinguish different Bloch states, use tilting angle Φ to indicate how much the spin S is pointing out of the 2D plane, and use the azimuth angle Ψ which is between the spin direction along tangent of the circle and k, to quantify the spin-momentum locking. If spin S of the helical state is completely in the 2D plane, Φ should be identically zero for all θ. Figure 2c plots Φ versus θ for the helical states on the top (red cross) and bottom (blue circle) surfaces of the slab: Φ 6¼ 0 and oscillates periodically with a 3-fold symmetry which reflects the rotational symmetry of the Bi2Se3 slab (space group D3d5). Therefore, as the electron momentum of the helical states traces the circle at EF, their spin—while almost lying in the 2D plane of the circle, wobbles from above to below the plane with a period of 120 going around the circle. The maximum value of |Φ| is found to be 4. Experimentally, a tilting angle of 5 was reported9 for TIs belonging to the Bi-family. We have also calculated the azimuth angle Ψ of spin direction along the circle for both surfaces and found Ψ = 90 ( 0.5 for any k, namely, the spin-momentum locking of the helical states is very well established. Quantitatively, the calculated value of the spin polarization at all k on EF, (Fv FV)/(Fv þ FV), is ∼0.86, namely, the helical states are highly polarized in k-space even though it is not 100% polarized, consistent with that reported previously.17 It is important to point out that the total spin, obtained by integrating over the entire slab and entire k, is identically zero. This is expected since SOI does not break time reversal symmetry thus there can be no spontaneous magnetism in this system.29 The calculated conductance20 G of the Bi2Se3 slab versus electron energy E is shown in Figure 2d. It can be seen that G = G(E) exhibits a “V-shape” whose minimum is at the Fermi level. For E < EF, G increases rapidly as E is reduced; this is actually
Figure 3. The charge distribution of the surface helical states calculated at EF and θ = 90. Clearly, the 2-fold degenerate (two surfaces) helical states are spatially rather localized to the individual surfaces: (a) for top surface and (b) bottom surface. The white arrows show the calculated spin-polarization direction of the helical states projected onto each atom.
expected from the band structure in Figure 1b: there are more bulk bands as E is lowered from the EF. As E is increased from EF to 250 meV, conductance increases linearly which reflects the nearly linear dispersion of the Dirac-like Hamiltonian in this range of energy. The increase of G is due to the increase of the circular BZ as E is increased from EF. When E reaches and goes beyond 250 meV, more bands contribute to G as shown in Figure 1b and G increases rapidly. The results are qualitatively very similar to the measured tunneling conductance by scanning tunneling microscopy.13 It is very interesting look at the surface helical states in real space obtained by eq 4. Figure 3 shows the calculated charge density of the helical states for the top and bottom surfaces at Fermi energy EF and at θ = 90 (see Figure 2d). Different colors indicate total contribution to the helical state by atoms in real space. Clearly, at EF, the helical states living in the top and bottom surfaces are largely separated spatially, and they decay into the Bi2Se3 slab. We found that at least five QLs are needed to spatially separate the two degenerate surface helical states living on the top and bottom surfaces in the slab calculation.29 In Figure 3, the white arrows show the spin polarization (∼Fv FV) projected on each atom. The first QL contributes most significantly to the charge and spin of the helical states. In particular, the main contribution to the spin is found from the Bi atoms in the first QL. Inner QLs in the slab contribute weakly to both charge and spin; thus the projected spin is too small to show in Figure 3. Finally, the projection shows that the top and bottom surface helical states have opposite spin direction, as expected from the spin-momentum locking results of Figure 2a,b. Quantitatively, the charge densities of the helical states on both surfaces are contributed 70% by the first QL, 22.3% by the second QL, and 5.6%, 1.4%, 0.45%, and 0.25% by the other inner QLs, respectively. In summary, for the slab model of Bi2Se3, we have calculated properties of the surface helical states from atomic principles selfconsistently without any phenomenological parameter. The calculated band structure has a band inversion and a Dirac cone at the Γ-point of the BZ due to SOI. While the spin-momentum locking is confirmed, the spins of the helical states tilt out of the surface plane by a maximum angle of about 4. The tilting angle oscillates with a 3-fold symmetry going around the 2D BZ, which is commensurate with the rotational symmetry of the crystal. This result should be experimentally testable by spin-resolved ARPES measurements. The topological helical states penetrate into the Bi2Se3 slab and are contributed mainly by atoms in the first quintuple layer. The calculated conductance of the helical states versus electron energy has a V-shaped spectrum, it increases slowly but linearly as a function of energy above EF due to the Dirac dispersion, and becomes rapidly increasing when bulk bands start to contribute. 2090
dx.doi.org/10.1021/nl200584f |Nano Lett. 2011, 11, 2088–2091
Nano Letters
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected];
[email protected].
’ ACKNOWLEDGMENT H.G. thanks Professor Zhong Fang for a communication concerning the thickness needed in the slab model. H.G. is financially supported by NSERC of Canada, FQRNT of Quebec and CIFAR. Y.Z. is supported by National Natural Science Foundation of China (Grant No. 11047104). We gratefully acknowledge RQCHP and CLUMEQ for providing computation facilities. ’ REFERENCES (1) Kane, C. L.; Mele, E. J. Phys. Rev. Lett. 2005, 95, 146802. (2) Bernevig, B. A.; Hughes, T. L.; Zhang, S. C. Science 2006, 314, 1757. (3) Dang, W.; Peng, H.; Li, H.; Wang, P.; Liu, Z. F. Nano Lett. 2010, 10, 2870. Kong, D.; Dang, W.; Cha, J. J.; Li, H.; Meister, S.; Peng, H.; Liu, Z.; Cui, Y. Nano Lett. 2010, 10, 2245. Zuev, Y. M.; Lee, J. S.; Galloy, C.; Park, H.; Kim, P. Nano Lett. 2010, 10, 3037. (4) Kong, D.; Randel, J. C.; Peng, H.; Cha, J. J.; Meister, S.; Lai, K.; Chen, Y.; Shen, Z.-X.; Manoharan, H. C.; Cui, Y. Nano Lett. 2010, 10, 329. Hong, S. S.; Kundhikanjana, K.; Cha, J. J.; Lai, K.; Kong, D.; Meister, S.; Kelly, M. A.; Shen, Z.-X.; Cui, Y. Nano Lett. 2010, 10, 3118. Liu, H. W.; Yuan, H. T.; Fukui, N.; Zhang, L.; Jia, J. F.; Iwasa, Y.; Chen, M. W.; Hashizume, T.; Sakurai, T.; Xue, Q. K. Cryst. Growth Des. 2010, 10, 4491. (5) Teweldebrhan, D.; Goyal, V.; Balandin, A. A. Nano Lett 2010, 10, 1209. (6) Hasan, M. Z.; Kane, C. L. Rev. Mod. Phys. 2010, 82, 3045. (7) Qi, X. L.; Zhang, S. C. Phys. Today 2010, 63, 33. (8) Xia, Y.; Qian, D.; Hsieh, D.; Wray, L.; Pal, A.; Lin, H.; Bansil, A.; Grauer, D.; Hor, Y. S.; Cava, R. J.; Hasan, M. Z. Nat. Phys. 2009, 5, 398–402. (9) Hsieh, D.; Xia, Y.; Qian, D.; Wray, L.; Dil, J. H.; Meier, F.; Osterwalder, J.; Patthey, L.; Checkelsky, J. G.; Ong, N. P.; Fedorov, A. V.; Lin, H.; Bansil, A.; Grauer, D.; Hor, Y. S.; Cava, R. J.; Hasan, M. Z. Nature 2009, 460, 1101–1106. (10) Zhang, Y.; He, K.; Chang, C. Z.; Song, C. L.; Wang, L. L.; Chen, X.; Jia, J. F.; Fang, Z.; Dai, X.; Shan, W. Y.; Shen, S. Q.; Niu, Q.; Qi, X. L.; Zhang, S. C.; Ma, X. C.; Xue, Q. K. Nat. Phys. 2010, 6, 584. (11) Fu, L.; Kane, C. L. Phys. Rev. B 2007, 76, 045302. (12) Seo, J.; Roushan, P.; Beidenkopf, H.; Hor, Y. S.; Cava, R. J.; Yazdani, A. Nature 2010, 466, 343. (13) Hanaguri, T.; Igarashi, K.; Kawamura, M.; Takagi, H.; Sasagawa, T. Phys. Rev. B 2010, 82, 081305. (14) Zhang, H.; Liu, C. X.; Qi, X. L.; Dai, X.; Fang, Z.; Zhang, S. C. Nat. Phys. 2009, 5, 438. (15) Zhang, W.; Yu, R.; Zhang, H. J.; Dai, X.; Fang, Z. New J. Phys. 2010, 12, 065013. (16) Liu, C. X.; Zhang, H. J.; Yan, B.; Qi, X. L.; Frauenheim, T.; Dai, X.; Fang, Z.; Zhang, S. C. Phys. Rev. B 2010, 81, 041307(R). (17) Yazyev, O. V.; Moore, J. E.; Louie, S. G. Phys. Rev. Lett. 2010, 105, 266806. (18) Blaha, P.; Schwarz, K.; Sorantin, P.; Trickey, S.B. Comput. Phys. Commun. 1990, 59, 399–415. (19) In the structure relaxation using the WIEN2k electronic package, the radii of the muffin-tin spheres are chosen to be 2.4 bohr. A plane-wave cutoff of Rmt Kmax has been chosen as 8.0 and the angular expansion is made up by 10 spheric harmonic functions in the muffin tins. (20) Taylor, J.; Guo, H.; Wang, J. Phys. Rev. B 2001, 63, 245407. Maassen, J.; Ji, W.; Hong, Guo Nano Lett. 2011, 11, 151.
LETTER
(21) Wyckoff, R. W. G. Crystal Structures, 2nd ed.; Interscience Publishers: New York, 1964. (22) For more details of the NEGF-DFT method used in the nanodcal tranpsort package, see ref 20; see also http://www.nanoacademic.ca. In this work, we focus on equilibrium tranpsort property at zero external bias; hence NEGF is reduced to the usual equilibrium Green’s functions. (23) Kleinman, L.; Bylander, D. M. Phys. Rev. Lett. 1982, 48, 1425. (24) Theurich, G.; Hill, N. A. Phys. Rev. B 2001, 64, 073106. (25) Fernandez-Seivane, L.; Oliveria, M. A.; Sanvito, S.; Ferrer, J. J. Phys: Condens. Matter 2006, 18, 7999. (26) Perdew, J. P.; Wang, Y. Phys.Rev. B 1992, 45, 13244. (27) For bulk Bi2Se3, nanodcal20,22 gives a band gap 0.19 eV which is close to the 0.22 eV obtained by WIEN2k. These LSDA bulk values are smaller than the experimental gap of 0.35 eV; see Table 3 in:Black, J.; Conwell, E.M.; Seigle, L.; Spencer, C.W. J. Phys. Chem. Solids 1957, 2, 240. From Figure 1b, the 6 QL slab has a ‘‘bulk band gap’’ of 0.26 eV. We have also calculated 3QL and 2QL; their ‘‘bulk gaps’’ are 0.44 and 1.2 eV, respectively. The gaps thus shrink toward the bulk value as the slab becomes thicker. (28) Chang, C. Z.; He, K.; Liu, M. H.; Zhang, Z. C.; Chen, X.; Wang, L. L.; Ma, X. C.; Wang, Y. Y.; Xue, Q. K. Preprint http://arxiv.org/abs/ 1012.5716. (29) In Figure 2c the tilting angle Φ(θ) does not average to zero when θ goes from zero to 360, indicating a small þSz value on the bottom surface and an equal but negative Sz value on the top surface. The total Sz of the slab is still identically zero. The nonvanishing average Φ on the two surfaces is an indication that even at the thickness of six QLs, there is still a very weak but nonzero interaction between the two surfaces across the slab. Investigating even thicker slabs by DFT requires prohibitively large computer resources beyond our ability.
2091
dx.doi.org/10.1021/nl200584f |Nano Lett. 2011, 11, 2088–2091
