# Impact of complex adatom–induced interactions on quantum spin Hall phases

###### Abstract

Adsorbate engineering offers a seemingly simple approach to tailor spin–orbit interactions in atomically thin materials and thus to unlock the much sought-after topological insulating phases in two dimensions. However, the observation of an Anderson topological transition induced by heavy adatoms has proved extremely challenging despite substantial experimental efforts. Here, we present a multi-scale approach combining advanced first-principles methods and accurate single-electron descriptions of adatom–host interactions using graphene as a prototypical system. Our study reveals a surprisingly complex structure in the interactions mediated by random adatoms, including hitherto neglected hopping processes leading to strong valley mixing. We argue that the unexpected intervalley scattering strongly impacts the ground state at low adatom coverage, which would provide a compelling explanation for the absence of a topological gap in recent experimental reports on graphene. Our conjecture is confirmed by real-space Chern number calculations and large-scale quantum transport simulations in disordered samples. This resolves an important controversy and suggests that a detectable topological gap can be achieved by increasing the spatial range of the induced spin–orbit interactions on graphene, e.g., using nanoparticles.

The attachment of adsorbates to two-dimensional materials has attracted much interest in recent years, as a route to tailoring material properties and realising novel phenomena Adorb_G_abinitio_Peeters_08 (); PhotoLum_TMD_engineering (); Adsorb_G_Review_Lee_12 (); Adsorb_Various_TMD_Voiry_15 (); Adsorb_2D_Review_Ryder_16 (). In graphene, adatoms have been shown to induce band gaps BGap1_Fluorographene_Zhu_10 (); BGap2_Fluorographene_Nair_10 (); BGap3_Fluorographene_Snow_10 (); BGap4_Fluoride_DFT_Peeters_10 (); BGap5_Aryl_Lau_11 (), magnetic moments G_MAGN_Ochoa_08 (); G_MAGN_Nair_12 (); G_MAGN_Herrero_16 () and even superconductivity Supercond_G_Profefa_12 (); Supercond_G_Li_13 (); Supercond_G_Ludbrook (); Supercond_G_Chapman_16 ().

Adsorbate engineering could likewise provide atomic control over fundamental spin–orbit phenomena, such as spin relaxation SpinRelax_G_Hernando_09 (); SpinRelax_G_Pi_10 (); SpinRelax_G_Federov_13 (); SpinRelax_G_Tuan_14 () and Mott (skew) scattering SHE_G_Pachoud_14 (); SHE_Graphene_Ferreira_14 (); SHE_Graphene_Ozyilmaz_Ferreira_14 (); SHE_Huang_16 (). Recent studies have predicted that the dilute assembly of heavy adatoms can massively enhance the weak spin–orbit energy gap of graphene QAHE_G_Adatom_Zhang (); QSH_G_Adatom_Weeks_11 (); QSH_G_Adatom_Jiang_12 (); QSH_G_Adatom_Brey (), opening a promising route towards the realisation of nontrivial topological insulating phases, including the quantum spin Hall (QSH) state QSH_G_Kane_Mele_05 (). However, transport measurements on samples decorated with heavy species, including In and Ir, have yet to show any signature of topological gap opening QSH_Exp_Adatom_1 (); QSH_Exp_Adatom_2 (); QSH_Exp_Adatom_3 (); QSH_Exp_Adatom_4 (); QSH_Exp_Adatom_5 (). In this work, we show that a thorough treatment of disorder—combining accurate model Hamiltonians with quantum transport simulations—is essential to predict the topological character of adatom-engineered systems and reconcile this contradiction. Our approach reveals that randomly distributed heavy adatoms on graphene give rise to scattering between inequivalent valleys in the band structure, hindering the emergence of topologically protected edge states even in the absence of extrinsic factors, such as adatom clustering QSH_clustering_2014 (). This resolves a controversy regarding the nature of spin–orbit interactions in adatom-decorated graphene and indicates that decoration with small clusters or nanoparticles, for which intervalley scattering is strongly reduced, may offer a route towards the realisation of the much sought-after QSH phase.

Pristine graphene is a QSH insulator, but the smallness of its intrinsic spin–orbit interaction (in the range of 25–50 eV Intrinsic_SOC_G_Huertas-Hernando_06 (); Intrinsic_SOC_G_Min_06 (); Intrinsic_SOC_Konschuh_10 (); Intrinsic_SOC_Sichau_17 ()) has thusfar precluded achieving the dissipationless quantum transport regime QSH_QWells_Konig_07 (); QSH_QWells_Inverted_Liu_08 (); QSH_QWells_Inverted_Du_15 (). The opening of a detectable QSH gap requires a massive enhancement of graphene’s characteristic spin–orbit coupling (SOC), which preserves spin angular momentum QSH_G_Kane_Mele_05 (). Previous work suggested that this can be achieved via decoration with nonmagnetic adatoms with a -outer electron shell QSH_G_Adatom_Weeks_11 (), as they induce spin–conserving ’intrinsic-like’ SOC SHE_G_Pachoud_14 (). As a prototypical heavy element, we consider thallium (Tl) QSH_G_Adatom_Weeks_11 (). The quasiparticle band structure of the thallium adatom on graphene was obtained employing a fully-relativistic ab initio GW approach; see Supplemental Material (SM) SM () for details.

First-principles multi-scale approach.—We use our first-principles supercell calculations to parameterize a single-electron Hamiltonian capturing all the relevant interactions mediated by (dilute) adatoms embedded in large area graphene samples. The first step is to derive a graphene–single-adatom tight-binding (TB) model that faithfully reproduces the ab initio band structure. Enforcing time-reversal symmetry and invariance with respect to the point group, one easily finds SHE_G_Pachoud_14 (); QSH_G_Adatom_Weeks_11 (), where

(1) | ||||

(2) | ||||

