Topological Protected Dirac Cones in Compressed Bulk Black Phosphorus
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n Topological Protected Dirac Cones in Compressed Bulk Black Phosphorus
Ruixiang Fei, Vy Tran and Li Yang Department of Physics, Washington University in St. Louis, St. Louis, MO 63130, USA (Dated:)Using the k · p theory and first-principles simulations, we report that applying a moderate pressure( > PACS numbers:
To date Dirac materials, such as graphene [1, 2],topological insulators (TIs) [3–5] or spin Hall insulators(SHIs) [6, 7], topological crystalline insulators (TCIs)[8, 9] and Weyl semimetals (WSs) [10, 11], have drawntremendous research interest. Graphene hosts two-dimensional (2D) Dirac fermions [12, 13]; TIs are mate-rials with a bulk energy gap but possess gapless 2D (1D)Dirac surface states that are protected by time-reversalsymmetry [14, 15]; TCIs exhibit metallic surface statesprotected by the mirror symmetry of the crystal [9]; WSs,which are protected by crystal symmetries, exhibit en-ergy overlaps that occur only at a set of isolated pointsin momentum space where the linearly dispersed bandscross at the Fermi level [10]. The impact and implicationsof these materials are difficult to overstate. The discov-ery of the anomalous Hall effect [16, 17], the fractionalHall effect [18], chiral anomaly [19, 20], and Majoranafermions [21, 22], to name a few, have inspired enormousefforts in search for new Dirac materials.Recently a new type of 2D semiconductor, few-layerblack phosphorus (BP) (phosphorene), has been fabri-cated [23–25]. It exhibits promising carrier mobilitiesand unique anisotropic transport, optical, and thermalproperties [23–29]. In particular, few-layer BP posses aunique band structure, whose dispersion is nearly linearalong the armchair direction but parabolic along otherdirections [30, 31]. An obvious idea follows from thisunique structure: if one can significantly reduce the bandgap or even achieve the band inversion, novel featuresmay result, e.g., the formation of Dirac cones. However,the band gap of monolayer and few-layer BP is signifi-cant (up to 2.2 eV) [26, 29, 32], making it impractical toclose. Alternatively, bulk black phosphorus has a muchsmaller band gap, only around 300 meV [26, 33], whichis more tractable.In this work, we show that the k · p Hamiltonian ofbulk BP can be written as the free-electron-gas Hamil-tonian with off-diagonal terms for describing interbandinteractions, which are of the same form as the Rasshbaand Dressehaus spin-orbit coupling (SOC). Thus we callthis off-diagonal interaction the ”pseudo spin-orbit” cou-
FIG. 1: (Color online) (a) Top view of the atomic structureof bulk BP. The unit cell and lattice vectors are marked. Thezigzag and armchair directions are illustrated as well. (b) Sideview of BP. (c)The first BZ. pling (PSOC). If the band gap is diminished and the bandinversion is realized, PSOC tends to open a gap in theenergy spectrum. However, similar to TIs, the band sym-metry requires band connections. Competition betweenPSOC and band symmetry causes the Fermi surface toundergo a topological transition from a sphere to a ringof 2D Dirac cones. Using first-principles simulations, weshow that pressure can be an efficient tool to producethese 2D Dirac cones and even 1D Dirac cones, whilealso tuning the characteristic Fermi velocity. Interest-ingly, because our predicted Dirac cone feature in com-pressed BP does not require the time reversal symmetry,it is not destroyed by magnetic fields. Thus we providea means for detecting this unique electronic structure bydemonstrating the unusual energy scaling law exhibitedby Landau levels.The atomic structure of bulk BP is plotted in Figures1 (a) and (b). Here we use a simple orthorhombic con-ventional unit cell, which is twice the size of the primi-tive cell. The corresponding first Brillouin Zone (BZ) isshown in Figure 1 (c). The double-cell procedure sim-
FIG. 2: (Color online)(a) and (b) are the band structures ofintrinsic BP. The red dots are from the HSE06 simulation andthe blue line is from the k · p calculations. (c) and (d) are theband structures calculated from the k · p calculations underthe band-inversion condition. The Fermi level is set to bezero. ply folds the first BZ and does not affect any essentialphysics; this is convenient for presenting the topologyof the Dirac cones. The bulk BP atomic structures isfully relaxed according to the forces and stresses calcu-lated using density functional theory (DFT) with thePerdew, Burke and Ernzerhof (PBE) functional associ-ated with van der Waals corrections [34, 35]. BecauseDFT is known to underestimate the band gap, we usea hybrid functional theory (HFT), HSE06 [36, 37], tocalculate electronic band structures, which is in goodagreement with the GW calculations and experimentalmeasurements [26, 33] of bulk BP.We begin with the widely used k · p model to capturethe essential BP’s electronic structure. The D h pointgroup invariance allows us to describe BP using a two-band model similar to the one used for monolayer BP[38, 39]. Thus the low-energy effective Hamiltonian isgiven by H = (cid:18) E c + αk x + βk y + γk z v f k y v f k y E v − λk x − µk y − νk z , (cid:19) (1)where k x , k y , and k z are components of the crystal mo-mentum along the x , y , and z directions, respectively.The conduction band and the valance band energies, E c = 0 . eV and E v = − . eV , are set to fit thebulk band gap (300 meV). α , β , γ , λ , µ , and ν are fitted FIG. 3: (Color online) The 3D band structures of BP of theX-Y plane (a), the Y-Z plane (b), and the X-Z plane (c). Theband overlap in (c) forms a ring structure. (d) The evolutionof the band structure from the X-Z plane to the X-Y planevia changing of the tilted angle. parameters that describe the band edge curvature, and v f is the Fermi velocity of the nearly linear band disper-sion along the armchair ( y ) direction. By appropriatelychoosing these parameters, as shown in Figures 2 (a) and(b), the k · p result nearly perfectly matches that of thefirst-principles HSE06 calculations of intrinsic BP.The special characteristics of Eq. 1 are the off-diagonalterms of the effective Hamiltonian, which produce theunique linear dispersion along the armchair ( y ) direction,as shown in Figure 2 (a). In fact, This off-diagonal inter-action can be split into two terms H off = H DSOC + H RSOC = 1 / v f ( k x σ y + k y σ x ) − / v f ( k x σ y − k y σ x ) , (2)where σ x and σ y are Pauli matrices. Interestingly, thesetwo terms ( H DSOC and H RSOC ) have the forms of theDressehaus and Rashba SOC [40], respectively. Realis-tically, Dressehaus and Rashba SOC cannot exist in in-trinsic bulk BP, because its lattice structure possessesthe inversion symmetry. Hence we call these interactionsPSOC. Because of their similar formalisms, we expectequivalent ”SOC effects” to appear in bulk BP. In par-ticular, considering the dramatic effect of SOC in TIs,it is natural to expect that diminishing or inverting theband gap in bulk BP will produce novel electronic char-acteristics.It is easy to mimic the band inversion in the k · p model; one simply assigns the valence band edge energy E v a higher value than the conduction band edge energy (d)(c) (b) Pressure (x103 atmosphere) F e r m i v e l o c i t y ( x m / s ) Compress (%)
Armchair Hole/Electron Zigzag Hole Zigzag Electron (a) E ne r g y ( e V ) a a a a -0.2 -0.1 0.0 0.1 0.201234 Energy (eV) D en s i t y o f S t a t e ( D O S ) X Y
X Y
FIG. 4: (Color online)(a) The HSE06-calculated band struc-ture of BP under the critical pressure (0.6 GPa). (b) TheHSE06-calculated band structure of BP under the 0.9 GPapressure. (c) The Anisotropic fermi velocity of electrons andholes of 2D Dirac cones varied with the pressure. E c in Eq. 1. Surprisingly, instead of producing simpleband overlaps, we observe an unusual topological transi-tion at the Fermi surface. In Figures 2 (c) and (d), theband crossing at the Γ point is opened by the PSOC,while the band crossing is preserved at the K point inthe Γ − X direction. As a result, 2D graphene-like Diraccones are formed and the Schematic 3D plots of the elec-tronic structure are presented in Figures 3 (a), (b), and(c) for three specific cutting planes, X-Y, Y-Z and X-Z,respectively. The linear band dispersions and Dirac conesare observed in both X-Y and Y-Z planes and a ring-likeband crossing occurs in the X-Z plane. Actually, we canregard the Dirac cones observed in Figure 3 (a) and (b)are those points marked in the ring structures of Figure3 (c). In this sense, any point on the ring of Figure 3(c) corresponds to a Dirac cone in other 2D planes. InFigure 3 (d), we show that this ring-like band crossingcontinuously evolves to a pair of Dirac cones by varyingthe tilting angles from the X-Z plane to the Y-Z plane.In other words, the set of tilting angles that do not ex-hibit these Dirac cones is of measure zero, and thus ourfindings are robust.The mechanisms responsible for forming the ring ofDirac cones are similar to those in TIs, but there are somecrucial differences. First, similar to the band inversionand gap-opening that occurs in TIs due to SOC interac- tions, the gap-opening at the Γ point of Figure 2 (c) fromthe PSOC of states residing on the linear band along theΓ − Y direction. Second, the gap of BP is not completelyopened; the valence and conduction bands still touch atthe K point (Figure 2 (c)) or by the ring (Figure 3 (c)).This gapless feature is guaranteed by a symmetry con-straint based on two facts: the bottom conduction bandis formed by what was originally top valence band priorto the band inversion; the parity symmetries of the topof the valence band and the bottom of the conductionband are inverse to each other, which is evidenced by thebright dipole optical transitions [26]. Therefore, as shownin Figures 2 (c) and 3 (c), contact between the conduc-tion and valence band is required in order to provide acontinuous symmetry evolution. This mechanism is sim-ilar to the one responsible for forming gapless surfacestates in TIs [41]. However, there are no surface statesin our calculated bulk BP, thus it is the bulk valence andconduction bands that must touch. As a result, thesetwo factors force the formation Dirac cones. Moreover,as proposed by a recent work [42], further breaking theinversion symmetry and including the spin-orbit couplingin few-layer BP will open the bulk band gap and producesurface topological states. Finally, another special char-acteristic of these bulk Dirac cone states in BP is thatthey are robust even in the face of broken time reversalsymmetry. Therefore, these Dirac cones can uniquely beobserved in experiments that employ magnetic fields andare potentially useful for magnetic devices.Following the above analysis of the k · p model and be-yond 2D Dirac cones, 3D Weyl points could be achievedif an extra direction, such as the x or z , also exhibits theoff-diagonal linear interband interaction (PSOC) term.Therefore, one may search for WSs by looking for re-alizing the band inversion of the materials that exhibitnearly linear dispersions along two orthogonal directions.Beyond the k · p model, an obvious question is howone can practically reduce the band gap or achieve bandinversion in BP. Towards this end, we turn to the first-principles simulations for quantitative answers. OurHSE06 calculations show that applying an external pres-sure along the y (armchair) or x (zigzag) direction pro-duces promising effects. Note that, BP is much softeralong the y (armchair) direction because of its anisotropicatomic structure [43]. We find that applying a pressureof 0.75 GPa can shrink the lattice constant along the y direction by 4% but only by 1% and 2% for the x and z directions, respectively. Therefore, we will consider thechanges made by applying pressure along the y direc-tion because it can change the atomic structure and cor-responding electronic structure more significantly via arelatively small pressure.First, as shown in Figure 4 (a), an applied critical pres-sure (0.6 GPa) diminishes the band gap and a 1D Diraccone is produced along the Γ − X direction. This pres-sure is well within the capability of current high-pressureexperiments [44, 45].Second, larger pressures create a band inversion. Ex- nn E ne r g y ( e V ) n Intrinsic B// Z 1T 5T 10T 2D Dirac cone B// Z 1T 5T 10T 1D Dirac cone B// Z 1T 5T 10T
FIG. 5: (Color online) (a) The energy spectrum of Landaulevels (index n) of intrinsic BP with magnetic field; (b) thatunder the critical pressure (0.6 GPa), in which the 1D Diraccone is formed; (c) That under the pressure of 0.9 GPa, inwhich 2D Dirac cones are formed. citingly, similar to the k · p theory predictions, a pressureof 0.9 GPa yields a band crossing at the Γ point and apair of 2D Dirac cones are formed, as shown in Figure 4(b). The corresponding density of states (DOS) is shownin Figure 4 (c). We can clearly see that a linear DOSin the low-energy regime, which is similar to graphene.From these first-principles results, we also find that thatthe Fermi velocities of compressed BP is smaller thanthat of graphene and that these generated Dirac conesare highly anisotropic in momentum space. Obviously,the Fermi velocities must be anisotropic as well and canbe efficiently tuned by varying the applied pressure, asconcluded in Figure 4 (d).We propose that the evolution of the electronic struc-ture from conventional semiconductors to one possessing1D Dirac cones and 2D Dirac cones can be observed bymeasuring the energy spectrum of Landau levels; thismethodology is widely used for identifying Dirac bandstructures [2]. For example, considering the X-Y plane,before reaching the critical point for closing the bandgap, the energy spectrum of the Landau levels is that ofa typical semiconductor, i.e., it is linear with the levelindex, as shown in Figure 5 (a). Here the spectrum isgiven by E = E c + ( n + 1 / hω c E = E v − ( n + 1 / hω v , (3)where the cyclotron frequencies of electrons and holes are ω c = eB/ √ m cx m cy and ω v = eB/ √ m vx m vy . m cx , m cy , m vx and m vy are effective mass of electrons and holesalong the X and Y directions, respectively. Importantly, there is no zero-energy solution for this case.At the critical pressure, a perfectly 1D linear band dis-persion along the Γ − Y direction simplifies Eq. 1 by set-ting the parabolic term β and µ to be zero. The effectiveHamiltonian with the magnetic field introduced by theLandau gauge −→ A = ( − By, ,
0) is12 m c (¯ hk x − eBy ) φ A − i ¯ hv f ddy φ B = Eφ A , − m v (¯ hk x − eBy ) φ B − i ¯ hv f ddy φ A = Eφ B , (4)where v f = 3 . m/s is the Fermi velocity of the 1Ddirac cone obtained from HSE06 simulations. The corre-sponding energy spectrum of Landau Levels in magneticfield is presented in Figure 5 (b), which is deviated fromthe linear relation of Figure 5 (a). Finally, for higherpressures that form 2D Dirac cones, the Landau levelsfollows the well-known square-root dependence [2, 12], E n ∝ √ n , and exhibits a unique zero-energy solution,as shown in Figure 5 (c).In conclusion, we predict that applying a moderatepressure ( > [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, andA. A. Firsov. Nature , 197-200 (2005).[2] Y. Zhang, Y. W. Tan, H. L. Stormer, and Philip Kim.Nature , 201-204 (2005).[3] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K.Mo, X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang,S. C. Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen. Science , 178-181(2009).[4] D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, F. Meier,J. Osterwalder, L. Patthey, J. G. Checkelsky, N. P. Ong,A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor,R. J. Cava, and M. Z. Hasan. Nature , 1101-1105(2009).[5] L. Fu, C. L. Kane, and E. J. Mele. Phys. Rev. Lett. ,106803(2007). [6] S. Murakami, N. Nagaosa, and S. C. Zhang. Phys. Rev.Lett. , 156804 (2004).[7] M. Knig, S. Wiedmann, C. Brne, A. Roth, H. Buhmann,L. W. Molenkamp, X. L. Qi, and S. C. Zhang. Science , 106802(2011).[9] Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T.Takahashi, K. Segawa, and Yoichi Ando. Nat. Phys. ,800-803 (2012).[10] S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane,E. J. Mele, and A. M. Rappe. Phys. Rev. Lett. , 864-867(2014).[12] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim. Reviews of Modern Physics , 109 (2009).[13] S. Das Sarma, Shaffique Adam, E. H. Hwang, and EnricoRossi. Reviews of Modern Physics , 407(2011).[14] M. Z. Hasan, and C. L. Kane, Reviews of Modern Physics , 3045(2010).[15] X. L. Qi, and S. C. Zhang, Reviews of Modern Physics , 1057(2011).[16] C. Z. Chang, J.Zhang, X. Feng, J. Shen, Z. Zhang, M.Guo, K. Li, Y. Ou, P. Wei, L. L. Wang, Z. Q. Ji, Y. Feng,S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S. C. Zhang, K.He, Y. Wang, L. Lu, X. Ma, and Q. K. Xue, Science ,167-170(2013).[17] W. K. Tse, Z. Qiao, Y. Yao, A. H. MacDonald, and QianNiu. Phys.Rev. B. , 155447(2011) .[18] X. Du, I. Skachko, F. Duerr, A. Luican, and E. Y. Andrei.Nature , 192-195(2009).[19] A. A. Burkov. Phys. Rev. Lett. , 247203(2014).[20] A. A. Zyuzin, and A. A. Burkov. Phys. Rev. Lett. ,115133(2012).[21] L. Fu, and C. L. Kane. Phys. Rev. Lett. ,096407(2008).[22] A. R. Akhmerov, J. Nilsson, and C. W. J. Beenakker.Phys. Rev. Lett. , 372(2014).[24] H. Liu , A. T. Neal , Z. Zhu , Z. Luo , X. Xu , D. Tomnek, and P. D. Ye , ACS Nano , 4033 (2014).[25] F. Xia, H. Wang, and Y. Jia. Nat. Commun. ,5458(2014).[26] V. Tran, R. Soklaski, Y. Liang, and L. Yang. Phys. Rev.B , 235319 (2014).[27] Ruixiang Fei, and L. Yang. Nano Lett. , 2884 (2014).[28] Ruixiang Fei, A. Faghaninia, R. Soklaski, J. Yan, C. Lo,and Li Yang. Nano Lett. , 63936399 (2014).[29] X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia,H. Zhao, H. Wang, L. Yang, X. Xu, and F. Xia. arXivpreprint, arXiv:1411.1695 (2014).[30] Vy Tran and Li Yang, Phys. Rev. B 89, 245407 (2014).[31] A. S. Rodin, A. Carvalho, and A. H. Castro Neto. Phys.Rev. Lett. , 176801 (2014).[32] Liangbo Liang, Jun Wang, Wenzhi Lin, Bobby G.Sumpter, Vincent Meunier, and Minghu Pan, Nano Lett. , 6400 (2014).[33] Robert W. Keyes, Phys. Rev. , 580 (1953)[34] J. P. Perdew, K. Burke, and M. Ernzerhof. Phys. Rev.Lett. , 3865(1996).[35] S. Grimme. J. Chem. Phys. , 1463-1473(2004).[36] J. Heyd, G. E. Scuseria, and M. Ernzerhof. J. Chem.Phys. , 8207-8215(2003).[37] J. Heyd, G. E. Scuseria, and M. Ernzerhof. J. Chem.Phys. , 219906(2006).[38] X. Y. Zhou, R. Zhang, J. P. Sun, Y. L. Zou, D. Zhang, W.K. Lou, F. Cheng, G. H. Zhou, F. Zhai, and K. Chang.arXiv preprint arXiv:1411.4275(2014).[39] M. Ezawa. New Journal of Physics , 115004(2014).[40] I. utic, J. Fabian, and S. D. Sarma. Reviews of ModernPhysics , 323(2004).[41] H.Zhang, C. X. Liu, X. L. Qi, X. Dai, Z. Fang, and S. C.Zhang. Nat. Phys. , 438-442 (2009).[42] Q. Liu, X. Zhang, L. B. Abdalla, A. Fazzio, and A.Zunger. arXiv preprint arXiv:1411.3932(2014).[43] Q. Wei, X. Peng. Appl. Phys. Lett. , 251915(2014).[44] Y. Katayama, T. Mizutani, W. Utsumi, O. Shimomura,M. Yamakata, and K. Funakoshi. Nature, , 170-173(2000).[45] G. Monaco, S. Falconi, W. A. Crichton, and M. Mezouar.Phys. Rev. Lett. , 255701(2003).[46] G. Kresse, and J. Furthmller. Phys.Rev. B.54