Disorder effects of vacancies on the electronic transport properties of realistic topological insulators nanoribbons: the case of bismuthene
Armando Pezo, Bruno Focassio, Gabriel R. Schleder, Marcio Costa, Caio Lewenkopf, Adalberto Fazzio
DDisorder effects of vacancies on the electronic transport properties of realistictopological insulators nanoribbons: the case of bismuthene
Armando Pezo,
1, 2, ∗ Bruno Focassio,
1, 2
Gabriel R. Schleder,
1, 2
Marcio Costa,
3, 2
Caio Lewenkopf, and Adalberto Fazzio
2, 1, † Federal University of ABC (UFABC), 09210-580, Santo André, São Paulo, Brazil Brazilian Nanotechnology National Laboratory (LNNano),CNPEM, 13083-970, Campinas, São Paulo, Brazil Instituto de Física, Universidade Federal Fluminense, 24210-346 Niterói, Rio de Janeiro, Brazil (Dated: October 23, 2020)The robustness of topological materials against disorder and defects is presumed but has notbeen demonstrated explicitly in realistic systems. In this work, we use state-of-the-art densityfunctional theory and recursive nonequilibrium Green’s functions methods to study the effect ofdisorder in the electronic transport of long nanoribbons, up to
157 nm , as a function of vacancyconcentration. In narrow nanoribbons, even for small vacancy concentrations, defect-like localizedstates give rise to hybridization between the edge states erasing topological protection and enablingbackscattering events. We show that the topological protection is more robust for wide nanoribbons,but surprisingly it breaks down at moderate structural disorder. Our study helps to establish somebounds on defective bismuthene nanoribbons as promising candidates for spintronic applications.
I. INTRODUCTION
Topological materials have been intensively studied inrecent years [1–5] unveiling interesting new physics andopening new applications possibilities in spin-based elec-tronic devices [6]. Of particular interest are large bandgap topological insulators (TIs), that are good candidatesfor the realization of the quantum spin Hall (QSH) effectat room temperature [7, 8]. In two-dimensional (2D)QSH insulators, the edges of the sample carry metal-lic states that are protected by time-reversal symme-try (TRS) [9, 10] and decay exponentially into the bulk[11–13]. Moreover, theory predicts that in nanoribbonsthese edge states carry dissipationless helical spin cur-rents. The penetration depth, that quantifies the edgestates exponential decay rate, is (roughly) inversely pro-portional to the band gap [12, 14] and plays a key rolein the nanoribbon transport properties. As it happenswith 3D topological insulators [15], to fully display thefeatures of a TI, the nanoribbon width must be muchwider than the penetration depth ξ of the edge states,otherwise they can easily hybridize.Among several candidates for topological materials,bismuthene, also known as buckled or bilayer bis-muthene, is of special interest due to its large electronicband gap of . , along with its structural stabilityand large spin-orbit coupling (SOC) [16–18]. The exis-tence of charge puddles and other types of defects thatcould be detrimental for the formation of a topologicalphase do not play an important role in bismuthene, asexperimentally demonstrated [7]. Indeed, since its ex-perimental realization, it was proposed that even amor-phous structures of bismuthene occurring before the an- ∗ [email protected] † [email protected] nealing process [19], as well as strongly disordered sys-tems [20], support topological states. The robustness ofbismuthene non-trivial topology [19] makes it ideal formaterial design by tuning the lattice constant and effec-tive spin-orbit coupling (SOC) by epitaxial constraint,and lighter elements substitutional alloying with Sb orAs, without causing a topological transition [16, 21]. Thepresence of non-magnetic defects leaves the band topol-ogy unchanged in these materials [18–20].However, a recent study using a schematic tight-binding model [22] has shown that vacancy induced lo-calized states can give rise to local magnetic momentsand destroy the topological protection. Contrariwise, us-ing a different tight-binding toy model, it was numer-ically shown that small concentrations of vacancy de-fects do not eliminate the topological edge states, butcan cause edge state hybridization in certain energies in-tervals within the topological gap [23]. To settle suchkind of controversy and better understand the interplaybetween edge states and defects a systematic study using ab initio methods is in order.Vacancies on bismuthene show small formation en-ergies where sp bismuth dangling bonds can be dis-tributed around the vacancy leading to resonances in theband gap, modifying the electronic properties of the hostmaterial. Theoretical calculations on bismuthene showno evidence of magnetic moments induced from vacan-cies [17], preserving TRS and retaining its non-trivialtopological band structure. Recent ab initio calculationsdemonstrate that transport along the edges is insensitiveto vacancies created in ultra-narrow TI zigzag nanorib-bons [24]. Even though these vacancies allow the devel-opment of magnetic moments, the perfect conductancefor energies within the bulk gap is recovered providedthe vacancies are passivated by hydrogen. Our study, inturn, addresses more realistic system sizes, namely, bothmuch wider and much longer, unveiling a different kind a r X i v : . [ c ond - m a t . m e s - h a ll ] O c t of disordered-induced transport mechanism.In this work, we investigate the electronic transportproperties of buckled bismuthene nanoribbons of realisticsizes using the full Hamiltonian in the orbital representa-tion obtained from density functional theory (DFT) [25–27] calculations combined with recursive nonequilibriumGreen’s functions (NEGFs) [28–31]. First, we explore thebackscattering mechanism due to inter-edge hopping me-diated by vacancy localized states in different nanoribbonwidths and discuss the detrimental effects in transportproperties caused by single vacancies in narrow nanorib-bons. We proceed to investigate the robustness of theedge states for different ribbon widths, reaching extendedlengths (
157 nm ), as a function of vacancy concentration,demonstrating distinct behaviors in the transport prop-erties for bulk and edge states. We show that narrowribbons, whose widths w are comparable with the edgestates penetration depth ξ , present the onset of Ander-son localization effects already at low vacancy concentra-tions. In distinction, for wide nanoribbons, where w (cid:29) ξ ,the edge states maintain a quantized conductance in thetopological gap at low disorder concentrations. Surpris-ingly, we find that topological protection is destroyed al-ready at modest vacancy concentrations. II. COMPUTATIONAL METHODS
Our calculations combine the flexibility of a plane-wavebasis set to obtain the optimized structures with a local-ized basis for electronic transport. In both cases, we usethe Perdew-Burke-Ernzerhof (PBE) [32, 33] exchange-correlation functional. We perform the geometry opti-mizations with a plane-wave basis as implemented in the
VASP package [34, 35]. In these calculations, we employ
400 eV for the plane-wave expansion cutoff. Vacancies aremodeled by removing an atom and performing the geome-try optimization with force criterion of × − eV Å − .The ionic potentials are described using the projectoraugmented-wave (PAW) method [36].The transport calculations use full ab initio DFTHamiltonian matrices obtained directly from the
SIESTA code [37], employing atom-centered single- ζ plus polarization (SZP) basis sets. We use the energy cut-off for real-space mesh of
350 Ry , sampling the reciprocalspace with 10 k -points along the periodic direction ofthe presented nanoribbons, which edges were hydrogen-passivated. We add Å of vacuum in both non-periodicdirections to avoid spurious interactions between periodicimages. The self-consistent SOC is introduced via anon-site approximation [38] using fully-relativistic norm-conserving pseudopotentials [39]. The system Hamilto-nian and overlap matrices are obtained after performinga full self-consistent cycle.The electronic transport calculations are implementedfollowing the standard NEGF approach [28–30]. Weconsider a two-probe terminal setting, as illustrated inFig. 1(a). The left (L) and right (R) electrodes are mod-
L R xy (a)(b) FIG. 1. Schematic representation of the two-probe setup stud-ied in this work. (a) The electrodes are located at the left (L)and right (R) sides of the device. (b) Zoomed-in view of thescattering region (S), illustrating some of its building blockscontaining vacancies. eled by semi-infinite pristine zigzag bismuthene leads.Our approach employs a decimation technique [40–45]that allows us to address large system sizes. The scatter-ing region is partitioned into building blocks connectedby first neighbour interactions. Each building block iscomputed through DFT. The procedure takes into ac-count all degrees of freedom of the system comprised ofall the building blocks. The partition scheme is depictedin Fig. 1(b). Accordingly, the system Hamiltonian is ex-pressed in the localized basis by the block matrix H = H L H C H † C H S H C H † C H R (1)where H L and H R are the Hamiltonian matrices describ-ing the left and right electrodes, while H C is the couplingbetween the leads and the central region, and H S is thetridiagonal block matrix H S = H H C . . . H † C H . . . . . . ... ... . . . ... ... . . . . . . H N − H C . . . H † C H N (2)representing the scattering region (S).The scattering region spin resolved retarded Green’sfunction [30, 46] reads G rS ( E ) = (cid:0) E + S S − H S − Σ rL − Σ rR (cid:1) − (3)where E + = lim δ → + E + iδ , S S is the overlap matrix,and Σ rR/L ( E ) = ( E + S C − H C ) G r ,R/L ( E + S † C − H † C ) arethe embedding self-energies that account for the systemdecay width due to the coupling with the leads. Here G r ,R/L is the retarded surface Green’s function of the R/L electrode [40, 47]. The block structure of H S allowsfor a very efficient computation of G rS ( E ) using the recur-sive Green’s function method (see, for instance, ref. [44],for a review).In the linear response, at small bias, the zero-temperature conductance is given by the Landauer for-mula G ( E F ) = ( e /h ) T ( E F ) , where the transmission T reads [28–30] T ( E ) = Tr [ Γ L ( E ) G aS ( E ) Γ R ( E ) G rS ( E )] (4)where G aS = [ G rS ] † and the decay width matrices Γ L/R are given by Γ L/R = i ( Σ rL/R − Σ aL/R ) . III. RESULTS AND DISCUSSION
In this section we analyze the effect of vacancies onthe transport properties of bismuthene zigzag nanorib-bons. We begin by discussing the topological propertiesof pristine nanoribbons of different representative widths.Next, we investigate the effect of a single-vacancy on theelectronic and transport properties of these systems. Fi-nally, we study the conductance of these systems for dif-ferent vacancy concentrations and nanoribbon widths.Let us first investigate the effect of the ribbon width onpristine systems. In Fig. 2 we show the electronic bandstructure and the conductance for bismuthene nanorib-bons with 3 different representative widths. To help thediscussion, the bulk-topological gap ∆ TG , correspondingto − . < E < . , is indicated in grey. Figure2(a) shows the results for a narrow nanoribbon, w ofwidth Å, comparable to those of ref. [24]. The elec-tronic states bridging the bulk-topological gap ∆ TG aresplit and do not display the Kramer’s degeneracy. Hence,there is no manifestation of the topological edge states.Figure 2(b) shows an intermediate width nanoribbon.namely, w of width Å. Here, there is still a smallgap at the Γ -point, though the edge states are degener-ate as expected for a topological insulator. In addition,the top of the valence and bottom of the conduction triv-ial state bands approach the corresponding bulk energies,indicating that finite width quantization effects are muchsmaller than in the previous case. Finally, Fig. 2(c) showsa wide nanoribbon, w of width Å, with no gap anddegenerate edge states with helical texture, as expectedfor a TI. Here, the bulk-topological gap corresponds veryclosely to the energy interval where one finds only edge (a)(b)(c)
FIG. 2. Electronic band structure along the Γ − X direction(left panels) and the corresponding conductance G (right pan-els) for pristine bismuthene nanoribbons of widths (a) w =20 Å, (b) w = 65 Å, and (c) w = 110 Å. The gray energywindows correspond to the bulk topological gap. states. In all cases, the conductance G ( E ) is quantizedand the conductance steps are observed as expected forpristine systems [30].The lack of topological protection in narrow ribbons isa result of strong overlap between the states at the op-posite system edges. This can be understood in terms ofthe spatial localization of edge states or, more precisely,their so-called penetration length ξ . The latter can beroughly estimated from the mass term in the k · p DiracHamiltonian describing the inverted bulk band gap as[14, 48–50] ξ ≈ (cid:126) v F / ∆ TG , (5)where v F is the Fermi velocity corresponding to the topo-logical bands, namely, (cid:126) v F = dε k /dk . By estimating theFermi velocities from Figs. 2(b) and (c) we find: Forthe w nanoribbon v F = 4 . × m s − , that ren-ders ξ = 0 . , while for w the Fermi velocity is . × m s − and ξ = 0 . . These estimates of ξ y ( nm ) . . . . L DO S (˚ A − e V − ) × − Averaged over x,z-directions exponential fit
FIG. 3. Local density of states (LDOS) at the Fermi level( E F = 0 ) for pristine w nanoribbon. The LDOS is averagedover the x and z directions and the y-axis corresponds to thetransverse direction along the ribbon width. The gray dashedline shows the exponential fit. for buckled bismuthene are slighlty larger than the valuesreported in experiments [7, 51] that obtain ξ ≈ . in SiC(0001) supported flat bismuthene. We stress thatthese are different material systems. The discrepancycan be explained by recalling that ∆ TG of the buckledbismuthene is smaller than the band gap of the flat one.A more quantitative estimate of ξ is taken from thelocal density of states (LDOS) averaged over the x, z di-rections as a function nanoribbon transversal axis y , seeFig. 3. Fitting an exponential function, we estimate thepenetration length as ξ = 0 . for the w ribbon, ingood agreement with ξ obtained using Eq. (5). Using ξ ,we can also estimate a threshold length for the LDOS de-cay for which the interedge states hybridization becomesnegligible. For instance, at distances ∆ y (cid:39) . fromthe edges, the LDOS decays by ∼ . Hence, we expectthat narrow nanoribbons with w (cid:46) y lack topologicalprotection due to the strong overlap of states localized atopposite system edges. The w nanoribbons are at thecrossover between non-protected and topologically pro-tected phases.Next, we investigate ribbons containing a single-vacancy. For a slab geometry, we calculate the formationenergy E V of these single vacancies systems using thefollowing expression [52] E V = E tot − ( E pristine + µ V N V ) (6)where E tot is the total energy of the single-vacancy sys-tem given by a fully relaxed DFT calculation, E pristine is the energy for the pristine ribbon, N V is the numberof vacancies (in our case N V = 1 ), and µ V is the chem-ical potential to remove a bismuth atom. Here, we take µ V as the energy per atom of the pristine bismuthenemonolayer. For purposes of comparison, we simulate a × supercell containing a vacancy. By doing so, weavoid that theses vacancies interact with their periodicimages obtaining a formation energy of .
04 eV , which isin agreement with previous reports [17].Figure 4 shows the conductance of the w and w ribbons in the presence of a single vacancy placed close toone of the system edges. For this calculation, the scatter- FIG. 4. Effect of a single vacancy at the edge in bismuthene:(a) structure of a single vacancy; conductance for (b) w nanoribbon and (c) w nanoribbon. The gray region indi-cates the topological gap. The black dashed line correspondsto the pristine case. ing region corresponds to a single building block N = 1 containing the vacancy. For the w ribbon the quantizedconductance is destroyed, showing that the edge statesare not robust against disorder. In turn, the w ribbonshows no deviation from perfect conductance within thetopological gap.Several studies on a variety of 2D materials indicatethat single-vacancies give rise to localized states [22, 53].It has been further shown that such disorder-induced lo-calized states cause the formation of local magnetic mo-ments [53–57] that, if close to the system edges, canbe detrimental for the topological protection, differentlyfrom dual topological insulators [58]. Let us study therelevance of these findings to bismuthene.Figure 5 shows the local density of states (LDOS) forthe w , w and w nanoribbons calculated for theenergy window . < E < . . The strong en-hancement of the LDOS centered around the vacancycorresponds to exponentially decaying orbitals paired af-ter the atom relaxation. Figure 5(a) shows that thevacancy states increase the overlap between inter-edgestates, whereas for the w nanoribbon the overlap issmall even in the presence of a vacancy. For the w ribbon shown in Fig. 5(c) the overlap between the edgestates and the vacancy localized states is negligible. FIG. 5. Top ( x - y axis) and laterial ( y - z axis) projectionsof the local density of states (LDOS) for single vacancy in(a) w , (b) w and (c) w nanoribbons. The isosurfacevalue for the LDOS is of . Å − eV − , corresponding to anenergy range . < E < . . Our ab initio fully relativistic calculations show nomagnetic moments generated by the defects. Nonethe-less, when SOC is turned off, spin polarized calculationsshow a magnetic moment of . µ B around the vacancy.This result rules out vacancy-induced local moments [22]as a method to hinder the topological protection.We now study the conductance as a function of va-cancy concentration. For that purpose, we put together N building blocks (as described in Sec. II) of length (cid:96) = 17 . Å, with a given concentration n V of randomlyplaced vacancies. Figure 6 shows the conductance fora single disorder realization as a function of energy fordifferent vacancy concentrations and nanoribbon widths.In all cases, the scattering region has a total length of ≈
157 nm .Figure 6 shows that vacancies, even in small concen-trations, have a strong effect on the transport propertiesof bismuthene nanoribbons. In all cases, the pristine con-ductance quantization is destroyed outside the topologi-cal gap ∆ TG . For instance, for n V ≈ . , one observesthe onset of localization ( G/G (cid:28) ) for w and a verystrong suppression of G for w , when compared to thepristine case. For energies outside ∆ TG , bismuthene is anordinary material and the results can be interpreted interms of standard disorder-induced backscattering mech-anisms. In what follows we discuss the behavior of G forenergies within the topological phase.Figures 6(a) to 6(d) correspond to the narrow widthnanoribbon w with vacancy concentration n V ranging from .
22 % to .
72 % . As discussed above, these sys-tems are not topological insulators. Interestingly, for in-creasing n V the w ribbons behave as trivial insulatorsfor energies outside the topological gap, but the conduc-tance is less suppressed inside ∆ TG . A similar behaviorhas been studied in graphene [59, 60] and MoS nanorib-bons [61, 62]. It is explained by the large momentumtransfer necessary to enable backscattering processes, asit can be inferred by inspecting the pristine electronicband structure. See, for instance [59, 60] for detaileddiscussions.Figures 6(e) to 6(h) correspond to w nanoribbonsof intermediate width and n V ranging from .
06 % to .
36 % . For energies outside ∆ TG the conductance isstrongly suppressed with respect to the pristine case, asexpected for a trivial disordered system. In contrast tothe w case, for n V (cid:46) .
36 % the system does not dis-play any Anderson localization features. More interest-ingly, the conductance is close to G within the topo-logical bandgap. The edge states are not fully protectedby topology: there are several energy intervals where G is strongly suppressed. The size of such intervals andthe G suppression grow with increasing n V . Interedgebackscattering processes are induced by disorder due tothe hybridization of the edge states with the randomlydistributed vacancy-induced localized states in the rib-bon. For weak disorder, G/G < only in narrow energywindows that depend on the disorder configuration. Thismechanism has already been shown to destroy topologi-cal protection for such narrow nanoribbons [24] and forwide ones in a tight-binding toy model with much largervacancy concentrations, n V (cid:39) [23].Let us now address the case of wide nanoribbons,namely, w (cid:29) ξ . Figures 6(i) to 6(l) correspond to w nanoribbons with n V ranging from .
13 % to .
48 % . Asin the previous cases, G decreases with increasing n V for energies outside ∆ TG . In distinction, the edge statesare topologically protected by time-reversal symmetryand G = G over most of the bulk band gap energies.Although topological protection is more robust than inthe previous cases, when n V (cid:38) .
30 % , G/G is stronglysuppressed for certain energy intervals within ∆ TG , seeFigs. 6(k) and (l), indicating the presence of vacancy-induced inter-edge backscattering processes.These results suggest that the number of vacancy-induced states necessary for an effective inter-edge hy-bridization the system edge states increases with thenanoribbon width w . IV. CONCLUSIONS AND OUTLOOK
In this paper we have studied the robustness of the con-ductance quantization against vacancy disorder in largescale nanoscopic buckled bismuthene nanoribbons at theQSH phase.We have found that vacancies give rise to mid-gaplocalized states and are non-magnetic, ruling out local G ( e / h ) w , 0.22 %(a) w , 0.06 %(e) w , 0.13 %(i) G ( e / h ) w , 0.36 %(b) w , 0.12 %(f) w , 0.19 %(j) G ( e / h ) w , 0.51 %(c) w , 0.24 %(g) w , 0.35 %(k) − . − . − . . . . . Energy (eV) G ( e / h ) w , 0.72 %(d) − . − . − . . . . . Energy (eV) w , 0.36 %(h) − . − . − . . . . . Energy (eV) w , 0.48 %(l) FIG. 6. Conductance G (in units of e /h ) as a function of the energy E (in eV) for several vacancy concentrations for the w , w , and w nanoribbons. The ribbon widths and vacancy concentrations are indicated in each panel. The blue linecorresponds to G in the presence of disorder, while the black dashed one stands for G in the pristine case. The gray regionindicates the topological gap. In all cases the length is
157 nm . magnetic moments [22] as a mechanism to destroy topo-logical protection. We have shown that the vacancy-induced mid-gap states ca give rise to inter-edge scat-tering processes. These depend on the edge state pen-etration depth, vacancy concentration, and nanoribbonwidth. The interplay of these quantities has been quali-tatively discussed in QSH tight-binding models [23] andwithin DFT for ultra-narrow systems [24]. Here, we haveestablished the presence of vacancy-induced inter-edgebackscattering processes in bismuthene nanoribbons ofrealistic sizes using ab initio techniques.Our calculations show different transport behavior forbulk and edge states, the first demonstrating localizationeffects and the latter showing robust topological responsefor low vacancy concentrations. At moderate n V values,topological protection is destroyed even for wide ribbons w (cid:29) ξ .Our findings are also applicable to other materials inthe QSH regime. Since the penetration depth is mate-rial dependent, we conclude that it is possible to engineerdifferent samples with QSH electronic transport behaviorin the presence of disorder. Our findings suggest an in-teresting application for a spintronic device, the level ofdoping or the width of the device modules the edge states degeneracy, therefore providing the ON/OFF states fora transistor switch. This might be obtained by changingthe chemical potentials on different leads or by apply-ing a gate voltage. In particular, for the topological w and w nanoribbons, characterized by a pronounceddrop in the conductance becoming broader for energiesaround the valence band as the level of vacancy concen-tration increases, allowing the possibility of reaching highON/OFF ratios. The mechanism to modify the conduc-tion in these nanoribbons does not rely on changing thetopology of its band structure as it would happen by ap-plying an electric or magnetic field. ACKNOWLEDGMENTS
This work is partially supported by the Coordenaçãode Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) - Finance Code 001; FAPESP (Grants19/04527-0, 16/14011-2, 17/18139-6, and 17/02317-2);CNPq (Grant 308801/2015-6), INCT-Nanocarbono, andFAPERJ (Grant E-26/202.882/2018). The authors ac-knowledge the Brazilian National Scientific ComputingLaboratory (LNCC) for computational resources of theSDumont computer. [1] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010). [2] J. E. Moore, Nature , 194 (2010). [3] Y. Ando, J. Phys. Soc. Jpn. , 102001 (2013).[4] A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. , 1(2016).[5] F. Giustino, M. Bibes, J. H. Lee, F. Trier, R. Valentí,S. M. Winter, Y.-W. Son, L. Taillefer, C. Heil, A. I.Figueroa, B. Plaçais, Q. Wu, O. V. Yazyev, E. P. A. M.Bakkers, J. Nygård, P. Forn-Díaz, S. de Franceschi,L. E. F. Foa Torres, J. McIver, A. Kumar, T. Low,R. Galceran, S. O. Valenzuela, M. V. Costache, A. Man-chon, E.-A. Kim, G. R. Schleder, A. Fazzio, and S. Roche,Journal of Physics: Materials 10.1088/2515-7639/abb74e(2020).[6] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. , 323 (2004).[7] F. Reis, G. Li, L. Dudy, M. Bauernfeind, S. Glass,W. Hanke, R. Thomale, J. Schäfer, and R. Claessen, Sci-ence , 287 (2017).[8] A. Marrazzo, M. Gibertini, D. Campi, N. Mounet, andN. Marzari, Nano Lett. , 8431 (2019).[9] C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801(2005).[10] C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 146802(2005).[11] I. K. Drozdov, A. Alexandradinata, S. Jeon, S. Nadj-Perge, H. Ji, R. J. Cava, B. Andrei Bernevig, and A. Yaz-dani, Nat. Phys. , 664 (2014).[12] B. Zhou, H.-Z. Lu, R.-L. Chu, S.-Q. Shen, and Q. Niu,Phys. Rev. Lett. , 246807 (2008).[13] M. Wada, S. Murakami, F. Freimuth, and G. Bihlmayer,Phys. Rev. B , 121310 (2011).[14] S. Shen, W. Shan, and H. Lu, SPIN , 33 (2011).[15] W. Zhang, R. Yu, H.-J. Zhang, X. Dai, and Z. Fang, NewJournal of Physics , 065013 (2010).[16] E. Aktürk, O. Üzengi Aktürk, and S. Ciraci, Phys. Rev.B , 10.1103/physrevb.94.014115 (2016).[17] Y. Kadioglu, S. B. Kilic, S. Demirci, O. U. Aktürk, E. Ak-türk, and S. Ciraci, Phys. Rev. B , 245424 (2017).[18] X. Wang, G. Bian, C. Xu, P. Wang, H. Hu, W. Zhou,S. A. Brown, and T.-C. Chiang, Nanotechnology ,395706 (2017).[19] M. Costa, G. R. Schleder, M. Buongiorno Nardelli,C. Lewenkopf, and A. Fazzio, Nano Lett. , 8941 (2019).[20] X. Ni, H. Huang, and F. Liu, Phys. Rev. B , 125114(2020).[21] X. Wang, C. Xu, H. Hu, P. Wang, G. Bian, W. Tan, S. A.Brown, and T.-C. Chiang, EPL , 27002 (2017).[22] P. Novelli, F. Taddei, A. K. Geim, and M. Polini, Phys.Rev. Lett. , 016601 (2019).[23] S. Tiwari, M. L. V. de Put, B. Sorée, and W. G. Vanden-berghe, 2D Mater. , 025011 (2019).[24] L. Vannucci, T. Olsen, and K. S. Thygesen, Phys. Rev.B , 155404 (2020).[25] P. Hohenberg and W. Kohn, Phys. Rev. , B864(1964).[26] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[27] G. R. Schleder, A. C. M. Padilha, C. M. Acosta,M. Costa, and A. Fazzio, Journal of Physics: Materials , 032001 (2019).[28] C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James, J. Phys. C: Solid State Phys. , 916 (1971).[29] Y. Meir and N. S. Wingreen, Phys. Rev. Lett. , 2512(1992).[30] S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press, 1995). [31] F. D. Novaes, A. A. J. R. d. Silva, and A. Fazzio, Braz.J. Phys. , 799 (2006).[32] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson,M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev.B , 6671 (1992).[33] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[34] G. Kresse and J. Furthmüller, Comput. Mater. Sci. , 15(1996).[35] G. Kresse and J. Furthmüller, Phys. Rev. B , 11169(1996).[36] G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999).[37] J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Jun-quera, P. Ordejón, and D. Sánchez-Portal, J. Phys.: Con-dens. Matter , 2745 (2002).[38] L. Fernández-Seivane, M. A. Oliveira, S. Sanvito, andJ. Ferrer, J. Phys.: Condens. Matter , 489001 (2007).[39] N. Troullier and J. L. Martins, Phys. Rev. B , 1993(1991).[40] S. Sanvito, C. J. Lambert, J. H. Jefferson, and A. M.Bratkovsky, Phys. Rev. B , 11936 (1999).[41] A. R. Rocha, M. Rossi, A. Fazzio, and A. J. R. da Silva,Phys. Rev. Lett. , 176803 (2008).[42] A. R. Rocha, M. Rossi, A. J. R. da Silva, and A. Fazzio,J. Phys. D: Appl. Phys. , 374002 (2010).[43] J. M. de Almeida, A. R. Rocha, A. J. R. da Silva, andA. Fazzio, Phys. Rev. B , 085412 (2011).[44] C. H. Lewenkopf and E. R. Mucciolo, J. Comput. Elec-tron. , 203 (2013).[45] A. Kormányos, G. Burkard, M. Gmitra, J. Fabian, V. Zó-lyomi, N. D. Drummond, and V. Fal’ko, 2D Mater. ,022001 (2015).[46] Y. Xue, S. Datta, and M. A. Ratner, Chem. Phys. ,151 (2002).[47] M. P. Lopez Sancho, J. M. Lopez Sancho, and J. Rubio,J. Phys. F: Met. Phys. , 851 (1985).[48] S. Shen, Topological Insulators: Dirac Equation in Con-densed Matter , Springer Series in Solid-State Sciences(Springer Singapore, 2017).[49] G. Dolcetto, M. Sassetti, and T. L. Schmidt, Riv. NuovoCimento Soc. Ital. Fis. , 113–154 (2016).[50] L. Chen, G. Cui, P. Zhang, X. Wang, H. Liu, andD. Wang, Phys. Chem. Chem. Phys. , 17206 (2014).[51] R. Stühler, F. Reis, T. Müller, T. Helbig, T. Schwemmer,R. Thomale, J. Schäfer, and R. Claessen, Nat. Phys. ,47 (2020).[52] S. Haldar, R. G. Amorim, B. Sanyal, R. H. Scheicher,and A. R. Rocha, RSC Adv. , 6702 (2016).[53] O. V. Yazyev, Rep. Prog. Phys. , 056501 (2010).[54] J. J. Palacios, J. Fernández-Rossier, and L. Brey, Phys.Rev. B , 195428 (2008).[55] M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gómez-Rodríguez, Phys. Rev. Lett. , 096804 (2010).[56] V. G. Miranda, L. G. G. V. Dias da Silva, and C. H.Lewenkopf, Phys. Rev. B , 075114 (2016).[57] V. Miranda, E. Mucciolo, and C. Lewenkopf, J. Phys.Chem. Solids , 169 (2019), spin-Orbit Coupled Ma-terials.[58] B. Focassio, G. R. Schleder, A. Pezo, M. Costa, andA. Fazzio, Phys. Rev. B , 045414 (2020).[59] K. Wakabayashi, Y. Takane, and M. Sigrist, Phys. Rev.Lett. , 036601 (2007).[60] L. R. F. Lima, F. A. Pinheiro, R. B. Capaz, C. H.Lewenkopf, and E. R. Mucciolo, Phys. Rev. B , 205111 (2012).[61] E. Ridolfi, L. R. F. Lima, E. R. Mucciolo, and C. H.Lewenkopf, Phys. Rev. B , 035430 (2017). [62] A. Pezo, M. P. Lima, M. Costa, and A. Fazzio, Phys.Chem. Chem. Phys.21