(3) |

The first two terms are the Hamiltonians of -electrons on graphene and 6-states of a Tl atom, respectively. and are the corresponding fermionic creation operators, are Pauli matrices acting on the spin space and . The sites adjacent to the adatom define a hexagonal plaquette, . is the adatom–graphene hybridization term written as a function of the plaquette operator for states with definite angular momentum, Note_Minimal_Model (). Next-nearest and third-nearest neighbor corrections () are included in order to improve agreement to the first-principles results.

The minimal Hamiltonian [Eqs. (1)-(3)] contains 9 parameters: the C-C hoppings, , and , the local chemical potential change on C sites next to Tl , the Tl outer-shell energies, and , the Tl spin-orbit energy and the C-Tl hoppings, and . We adjust these parameters until the band structures quantitatively reproduce the first-principles calculations over a window of eV around the Dirac point (see SM SM ()). The initial guess for the parameters is informed by a direct evaluation of hopping integrals between atom-centred maximally localized Wannier functions SM (). The quasiparticle band structure obtained from density functional theory (DFT)- and GW-parameterized minimal TB models is shown in Fig. 1. Bands below the Dirac point () derive mostly from graphene states. The flat band with energy eV is a Tl 6 state. Interestingly, the GW corrections are seen to bring this band closer to the Dirac point. Moreover, the gap at the Dirac point is 23 meV significantly larger than the DFT value of 13 meV. To what extent the optimistic first-principles estimates signal a measurable topological gap in real samples will depend on a delicate competition between spin-conserving SOC and two other interactions mediated by disorder, which we unveil in what follows.

Adatom scattering potential.—In realistic conditions, dilute adatoms occupy random positions and thus act as scattering centres. The information on the adatom scattering potential is contained in the local density of states (LDOS) LDOS_graphene_Bena_08 (). The crucial step in our multi-scale approach consists of deriving a (graphene-only) TB model for the scattering potential compatible with the first-principles results for the adatom–graphene supercell (see Fig. 1). To this end, we trace out the adatom degrees of freedom in Eqs. (1)-(3) through a Löwdin transformation. Formally, the resulting graphene–only TB Hamiltonian is given by , where is the real-space self-energy generated by a single adatom SM (). This term breaks translational symmetry and thus acts as a bona fide disorder interaction. Finally, the Hamiltonian of graphene with a dilute adatom coverage is obtained by adding up independent contributions from adatoms located at random plaquettes, (. This procedure has two advantages. It captures all short-range interactions induced by the adatom (see below). Second, being based on a graphene-only description, it allows a straightforward interpretation of quantum transport calculations. Explicit evaluation of the self-energy gives rise to the following effective interaction Hamiltonian, , with

(4) |

where equals for circulation around the -th plaquette (anti-)clockwise and with . The effective hoppings are defined by and

(5) |

where , , and . The first terms ( and ) modify the hopping between nearest-neighbor (NN) atoms, while the second line describes next-nearest-neighbor (NNN) hoppings, including a chiral component (), which—in the absence of the other terms—transforms graphene into a QSH insulator QSH_G_Adatom_Weeks_11 (). The terms in the third line capture hoppings between C atoms on opposite sides of the adatom () and spin-flip processes between all pairs of sites in the impurity plaquette (). The latter is a Rashba-type interaction, which is vanishing small near the Dirac points and thus can be safely neglected QSH_G_Adatom_Weeks_11 (); SHE_Graphene_Ferreira_14 (). The relevant interactions are visualized in Fig. 2.

Remarkably, the effective Hamiltonian obtained here by a rigorous adatom-decimation procedure is far more complex than previous models QSH_G_Adatom_Weeks_11 (). Importantly, and for heavy species can be significantly larger than the chiral NNN hopping. From the ab initio parameters derived for Tl SM () we obtain at . To shed further light on the significance of hitherto neglected terms [(c)-(d) in Fig. 2], we derive a long-wavelength effective description. As customary, we introduce the field operators, , with describing the projection low-energy states on the () sublattice [at ] for spin . Substituting in Eq. (4) and performing a series expansion around the inequivalent Dirac points, , one obtains the effective interaction: Comment (), with

(6) |

to leading order in , and where we omitted a scalar term (see SM SM ()). and denote Pauli matrices acting on sublattice and valley degrees of freedom, respectively, and describes the spatial profile of the adatom potential. The first term derives from the NNN hopping, QSH_G_Kane_Mele_05 (). The remaining terms are scalar and spin–orbit interactions connecting valleys, with strengths and , respectively. For -outer-shell adatoms [see Eq. (5)] and thus , i.e., intra- and inter-valley spin–orbit scattering processes must be considered on equal footing. Based on general arguments for disordered zero-gap semiconductors, one expects that the mixing of states at inequivalent degeneracy points is detrimental for the topological phase Hasan_Kane_10 (); Fradkin_86 (); Suzuura_Ando_02 (). An estimation using DFT-optimized parameters gives eV and eV at . Such a dominance of intervalley processes in the coarse-grained description is a strong indication that the topological gap displayed by Tl–graphene supercells (Fig. 1) will be fragile in a disordered scenario, which would naturally explain the negative experimental results QSH_Exp_Adatom_1 (); QSH_Exp_Adatom_2 (); QSH_Exp_Adatom_3 (); QSH_Exp_Adatom_4 (); QSH_Exp_Adatom_5 (). This idea is reinforced by the fact that electrons in graphene are very sensitive to valley mixing processes with an origin in short-range impurities, as the respective Friedel oscillations are known to decay as at large distances (as opposed to the law from intravalley scattering LDOS_graphene_Bena_08 ()).

Real-space quantum transport study.—To investigate the implications of the complex structure of the effective adatom potential, we have carried out transport calculations within the Landauer-Büttiker framework. In the QSH regime, a pair of counter-propagating gapless edge states protected from elastic backscattering emerge at the interfaces to vacuum [Fig. 3 (a)] QSH_G_Kane_Mele_05 (). To probe the robustness of the extrinsic QSH insulating phase and its concomitant helical edge states, we calculated the two-terminal conductance of large armchair nanoribbons (width nm, and length nm) with randomly distributed adatoms connected to pristine graphene leads. The central channel contains in excess of 3.5 million atoms and efficient recursion techniques are employed to solve a system of this large size SM (). The smoking gun for the topologically-protected edge states is the emergence of quantized conductance with a plateau width proportional to the SOC strength supp15_Buttiker_98 (); Roth_09 (). To probe the experimentally relevant adatom coverages would require prohibitively large computational domains, in order to resolve the typically small spin–orbit gaps (eV), where is the adatom coverage. To overcome this difficulty, we rescaled uniformly the effective hoppings, with . The main findings are summarized in Fig. 3. When only the intrinsic-type SOC term is included (, the two-terminal conductance exhibits a plateau at small energies with [Fig. 3 (b) (black dots)]. The variance is found to be zero up to numerical accuracy. Such a perfectly quantized effect shows that the nanoribbon has been transformed into a QSH insulator. The plateau width precisely saturates the upper bound , which is the topological gap obtained by Weeks et al QSH_G_Adatom_Weeks_11 (). However, the plateau shrinks when is turned on (blue dots) indicating a closing of the topological gap. This effect is accompanied by significative fluctuations, , showing that states delocalized through the nanoribbon contribute now to the electronic transport. This striking decay of helical edge states into the bulk is confirmed by the numerical evaluation of the spin-polarized bond current density maps; see SM SM (). The remaining adatom interactions affect the QSH phase in two different ways: (i) the imaginary NN hopping shifts the conductance plateau (not shown) and (ii) the hopping enlarges the central plateau and increases the fluctuations.

We now provide compelling evidence that QSH phases induced by random dilute adatoms are especially fragile in truly disordered scenarios, where additional sources of scattering are unavoidable. To this end, we introduce a topological defect (vacancy) in the nanoribbon, by cutting all bonds adjacent to one carbon atom. Vacancies introduce (quasi)localized states at zero energy strongly impacting the graphene electronic properties ZEM_Pereira_06 (); ZEM_Ugeda_10 (); ZEM_Hafner_14 (); ZEM_Hafner_15 (); ZEM_Weik_16 (); ZEM_nanoribbon_15 (). Given our choice of a metallic armchair nanoribbon, we locate the vacancy on a site with finite density of states GNR_defect_middle (). This choice guarantees that the vacancy works as a resonant scatter (introducing mid-gap states), leading to a strong suppression of the conductance at low energy, as [Fig. 3(c)]. Adding a small coverage of idealized adatoms with only NNN hopping () gives rise to the expected quantization of the conductance, as the helical edge states resulting from the SOC enhancement can perfectly avoid the vacancy (no backscattering). Quite strikingly, when a NN hopping correction (typical of Tl adatoms in rows 5 and 6 of the periodic table) is turned on, the conductance acquires its basic shape prior to adatom decoration, unambiguously demonstrating the inherent fragility of the QSH phase due to the activation of intervalley processes. The dependence with the adatom coverage is shown in panel Fig. 3(d), which also shows the conductance at low energies is further degraded when and are taken simultaneously.

Chern number in real-space.—At this stage, we have firmly established that edge states in adatom-decorated graphene are intrinsically unstable due to hitherto neglected hoppings, and (Fig. 2). As shown by our multi-scale theory bridging advanced first-principles calculations and accurate TB models, such hoppings are typically one order of magnitude larger than the chiral term inducing the topological phase (). To investigate the onset of the topological phase transition in more detail, we evaluate the spin Chern number of bulk states, defined by , where is the first Chern integer QSH_Sheng_06 (). To compute for disordered configurations of adatoms, we employ an efficient gauge-invariant approach developed in Refs. Chern_Fukui_05 (); Chern_Zhang_13 (). This allow us to assess the topological order for realistic TB parameters (i.e., no rescaling is required). In Fig. 3(e), we show the dependence of the average Chern number on the adatom coverage. The robustness of with respect to and its quick suppression at low concentrations when and are turned on is in perfect agreement with the previous conclusions based on quantum transport simulations in nanoribbon geometry. When all adatom-induced interactions are considered on equal footing, the onset to the transition to a topologically trivial phase is found to occur around 5% for typical values of the parameters. When approaching 1-2% coverage, the fluctuations are increasingly larger predicting the closing of the topological gap in the relevant experimental regime. The strong dependence of on the adatom coverage is reminiscent of Anderson topological insulators Li_TAInsulators_09 (), for which the character of the insulating ground state is known to critically depend on the disorder strength.

Our findings have several important consequences. The absence of topological gap signatures on recent measurements in adatom-decorated graphene QSH_Exp_Adatom_1 (); QSH_Exp_Adatom_2 (); QSH_Exp_Adatom_3 (); QSH_Exp_Adatom_4 (); QSH_Exp_Adatom_5 () is naturally explained by valley mixing processes beyond simple model Hamiltonians. The adatom-induced intervalley scattering uncovered in this work can be mitigated by increasing the spatial range of the interactions mediated by the adsorbate, thereby providing a possible path towards the engineering of a quantum spin Hall insulator in graphene, i.e., its decoration with dilute heavy nano-particles. These results highlight the importance of seamless multi-scale approaches bridging first-principles parameterized model Hamiltonians and large-scale quantum transport calculations for the predictive modelling of adatom–host systems.

Data availability statement.—The data underlying this paper are available from the Figshare database figshare ().

Acknowledgments.—We acknowledge support the Engineering and Physical Sciences Research Council (EPRSC) under Grants No. EP/N005244/1 (J.L.), and EP/K003151/1 and EP/P023843/1 (K.M.), the Thomas Young Centre under Grant No. TYC-101 (J.L.), the Brazilian funding agency CAPES under Project No. 13703/13-7 (F.J.S.), and the European Research Council (ERC) under the ERC-consolidator Grant No. 681405-DYNASORE (F.J.S.). Via J.L. and K.M.’s membership of the UK’s HEC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work used the ARCHER UK National Supercomputing Service. DAB acknowledges the support from Mackpesquisa and FAPESP under grant 2012/50259-8. E.V.C. acknowledges partial support from FCT-Portugal through Grant No. UID/CTM/04540/2013. A.F. gratefully acknowledges the support from the Royal Society through a Royal Society University Research Fellowship and partial funding from EPSRC (Grant No. EP/N004817/1). F.J.S thanks Filipe S. M. Guimarães for the support and helpful discussions. E.V.C. and A.F. thank M.P. López-Sancho and M.A.H. Vozmediano for useful discussions.

## References

- (1) O. Leenaerts, B. Partoens, and F. M. Peeters, Phys. Rev. B 77, 125416 (2008).
- (2) T. Kuila, S. Bose, A. K. Mishra, P. Khanra, N. H. Kim, J. H. Lee, Prog. Mater. Sci. 57, 1061 (2012).
- (3) H. Nan, Z. Wang, W. Wang, Z. Liang, Y. Lu, Q. Chen, D. He, P. Tan, F. Miao, X. Wang, J. Wang, and Z. Ni. ACS Nano 8, 5738 (2014).
- (4) D. Voiry, A. Goswami, R. Kappera, C. C. C. Silva, D. Kaplan, T. Fujita, M. Chen, T. Asefa, and M. Chhowalla. Nature Chemistry 7, 45 (2015).
- (5) C. R. Ryder, J. D. Wood, S. A. Wells, and M. C. Hersam, ACS Nano, 10, 3900 (2016).
- (6) S. H. Cheng, K. Zou, F. Okino, H. Gutierrez, A. Gupta, N. Shen, P. Eklund, J. Sofo, and J. Zhu, Phys. Rev. B 81, 205435 (2010).
- (7) R. R. Nair, W. C. Ren, R. Jalil, I. Riaz, V. G. Kravets, L. Britnell, P. Blake, F. Schedin, A. S. Mayorov, S. J. Yuan, M. I. Katsnelson, H. M. Cheng, W. Strupinski, L. G. Bulusheva, A. V. Okotrub, I. V. Grigorieva, A. N. Grigorenko, K. S. Novoselov, and A. K. Geim, Small 6, 2877 (2010).
- (8) J. T. Robinson, J. S. Burgess, C. E. Junkermeier, S. C. Badescu, T. L. Reinecke, F. K. Perkins, M. K. Zalalutdniov, J. W. Baldwin, J. C. Culbertson, P. E. Sheehan, and E. S. Snow, Nano Lett. 10, 3001 (2010).
- (9) O. Leenaerts, H. Peelaers, A. Hernandez-Nieves, B. Partoens, and F. Peeters, Phys. Rev. B 82, 195436 (2010).
- (10) H. Zhang, E. Bekyarova, J. W. Huang, Z. Zhao, W. Bao, F. Wang, R. C. Haddon, and C. N. Lau, Nano Lett. 1, 4047 (2011).
- (11) B. Uchoa, V. N. Kotov, N. M. R. Peres, and A. H. Castro Neto. Phys. Rev. Lett. 101, 026805 (2008).
- (12) R. R. Nair, M. Sepioni, I-Ling Tsai, O. Lehtinen, J. Keinonen, A. V. Krasheninnikov, T. Thomson, A. K. Geim and I. V. Grigorieva, Nat. Phys. 8, 199 (2012).
- (13) H. G.-Herrero, J. M. G.-Rodríguez, P. Mallet, M. Moaied, J. J. Palacios, C. Salgado, M. M. Ugeda, J.-Y. Veuillen, F. Yndurain, and I. Brihuega, Science 352, 437 (2016).
- (14) G. Profeta, M. Calandra, and F. Mauri. Nature Physics 8, 131 (2012).
- (15) K. Li, X. Feng, W. Zhang, Y. Ou, L. Chen, K. He, L.-L. Wang, L. Guo1, G. Liu, Q.-K. Xue, and X. Ma. Appl. Phys. Lett. 103, 062601 (2013).
- (16) B. M. Ludbrook, G. Levy, P. Nigge, M. Zonno, M. Schneider, D. J. Dvorak, C. N. Veenstra, S. Zhdanovich, D. Wong, P. Dosanjh, C. Straßer, A. Stöhr, S. Forti, C. R. Ast, U. Starke, and A. Damascelli. Proc. Natl Acad. Sci. 112, 11795 (2015).
- (17) J. Chapman, Y. Su, C. A. Howard, D. Kundys, A. N. Grigorenko, F. Guinea, A. K. Geim, I. V. Grigorieva, and R. R. Nair. Sci Rep. 6, 23254 (2016).
- (18) D. Huertas-Hernando, F. Guinea, and Arne Brataas. Phys. Rev. Lett. 103, 146801 (2009).
- (19) K. Pi, W. Han, K. M. McCreary, A. G. Swartz, Y. Li, and R. K. Kawakami. Phys. Rev. Lett. 104, 187201 (2010).
- (20) D. V. Fedorov, M. Gradhand, S. Ostanin, I. V. Maznichenko, A. Ernst, J. Fabian, and I. Mertig. Phys. Rev. Lett. 110, 156602 (2013).
- (21) D. V. Tuan, F. Ortmann, D. Soriano, S. O. Valenzuela, and S. Roche. Nature Physics 10, 857 (2014).
- (22) A. Pachoud, A. Ferreira, B. Ozyilmaz, and A. H. Castro Neto. Phys. Rev. B 90, 035444 (2014).
- (23) A. Ferreira, T.G. Rappoport, M.A. Cazalilla, and A.H. Castro Neto, Phys. Rev. Lett. 112, 066601 (2014).
- (24) C. Huang, Y. D. Chong, and M. A. Cazalilla. Phys. Rev. B 94, 085414 (2016).
- (25) J. Balakrishnan, G.K.W. Koon, A. Avsar, Y. Ho, J.H. Lee, M. Jaiswal, S.-J. Baeck, J.-H. Ahn, A. Ferreira, M.A. Cazalilla, A.H. Castro Neto, B. Ozyilmaz, Nat. Comm. 5, 4748 (2014).
- (26) C. Weeks, J. Hu, J. Alicea, M. Franz, and R. Wu, Phys. Rev. X 1, 021001 (2011).
- (27) H. Zhang, C. Lazo, S. Blügel, S. Heinze, and Y. Mokrousov. Phys. Rev. Lett. 108, 056802 (2012).
- (28) H. Jiang, Z. Qiao, H. Liu, J. Shi, and Q. Niu. Phys. Rev. Lett. 109, 116803 (2012).
- (29) L. Brey. Phys. Rev. B 92, 235444 (2015).
- (30) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005).
- (31) U. Chandni, E. A. Henriksen, and J. P. Eisenstein. Phys. Rev. B 91, 245402 (2015).
- (32) Y. Wang, S. Xiao, X. Cai, W. Bao, J. R.-Robey, and M. S. Fuhrer. Sci. Rep. 5, 15764 (2015).
- (33) Z. Jia, B. Yan, J. Niu, Q. Han, R. Zhu, D. Yu, and X. Wu. Phys. Rev. B 91, 085411 (2015).
- (34) Y. Wang, X. Cai, J. R.-Robey, and M. S. Fuhrer. Phys. Rev. B 92, 161411(R) (2015).
- (35) J. A. Elias and E. A. Henriksen. Phys. Rev. B 95, 075405 (2017).
- (36) A. Cresti, D. V. Tuan, D. Soriano, A. W. Cummings, S. Roche. Phys. Rev. Lett. 113, 246603 (2014).
- (37) M. Z. Hasan and C. L. Kane. Rev. Mod. Phys. 82, 3045 (2010).
- (38) D. N. Sheng, Z. Y. Weng, L. Sheng, and F. D. M. Haldane. Phys. Rev. Lett. 97, 036808 (2006).
- (39) D. Huertas-Hernando, F. Guinea, and A. Brataas, Phys. Rev. B 74, 155426 (2006).
- (40) H. Min, J. E. Hill, N. A. Sinitsyn, B. R. Sahu, L. Kleinman, and A. H. MacDonald. Phys. Rev. B 74, 165310 (2006).
- (41) S. Konschuh, M. Gmitra, and J. Fabian, Phys. Rev. B 82, 245412 (2010).
- (42) J. Sichau, M. Prada, T. J. Lyon, B. Bosnjak, L. Tiemann, R. H. Blick. pre-print arxiv:1709.05705 (2017).
- (43) M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, S.-C. Zhang. Science 318, 766 (2007).
- (44) C. Liu, T. L. Hughes, X.-L. Qi, K. Wang, and S.-C. Zhang. Phys. Rev. Lett. 100, 236601 (2008).
- (45) L. Du, I. Knez, G. Sullivan, and R.-R. Du. Phys. Rev. Lett. 114, 096802 (2015).
- (46) See Supplemental Material at attached for methodology, derivation of the effective Hamiltonian and additional quantum transport and topological invariant calculations, which includes Refs. supp1_QuantumEspresso (); supp2_hamann2013optimized (); supp3_perdew1996generalized (); supp4_monkhorst1976special (); supp5_CRC (); supp6_LouieHybertsen (); supp7_BGWpaper (); supp8_StaticRemainder (); supp9_W_Mostofi_08 (); supp10_W_Marzari_97 (); supp11_W_Souza_01 (); supp12_W_Jung_13 (); supp13_GNR_Son_06 (); supp14_GNR_Zheng07 (); supp15_Buttiker_98 (); supp16_Landauer_70 (); supp17_Caroli_71 (); supp18_Datta_Book (); supp19_Fisher_81 (); supp20_Lewenkopf_13 (); supp21_AntiPekka ().
- (47) P. Giannozzi, et al. Journal of Physics: Condensed Matter 21, 395502 (2009).
- (48) D. Hamann, Phys. Rev. B 88, 085117 (2013).
- (49) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (50) H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
- (51) R. C. Weast and M. J. Astle, CRC Handbook of Chemistry and Physics.
- (52) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
- (53) J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, and S. G. Louie, Comput. Phys. Commun. 183, 1269 (2012).
- (54) J. Deslippe, G. Samsonidze, M. Jain, M. L. Cohen, and S. G. Louie, Phys. Rev. B 87, 165124 (2013).
- (55) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
- (56) N. Marzari and D. Vanderbilt, Phys. Rev. B 56 12847 (1997).
- (57) I. Souza, N. Marzari and D. Vanderbilt, Phys. Rev. B 65 035109 (2001).
- (58) J. Jung and A. H. MacDonald, Phys. Rev. B 87, 195450 (2013).
- (59) Y.-W. Son, M. L. Cohen, and S. G. Louie. Phys. Rev. Lett. 97, 216803 (2006).
- (60) H. Zheng, Z. F. Wang, T. Luo, Q. W. Shi, and J. Chen. Phys. Rev. B 75, 165414 (2007).
- (61) M. Büttiker. Phys. Rev. B 38, 9375 (1988).
- (62) R. Landauer. Philosophical Magazine, 21 863 (1970).
- (63) C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James. J. Phys. C: Solid State Phys. 4, 916 (1971).
- (64) S. Datta. Electronic Transport in Mesoscopic Systems. Cambridge University Press, Cambridge (1995).
- (65) D. S. Fisher and P. A. Lee. Phys. Rev. B, 23 6851 (1981).
- (66) C. H. Lewenkopf and E. R. Mucciolo. J Comp. Elect. 12, 203 (2013).
- (67) H. Haug and A.-P. Jauho. Quantum Kinetics in Transport and Optics of Semiconductors. Springer; 2nd edition (2007).
- (68) Spin-dependent tunneling processes (e.g., ) have negligible probability and thus are omitted. Moreover, crystal-field effects are neglected QSH_G_Adatom_Weeks_11 (); SM (); Cysne_18 ().
- (69) T. P. Cysne, A. Ferreira, and T. G. Rappoport. Phys. Rev. B 98, 045407 (2018).
- (70) C. Bena. Phys. Rev. Lett. 100, 076601 (2008).
- (71) The low-energy theory for generic physisorbed adatoms was first derived in Ref. SHE_G_Pachoud_14 ().
- (72) E. Fradkin. Phys. Rev. B 33, 3257; ibidem 3263 (1986).
- (73) H. Suzuura and T. Ando. Phys. Rev. Lett. 89, 266603 (2002).
- (74) A. Roth, C. Brüne, H. Buhmann, L. W. Molenkamp, J. Maciejko, X.-L. Qi, S.-C. Zhang, Science 325, 294 (2009).
- (75) V. M. Pereira, F. Guinea, J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 96, 036801 (2006).
- (76) M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gomez-Rodriguez, Phys. Rev. Lett. 104, 096804 (2010).
- (77) V. Häfner, J. Schindler, N. Weik, T. Mayer, S. Balakrishnan, R. Narayanan, S. Bera, and F. Evers. Phys. Rev. Lett. 113, 186802 (2014).
- (78) A. Ferreira and E. Mucciolo, Phys. Rev. Lett. 115, 106601 (2015).
- (79) N. Weik, Phys. Rev. B 94, 064204 (2016).
- (80) H.-Y. Deng and K. Wakabayashi, Phys. Rev. B 91, 035425 (2015)
- (81) J.-Y. Yan, P. Zhang, B. Sun, H.-Z. Lu, Z. Wang, S. Duan, and X.G. Zhao. Phys. Rev. B 79, 115403 (2009).
- (82) T. Fukui, Y. Hatsugai, and H. Suzuki, J. Phys. Soc. Jpn. 74, 1674 (2005).
- (83) Y. F. Zhang, Y. Y. Yang, Y. Ju, L. Sheng, D. N. Sheng, R. Shen, and D. Y. Xing, Chin. Phys. B 22, 117312 (2013).
- (84) J. Li, Rui-L.Chu, J. K. Jain, and S.-Q. Shen. Phys. Rev. Lett. 102, 136806 (2009).
- (85) DOI: https://doi.org/10.6084/m9.figshare.6934883.v1

## Appendix A First-Principles Calculations

### a.1 Framework

We carry out ab initio density-functional theory (DFT) calculations for thallium (Tl) atoms adsorbed on graphene. We employ a graphene supercell containing 32 carbon (C) atoms and a single Tl adatom located above a hollow site, i.e. above the center of a graphene hexagon. Our calculations are carried out within the planewave-pseudopotential approach as implemented in the QUANTUM ESPRESSO supp1_QuantumEspresso () software package. For the electron-ion interaction, we employ fully-relativistic, multiple-projector, normconserving pseudopotentials, that were generated with the ONCVPSP code supp2_hamann2013optimized (). For the Tl pseudopotential, the and semicore states were included in the valence states. We further use the the PBE exchange-correlation energy functional supp3_perdew1996generalized (), a Monkhorst-Pack grid to sample the Brillouin zone of the supercell supp4_monkhorst1976special (), a 60 Ry plane-wave cutoff, a Gaussian smearing of 0.01 Ry and a separation of 20 Å between periodic images of the graphene sheet, which we assume to be stacked in the z-direction.

### a.2 Fully-relativistic structure relaxation

We first relax the structure to determine the atomic ground state configuration. Interestingly, we observe significant differences for the height of the Tl atom above the graphene sheet for scalar-relativistic and fully-relativistic DFT calculations: while the scalar-relativistic approach yields at Tl height of 2.58 Å above the graphene plane (defined as the average z-coordinate of all C atoms), the fully-relativistic approach predicts a significantly larger height of 2.87 Å. Importantly, this increase of the Tl heights results in a reduction of the band gap at the Dirac point from 24 meV to 13 meV.

method | IP |
---|---|

SR-KS | 8.19 |

SR-GW | 5.31 |

SR-SCF | 5.17 |

FR-KS | 9.28 |

FR-GW | 5.95 |

FR-SCF | 5.95 |

Experiment | 6.1 |

### a.3 GW corrections

To overcome the well-known band gap problem of DFT supp6_LouieHybertsen ()
and accurately determine the electronic structure of the adatom-graphene
system, we carry out a fully-relativistic GW calculation using the
BerkeleyGW program package supp7_BGWpaper (). As a test, we
first carry out calculations on an isolated Tl atom and an unperturbed
graphene sheet. We place the Tl atom in a cubic unit cell with a linear
extent of 20 bohr. The dielectric matrix is calculated using a 15 Ry
plane-wave cutoff and 700 bands and then extended to finite frequencies
using a generalized plasmon-pole model supp6_LouieHybertsen ().
Moreover, we employ a spherically truncated Coulomb interaction to
avoid spurious interactions between periodic images and a static remainder
correction to include contributions from high-energy unoccupied states
in the self-energy supp8_StaticRemainder (). To obtain the ionization
potential (IP) of the neutral atom, we consider a Tl ion with one
less electron than the neutral atom ^{1}^{1}1Carrying out GW calculations for a neutral Tl atom is conceptually
difficult because it has a partially filled 6p shell. The negatively-charged
ion, on the other hand, has a closed-shell configuration and calculations
are straightforward. and compute its electron affinity. Table 1 summarizes
our results and shows that the IP from fully-relativistic GW agrees
well with the experimental data, while the estimate obtained from
the Kohn-Sham orbital energy of fully-relativistic DFT deviates by
more than 3 eV from experiment. Note also that fully-relativistic
GW predicts a significantly smaller value of 0.95 eV for the spin-orbit
splitting of the states compared to fully-relativistic DFT,
which predicts 1.43 eV. We also computed the GW correction to the
DFT IP of unperturbed graphene and find it to be only 0.1 eV. This
analysis clearly shows that a GW calculation is necessary to reliably
describe the level alignment of the Tl-graphene system as the accuracy
of DFT-PBE for the IP of localized systems, such as atoms and molecules,
is significantly worse than for extended systems, such as graphene.

Finally, we carry out a fully-relativistic GW calculation for the combined Tl-graphene system. We use a plane-wave cutoff of 15 Ry for the dielectric matrix, a sum over 1000 states, a slab-truncated Coulomb interaction, a static remainder correction, a generalized plasmon-pole model and a kpoint grid. Fig. 4(a) shows the ab initio GW band structure (red dots) and compares it to the ab initio DFT result. In the graph, the bands below the Dirac point are graphene states. The flat bands above the Dirac point, however, derive from Tl p-states. We observe that the GW correction reduces the separation between the Dirac point and the lowest Tl -state by several 100 meV. Most importantly, we find that GW increases the band gap at the Dirac point from its DFT values of 13 meV to 19 meV.

## Appendix B Tight-Binding Models

### b.1 Wannierization

We employ the WANNIER90 code to construct maximally localized Wannier functions (MLWFs) for Tl-decorated graphene supp9_W_Mostofi_08 (); supp10_W_Marzari_97 (); supp11_W_Souza_01 (). The initial band-structure calculations are performed as described below in QUANTUM ESPRESSO supp1_QuantumEspresso (). This approach has been applied previously to pristine graphene supp12_W_Jung_13 (). Here we extend the approach to non-collinear calculations on Tl-decorated graphene. A fully self-consistent calculation is performed on the supercell of graphene with a single Tl adatom using a Monkhorst-Pack grid to sample the Brillouin zone. With the self-consistent charge density held fixed we perform non-self consistent calculations using a Monkhorst-Pack grid on 986 bands. The overlap matrices and projections are calculated using the PW2WANNIER90 routine of QUANTUM ESPRESSO. For the MLWF calculation we consider atom-centred functions on all C sites and , and functions on Tl (i.e. 70 Wannier functions in total accounting for spin). WANNIER90 is used to determine the MLWFs with the following parameters for disentanglement: frozen inner window: -4.0 to 2.0 eV, outer window -10.50 to 14.00 eV and maximum number of iterative steps: 2000. Typical Wannier function spreads for the optimised orbitals are Å (C) and Å(Tl).

As shown in Fig. 5 the MLWFs obtained by the procedure outlined above describe the band structure in very good agreement with the plane wave basis especially near the region of interest (i.e., close to the Dirac point). Hamiltonian matrix elements between the atom-centred MWLFs are computed and used to obtain an initial guess for the parameters in our minimal model. The first, second and third nearest neighbor C-C hopping parameters are computed for a pristine graphene supercell (without the Tl adatom) as appropriate for the dilute limit. Table I summarises the parameters extracted from the Wannier fit to the DFT calculation as well as those obtained by fitting directly to the DFT minimal model (see next section). (the local shift of the -orbital energies on C atoms neighboring Tl relative to pristine graphene) is estimated by taking the difference between the site energies of first and third nearest C atoms to the Tl adatom. In addition, our Wannierization calculations show that spin–dependent hopping integrals between Tl and the adjacent C atoms are negligible, justifying the neglect of those terms in the minimal model (see next section).

DFT Wannier | -1.07 | -0.59 | 0.37 | 0.19 | 1.65 | 1.28 | 2.87 | -0.24 | 0.26 | 0.26 |

### b.2 Minimal graphene–adatom model

Following Weeks et al. QSH_G_Adatom_Weeks_11 (), we parameterize
a tight-binding Hamiltonian informed by the first-principles calculations.
The Hamiltonian is separated into a graphene contribution ()
and a local adatom () and hybridization () terms.
The accurate description of bands () requires a third-nearest
neighbor approximation supp12_W_Jung_13 (). For the graphene–adatom
interaction it suffices considering the nearest neighbors [Sec. (B.3)].
Time-reversal symmetry and invariance with respect to the point group
constrain the Hamiltonian to the following form^{2}^{2}2Symmetry also allows spin-dependent tunneling processes, e.g.,
and . However, these corresponding
hopping integrals are extremely small (Sec. B.1),
and thus such terms can be safely omitted.:

(7) | ||||

(8) | ||||

(9) |

There are parameters in (9): the graphene hoppings , and , the chemical potential change on C atoms next to the Tl , the Tl p-state energies and , the Tl spin-orbit parameters (we assume ) and the C-Tl hoppings and . We adjust these parameters until the tight-binding band structures reproduce the ab initio DFT and GW results. Our final parameter values are given in Table 3 and the corresponding band structures are shown in Fig. 4(b).

DFT minimal | DFT Wannier | GW minimal | |
---|---|---|---|

0.83 | 0.26 | 0.28 | |

2.95 | 2.87 | 3.36 | |

-0.12 | -0.24 | 0.10 | |

-0.25 | 0.26 | -0.19 | |

1.59 | 1.65 | 1.30 | |

1.45 | 1.28 | 0.82 | |

-1.27 | -1.07 | -2.36 | |

-0.59 | -0.59 | -0.44 | |

0.33 | 0.37 | 0.51 |

The DFT optimized parameters are in excellent agreement with the hopping integrals obtained through the Wannierization procedure (Table 3). This is remarkable given that the Wannierization scheme considers all hopping integrals between any pair of orbitals, whereas the minimal model [Eq. (7)-(9)] is truncated. The discrepancy in the values of is expected given that the minimal model restricts the on-site potential to the six C atoms adjacent to Tl, which then overestimates the depth of the potential well.

### b.3 Decimation of the adatom degrees of freedom

We use standard degenerate perturbation theory to derive an accurate graphene-only Hamiltonian parameterized by our first-principles calculations. The effective interaction induced by a single adatom on graphene is

(10) |

with , and where () spans the adatom (graphene) subspace of states with definite angular momentum SHE_G_Pachoud_14-1. From Eqs. (8)-(9) we find:

(11) |

where the states have been arranged in ascending total angular momentum, i.e., and . Evaluation of the effective Hamiltonian Eq. (10) results in:

(12) |

with

(13) |

and

(14) |

Changing Eq. (12) into the honeycomb site basis yields the effective graphene–only coupling Hamiltonian reported in the main text (omiting on site terms):

(15) |

where

(16) |

The total effective tight-binding Hamiltonian for a graphene sheet with an adatom at the -th hollow site (restoring on-site terms) reads

(17) |

The on-site energy can be trivially absorbed in the definition of the local chemical potential in Eq. (7). Numerical studies showed that moderate values of do not impact the topological phase QSH_G_Adatom_Weeks_11 (). The underlying reason is that uniform on-site energies lead to a scalar (intra-valley) potential, to leading order in , in the long wavelenght description (see main text). The main effect of a moderate on site energy correction is thus a trivial shift in the spectrum: , where is the adatom coverage.

## Appendix C Real-Space Quantum Transport Calculations

### c.1 Methodology

We consider an adatom-decorated graphene nanoribbon connected to pristine
armchair-graphene leads^{3}^{3}3The popular choice of armchair termination for studies of the QSH
phase in graphene is meant to eliminate the possibility of edge states
occuring already in the absence of SOC (albeit with a different origin).
Our armchair nanoribbons are made deliberately metallic (;
with ) and hence their conductance exibit a natural
plateau at low energies of width associated with
states delocalized through the nanoribbon. This confinement effect
puts a constraint on the smallest SOC gaps that can be probed by the
transport calculations: (Refs. supp13_GNR_Son_06 (); supp14_GNR_Zheng07 ()). as depicted in Fig. 6. To determine the quantum
transport characteristics, we employ the Landauer-Büttiker formalism supp15_Buttiker_98 (); supp16_Landauer_70 (),
where the longitudinal conductance is expressed as supp17_Caroli_71 (); supp18_Datta_Book (); supp19_Fisher_81 ():

(18) |

where is Green’s function connecting the left and right interfaces of the central device and is the self-energy of the left (right) contact in the retarded (R)—advanced (A) sector. The central device with adatoms is modelled according to the following Hamiltonian:

(19) |

The calculation is carried out in two steps. Firstly, the pristine semi-infinite armchair graphene-contacts are treated analytically. Secondly, an efficient, highly-parallelizable recursive numerical approach is employed to extract the Green’s function of the central device. The method consists of subdividing the central region in smaller parts, whose Green’s functions can be evaluated via direct inversion supp20_Lewenkopf_13 (). The numerical evaluation of the Green’s function for the central region, together with the analytic treatment of the leads, allows to efficiently compute the conductance of sizable graphene nanoribbons with several millions of atoms.

### c.2 Bond current maps

The bond current density for up and down spins () is obtained from the variation rate of the charge density at a given site, which is written in terms of the flux of electrons through all of its bonds. For that, one needs to calculate the variation of the expectation value of the number operator over time, leading to the following expression for the current flowing from site to supp21_AntiPekka ():

(20) |

where is the correlation function and is the hopping between sites and . The density maps borne out the crucial role played by the new adatom-induced hoppings studied here (Fig. 7).

## Appendix D Real-Space Chern Number Calculations

### d.1 Methodology

We calculate the first Chern number using a definition which is valid for systems which are not translationaly invariant. To that purpose we use twisted boundary conditions,

(21) |

with for the twist angles , and for the system Halmiltonian with twisted boundary conditions , where is a basis vector of the direct lattice with lattice points in the direction of . The ground state wave function is obtained from single particle states ,

(22) |

with , and we assume that only single particle states are occupied.

If the ground state is non-degenerate and there is a finite energy gap between the ground state energy and the excited states in the bulk then the Chern number is well defined,

(23) |

with

(24) |

where is the surface . To evaluate Eq. (23), we apply Fukui’s method Chern_Fukui_05 () by discretizing the surface into points and then suming over the flux of individual plaquettes,

(25) |

where is a vector in the direction with magnitude , with the area of the surface . The advantage of this approach is that is strictly an integer for arbitrary lattice spacing, and it should rapidly converge to as one increases the size of the lattice . We adopted the implementation discussed in Ref. Chern_Zhang_13 ().

### d.2 Results

In order to assess quantitatively the transition between topologically trivial and non-trivial ground states, we have carried on a detailed calculation of the topological index as a function of adatom coverage. The results are shown in Fig. 8 from the lowest densities we can simulate up to 10% coverage. Since we calculate the Chern number on a finite size lattice with = unit cells, the lowest density corresponds to a single adatoms and takes the value .