PAPER • OPEN ACCESS You may also like - Realization of versatile electronic, Tuning the topological band gap of bismuthene magnetic properties and new topological phases in hydrogenated bismuthene with silicon-based substrates Ming-Yang Liu, Qing-Yuan Chen, Chao Cao et al. - Tunable electronic properties in To cite this article: Nils Wittemeier et al 2022 J. Phys. Mater. 5 035002 bismuthene/2D silicon carbide van der Waals heterobilayer Joy D. Sarker, Md. Sherajul Islam, Naim Ferdous et al. - Anisotropic thermoelectric effect on View the article online for updates and enhancements. phosphorene and bismuthene: first- principles calculations based on nonequilibrium Green’s function theory Yuto Tanaka, Mineo Saito and Fumiyuki Ishii This content was downloaded from IP address 158.109.48.66 on 02/11/2022 at 14:49 J. Phys. Mater. 5 (2022) 035002 https://doi.org/10.1088/2515-7639/ac84ad Journal of Physics: Materials PAPER Tuning the topological band gap of bismuthene with silicon-based OPEN ACCESS substrates RECEIVED 4 May 2022 Nils Wittemeier1,∗, Pablo Ordejón1 and Zeila Zanolli2 REVISED 1 7 July 2022 Catalan Institute of Nanoscience and Nanotechnology—ICN2 (CSIC and BIST), Campus UAB, Bellaterra, 08193 Barcelona, Spain 2 Chemistry Department, Debye Institute for Nanomaterials Science, Condensed Matter and Interfaces, Utrecht University and ACCEPTED FOR PUBLICATION European Theoretical Spectroscopy Facility, PO Box 80.000, 3508 TA Utrecht, The Netherlands 27 July 2022 ∗ Author to whom any correspondence should be addressed. PUBLISHED 9 August 2022 E-mail: nils.wittemeier@icn2.cat Keywords: bismuthene, topological insulator, 2D materials, quantum spin Hall, density functional theory, first principles methods Original Content from this work may be used Supplementary material for this article is available online under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution Abstract of this work must maintain attribution to Some metastable polymorphs of bismuth monolayers (bismuthene) can host non-trivial the author(s) and the title topological phases. However, it remains unclear whether these polymorphs can become stable of the work, journal citation and DOI. through interaction with a substrate, whether their topological properties are preserved, and how to design an optimal substrate to make the topological phase more robust. Using first-principles techniques, we demonstrate that bismuthene polymorphs can become stable over silicon carbide (SiC), silicon (Si), and silicon dioxide (SiO2 ) and that proximity interaction in these heterostructures has a significant effect on the electronic structure of the monolayer, even when bonding is weak. We show that van der Waals interactions and the breaking of the sublattice symmetry are the main factors driving changes in the electronic structure in non-covalently binding heterostructures. Our work demonstrates that substrate interaction can strengthen the topological properties of bismuthene polymorphs and make them accessible for experimental investigations and technological applications. 1. Introduction Topological materials (TM) are among the most promising candidates for next-generation nanoscale electronic devices and can host a variety of unconventional phenomena, such as quantum spin Hall states, Majorana, Weyl, or Dirac fermions [1]. TMs are characterized by material properties invariant under topological transformations and include topological insulators (TIs), topological semi-metals, and topological superconductors. TIs, in particular, can be exploited in applications due to their ability to conduct currents through edge states without dissipation [2–4]. However, the ability to generate robust TIs and control their band gap remains a critical issue for applications. Elements with strong intrinsic spin–orbit coupling (SOC), like Bismuth, are likely to be good TI candidates. Atomically thin Bi(111) layers are topological [5–9] but present an electronic band gap (0.08 eV) [9] too small for room-temperature operation. A possible solution is to find a substrate that enhances the gap while preserving the topological properties. Substrates can also induce a TI phase in topologically trivial monolayers: flat hexagonal (f-hex) bismuth monolayer is trivial in the free-standing form but becomes topological when grown on SiC(0001) [10], or partially passivated Si(111) [8]. Singh et al [9] recently predicted that there are other Bi monolayer phases besides the ground state puckered monoclinic (Bi(110)): buckled hexagonal (b-hex, Bi(111)), α, β, γ, f-hex, in order of formation energy. In free-standing form, only the b-hex and γ phases are topologically non-trivial [9]. The interaction with a suitable substrate could stabilize these metastable phases. However, only the b-hex and f-hex present the in-plane hexagonal symmetry as graphene and, therefore, can be Haldane-type TIs. In this paper, we explore the effect of silicon carbide (SiC), silicon (Si), and silicon dioxide (SiO2 ) substrates on the topological features of the two hexagonal phases of bismuthene. Thin layers of Bi have been synthesized on SiC [10] and Si [11]. The SiO2 substrate was instrumental for the experimental observation of © 2022 The Author(s). Published by IOP Publishing Ltd J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al mechanically exfoliated graphene [12] and has been widely used as a dielectric medium in integrated circuits. We demonstrate that buckled hexagonal bismuthene binds to silicon-based substrates via van der Waals interactions and retains its topological phase. F-hex bismuthene undergoes a trivial-to-topological phase transition only when placed on top of SiC. We further demonstrate that the interaction between b-hex and SiO2 is highly dependent on the surface structure: a hydroxylated SiO2 surface can enhance the topological band gap, and a cleaved silicon terminated (SiT) surface can reduce the band gap significantly. 2. Methods We carry out simulations using density functional theory (DFT) as implemented in the SIESTA code [13, 14], including SOC in the fully relativistic pseudopotential formalism [15, 16]. We employ fully relativistic, optimized norm-conserving Vanderbilt pseudopotentials [17] in PSML format [18] from PseudoDojo database [19] generated with PBE [20] exchange-correlation functional. The Kohn–Sham equations are solved using standard double-zeta polarized basis sets, a real space grid with a cut-off of 600 Ry, an electronic temperature of 5 meV, and a 15 × 15 × 1 Monkhorst–Pack grid for a single bismuthene unit cell. For larger cells, the k-point sampling is scaled accordingly. All structures are relaxed with a force threshold of 0.01 eV Å and a maximum stress tolerance of 0.006 eV Å3 . For the isolated substrate slabs, we employ the PBE functional [20] and, for the heterostructures, the van der Waals density functional of Dion et al [21] (vdW-DRSSL). We employ both the PBE and the vdW-DRSSL functional for the free-standing bismuthene phases and compare the difference in the resulting electronic structure. The use of a vdW functional is crucial for the heterostructures because GGA-functionals systematically underestimate binding energies [22]. The GGA-functional predicts that all heterostructures are not binding, except f-hex@SiC(0001). We use the vdW-DRSSL to avoid this systematic error and predict the binding energy of the vdW heterostructure more accurately. We ensure the decoupling of the upper and lower slab surfaces by imposing a 1 meV convergence threshold on the change in surface energy for any additional layer. We use a dipole correction to prevent exaggeration of the slab dipole due to interaction between periodic images. The Z2 topological invariant is determined using hybrid Wannier charge centers with Z2Pack [23]. We fix the lattice vectors of the substrate and strain the bismuthene phase. This approach recreates experimental conditions under which the monolayer will adjust to the more rigid substrate. We apply a counterpoise correction for the basis set superposition error [24] to evaluate the binding energies (EB ) between the monolayer and the substrate. 3. Results and discussion 3.1. Free-standing bismuthene The flat hexagonal phase of bismuthene (figures 1(a) and (b)) refers to the arrangement of bismuth atoms in a flat honeycomb lattice (space group: P6/mmm). We relax the structure and determine the length of the Bi–Bi bonds in this phase to be 3.09 Å, which corresponds to a lattice constant of 5.35 Å. We find that this phase of bismuthene is a semiconductor with an indirect band gap of 0.51 eV and trivial topology (Z2 = 0). In contrast, the buckled hexagonal phase (figures 1(c) and (d)) possesses a significantly smaller GGA band gap (0.13 eV) and features a non-trivial topology (Z2 = 1). In this phase, the Bi–Bi bonds are slightly shorter (3.04 Å) and form an angle with the lattice vectors. As a result, the atoms of the two sublattices are arranged in two parallel planes with a distance of 1.78 Å (buckling height). The lattice constant shrinks to 4.28 Å. In free-standing form, the buckled hexagonal phase is more stable than the flat hexagonal one (∆E = 0.4 eV). We observe a slightly larger band gap than Singh et al (0.08 eV [9]). When considering the use of different DFT codes, pseudopotentials, and basis sets, this degree of deviation in the band gap is reasonable. The band gap further increases (0.38 eV) when we use the van der Waals (vdW) functional instead of PBE. In the presence of vdW interaction, the buckling height in the b-hex phase decreases (1.65 Å) and small in-plane magnetic moments (0.01 µB ) emerge. These magnetic moments break the time-reversal symmetry and cause the splitting of the Kramer pairs (figure 1(f)). We determine the Z2 invariant and found that the topological phase of b-hex bismuthene is robust enough to persist despite the small spontaneous magnetization. 3.2. SiC(0001) In this section, we show that both metastable hexagonal bismuthene phases bind to the H-passivated SiC(0001) substrate, become stable (the structures could be experimentally synthesized), and are robust TIs. F-hex bismuthene can be grown epitaxially on SiC(0001) and becomes topologically√nontrivial √ [10]. The f-hex@SiC(0001) heterostructure consists of a 1% strained f-hex supercell on top of a ( 3 × 3)R30◦ 2 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Figure 1. Crystal structure of free-standing, flat hexagonal bismuthene (top view (a); side view (b)) and buckled hexagonal bismuthene (top view (d); side view (e)) and corresponding DFT electronic band structure (c), (f). supercell of SiC(0001) with silicon termination (figures 2(a) and (b)). The Bi atoms bind to the surface Si atoms leaving only the Si atom in the center of each Bi hexagon under-coordinated. which is passivated with a hydrogen atom. The monolayer binds to the substrate with a binding energy of 1.5 eV per Bi atom. The Bi–Si bond primarily involves pz orbitals (figure 3) and effectively removes the Bi-pz orbitals from the bands near the Fermi level. This removal of the pz orbitals drives the phase transition from trivial to topological [8, 10]. Our simulations reproduce the structural parameters and electronic structure reported by Reis et al [10]. We confirm the nontrivial topology of this combined structure and the indirect band gap, with the maximum of the valence bands at K and the minimum of the conduction bands of Γ. We find only small differences in the electronic structure calculated with the GGA-PBE and vdW-DRSSL functionals. For both functionals, we find a band gap of 0.6 eV, which is slightly larger compared to Reis et al [10] (0.5 eV). We attribute the small difference between vdW and GGA functional to the fact that the f-hex monolayer covalently binds to the substrate. Epitaxial growth of b-hex bismuthene on SiC(0001) requires surface passivation due to the large lattice mismatch several Si surface atoms have unsaturated bonds, which results in an unstable heterostructure. B-hex bismuthene can bind to passivated SiC(0001) via van der Waals interactions. The latter introduces only weak perturbations to the electronic structure, and we expect it to preserve the topological phase. A passivated surface will favor the b-hex phase or the puckered monoclinic phase over the other bismuthene phases with a higher free-standing formation energy [9] because the van der Waals interactions are too weak to change the energetic ordering. This is different from the case of f-hex in [10], where bismuthene is grown directly on top of silicon-terminated SiC(0001). The Bi atoms form covalent bonds with the surface Si atoms, and only the Si atoms that do not participate in the binding are passivated. The b-hex@H-SiC(0001) heterostructure consists of a −3% strained 3 × 3 b-hex supercell and a 4 × 4 supercell of SiC(0001) with silicon termination (figures 4(a) and (b)). The silicon atoms on both sides of the surface are passivated with hydrogen. As a result, the bonding between the monolayer and the substrate is weak (binding energy per Bi atom EB = 0.08 eV) and results in a van der Waals heterostructure. The proximity interaction in b-hex@H-SiC(0001) leads to a slight deformation of the first SiC layer, despite the weak binding energy. We computed the Z2 invariant and find that the topological phase of b-hex bismuthene is preserved on the passivated SiC substrate. The proximity interaction and the strain (SM figure 3) reduce the band gap (0.18 eV) with respect to the free-standing case (0.38 eV) and create a more pronounced Mexican-hat profile of the highest valence band. The robustness of the topological phase makes SiC a reasonable candidate for the growth of b-hex bismuthene irrespective the slightly reduced band gap. 3 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Figure 2. Crystal structure of flat hexagonal bismuthene@SiC(0001) (top view (a); side view (b)). Color code and schematic of the individual unit cells are displayed in the inset. Orbital projected DFT electronic band structure (c) in each panel the contribution of different orbitals is proportional to the line width. Figure 3. Flat hexagonal bismuthene@SiC(0001): change in charge density due to interaction between monolayer and substrate. The shape of the isosurface reveals, that bonds between Bi and Si are primarily formed by pz orbitals. The teal and yellow colors indicate a decrease and increase in the charge density respectively. 4 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Figure 4. Crystal structure of buckled hexagonal bismuthene@SiC(0001) (top view (a); side view (b)). Color code and schematic of the individual unit cells are displayed in the inset. Orbital projected DFT electronic band structure (c) in each panel the contribution of different orbitals is proportional to the line width. 3.3. SiO2 (0001) α-quartz Silicon dioxide occurs as an amorphous material, making it difficult to simulate with DFT methods that rely on periodic boundary conditions. We circumvent this limitation by simulating and comparing bismuthene monolayers on three different surface terminations of SiO2 α-quartz(0001): cleaved SiT, reconstructed oxygen terminated (ROT), and hydroxylated silicon terminated surface (OHT) (SM figure 1). We omit a simply oxygen terminated surface due to its structural instability. B-hex bismuthene binds to SiO2 via van der Waals interactions like SiC and becomes stable (could be experimentally synthesized). The topological phase of b-hex bismuthene is preserved. However, the band gap is highly sensitive to the surface structure of the SiO2 substrate. The hydroxylated silicon-terminated surface enhances the band gap, while a cleaved silicon-terminated surface reduces the band gap significantly. √ √ The b-hex@SiO2 heterostructures consist of a 2 × 2 b-hex supercell on top of a ( 3 × 3)R30◦ substrate supercell, imposing 2% (figure 5). B-hex bismuthene weakly binds to the three SiO2 surfaces with binding energies ⩽0.11 eV per bismuth atom (table 1). b-hex@OHT-SiO2 possesses the largest band gap (0.40 eV, figure 5(c)) among the three heterostructures, making it the most promising candidate for the realization of quantum spin Hall (QSH) states. The band gap is increased compared to that of the free-standing b-hex monolayer, and a rehybridization of the states near the Fermi level occurs. In the free-standing monolayer, all three Bi-p orbitals contribute equally to valence and conduction bands. In b-hex bismuthene@OHT-SiO2 the symmetry D3d of the free-standing b-hex bismuthene is broken (the heterostructure point group is P1) lifting the px –py symmetry. The Bi-px and Bi-py orbitals contribute mainly to the conduction and valence bands, respectively. As result, the valence band minimum of b-hex at Γ is pushed to slightly higher energies and the valence bands near M and K are lowered (figure 5(c)) and the direct band gap at Γ is slightly increased. The hybridization reveals a clear signature of band inversion at Γ, aside from increasing the size of the topological gap. These two characteristics indicate a 5 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Figure 5. Crystal structure of buckled hexagonal bismuthene on SiO2 (top view and side view) and orbital projected DFT electronic band structure for three different surface terminations: hydroxylated (a)–(c), reconstructed (d)–(f), and Si terminated (g)–(i). In each panel the contribution of different orbitals is proportional to the line width. On the hydroxylated and reconstructed surfaces the signature of a band inversion is visible in the projection of Bi 6py orbitals. Bi 6s orbitals do not contribute to the bands near the Fermi level. very robust topological phase of b-hex@OHT-SiO2 . We find that the Z2 invariant of b-hex@OHT-SiO2 is 1 (non-trivial). The b-hex@ROT-SiO2 displays a similar behavior: the topological phase is preserved, interaction with the substrate causes a rehybridization of the states near the Fermi level, and the signature of a band inversion is revealed (figure 5(f)). However, the band gap remains the same as in the free-standing b-hex bismuthene phase (0.38 eV). In the b-hex@SiT-SiO2 heterostructure the spacing between monolayer and substrate is significantly smaller (table 1) causing a stronger interaction between the two materials. In 6 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Table 1. Summary of basic properties of buckled hexagonal and flat hexagonal Bi monolayer phases in free-standing format and supported by different substrates. The positive positive binding energies indicate that the heterostructures can bind. The topological properties of the b-hex phase are preserved independently of the strain and the details of the interaction with the substrate. Similarly, the f-hex phase remains trivial, except on SiC where is covalently binds to the surface. Lattice Lattice EB per Interlayer Bi phase Substrate Z2 constant (Å) mismatch Bi (eV) distance (Å) Eg (eV) b-hex Free-standing 1 4.28 — — — 0.38 SiC(0001) Si term. & 1 3.12 3% 0.08 2.73 0.18 H passivated Si(111) H passivated 1 3.88 4% 0.07 2.66 0.12 SiO2 (0001) Hydroxylated 1 5.03 −2% 0.06 2.90 0.40 SiO2 (0001) Reconstructed 1 5.03 −2% 0.07 2.96 0.38 SiO2 (0001) Si terminated 1 5.03 −2% 0.11 2.34 0.08 f-hex Free-standing 0 5.35 — — — 0.51 SiC(0001) Si terminated 1 3.12 1% 1.49 2.76 0.62 SiO2 (0001) Hydroxylated 0 5.03 −2% 0.11 3.22 0.19 SiO2 (0001) Reconstructed 0 5.03 −2% 0.09 2.76 0.80 SiO2 (0001) Si terminated 0 5.03 −2% 0.19 — — the relaxed structure, the Bi-hexagons are skewed and the buckling height varies between the atomic sites (figure 5(g)). The proximity to the substrate also causes hybridization between Bi-6p and Si-3p orbitals. As a result of rehybridization and symmetry break, the band gap decreases and the splitting of the valence bands increases (figure 5(i)). The computed Z2 invariant shows that the topological phase of b-hex bismuthene is retained despite the distortion of the hexagonal structure and the hybridization with the substrate. The reduction of the band gap of b-hex on SiT-SiO2 prevents its application for quantum spin Hall devices. We predict that the optimal experimental realization of topological b-hex requires a hydroxylated SiO2 surface. F-hex bismuthene also binds to the SiO2 surfaces: weakly to the OHT-SiO2 and ROT-SiO2 surfaces (van der Waals), more strongly to SiT-SiO2 . The latter distorts the monolayer and breaks its symmetry. The topologically trivial phase of the f-hex monolayer is retained on all three surfaces. The computed band structures are reported in the supplementary materials (SM figure 2). 3.4. Si(111) The topologically non-trivial phase of b-hex bismuthene binds to Si(111) and becomes stable. However, the proximity interaction slightly lowers the band gap, making Si(111) a poor choice for TI devices based on b-hex bismuthene. √ √ The b-hex@Si(111) heterostructure consists of a 4% strained ( 3 × 3)R30◦ b-hex supercell and 2 × 2 Si(111) supercell (figures 6(a) and (b)), with a binding energy of 0.07 eV per Bi atom. Passivation of the Si(111) substrate is required to prevent the heterostructure from being metallic and results in a van der Waals heterostructure. Similar to the b-hex@SiO2 heterostructures the Si(111) substrate breaks the symmetry of the monolayer and as a result the orbital contribution of Bi-px and Bi-py to the bands are not symmetric. In this heterostructure no hybridization between the monolayer and the substrate occurs and the electronic structures of the two components are largely independent. The Fermi levels in this heterostructure are aligned such that the valence band minimum of bismuthene lies just above the conduction band minimum, and the band gap closes. However, the topological band gap should not be measured between the conduction band minimum of the substrate and the top of the Bi valence, because of the independence of the components. Instead, the relevant gap for the realization of quantum spin Hall states should be measured between Bi bands: 0.12 eV (figure 6(c)). This gap is significantly reduced compared to the free-standing phase due to strain and substrate interaction. This makes Si a poor choice for a substrate for b-hex bismuthene. 3.5. Discussion All bismuthene monolayers form van der Waals heterostructures with the silicon-based substrates, except f-hex@SiC(0001) in which covalent binding occurs. Despite the interaction of strain and substrate, the topological properties of the b-hex phase are preserved on all substrates. Similarly, the f-hex phase remains trivial. We demonstrate that all heterostructures have positive binding energies (table 1) but we do not predict thermodynamic stability. Nevertheless, our work can help guide experimental studies in choosing potential substrates. While our simulations of bismuthene on SiO2 apply directly to α-quartz, amorphous SiO2 (a-SiO2 ) is a more commonly used substrate. Our results could still apply to a-SiO2 despite a non-crystalline surface because (a) the local atomic structures of a-SiO2 are very similar to the α-quartz, and (b) the roughness of a-SiO2 is small (σ = 2 Å), and (c) monolayers can follow the roughness quite closely [25]. 7 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Figure 6. Crystal structure of buckled hexagonal bismuthene@passivated Si(111) (top view (a); side view (b)). DFT electronic band structure (c) in each panel the contribution of different orbitals is proportional to the line width. To ultimately observe the quantum spin Hall effect, the presence of geometrical edges of the topological monolayers is required. For example, in [10] the topological properties are unperturbed by the presence of the edges, despite the covalent bonding between the substrate and the monolayer. This has also been seen in [26, 27]. Therefore, we expect that the topological character of the heterostructures presented here will be preserved if an edge is introduced in the bismuth monolayer. 4. Conclusion In this work, we study the proximity interaction between hexagonal Bi monolayer phases and silicon-based substrates to identify potential substrates for room-temperature TI applications based on Bi monolayers. We demonstrate that buckled hexagonal and flat hexagonal bismuthene could be stabilized on H-SiC(0001), H-Si(111), and α-quartz SiO2 (0001). We show that the proximity interaction in the heterostructure has a significant effect on the electronic structure of the monolayer, even when no bonding occurs. We further identify the structural realignment and the breaking of the sublattice symmetry of the b-hex monolayer as the main factors driving the change in the electronic structure. The exact nature of this proximity interaction is sensitive to the substrate material and its surface composition. In particular, the size of the topological gap of b-hex bismuthene varies drastically depending on the chosen substrate. Importantly, the topology of the hexagonal Bi monolayers (f-hex: trivial, b-hex: non-trivial) is unaffected by the interaction with the substrate, except f-hex@SiC(0001). Hydroxylated SiO2 emerges as an especially good substrate choice for b-hex bismuthene due to an increase in the topological gap and could enable room-temperature application. Data availability statement The data that support the findings of this study are available upon reasonable request from the authors. 8 J. Phys. Mater. 5 (2022) 035002 N Wittemeier et al Acknowledgments We acknowledge the CERCA programme of the Generalitat de Catalunya (Grant 2017SGR1506), and by the Severo Ochoa programme (MINECO, SEV-2017-0706). Z Z acknowledges financial support by the Ramon y Cajal program (RYC-2016-19344), the Netherlands Sector Plan program 2019–2023, and support from the Dutch Gravity program “Materials for the Quantum Age (QuMat)”. P O acknowledges support from Spanish MICIU, AEI and EU FEDER (Grant No. PGC2018-096955-B-C43) and the European Union MaX Center of Excellence (EU-H2020 Grant No. 824143). N W acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 754558, and the ICN2 Severo Ochoa Outbound Mobility Programme. The work has been performed under the Project HPC-EUROPA3 (INFRAIA-2016-1-730897), with the support of the EC Research Innovation Action under the H2020 Programme. We acknowledge computing resources on MareNostrum4 at Barcelona Supercomputing Center (BSC), provided through the PRACE Project Access (OptoSpin Project 2020225411) and RES (Activity FI-2020-1-0014), resources of SURFsara the on National Supercomputer Snellius (EINF-1858 Project) and technical support provided by the Barcelona Supercomputing Center. ORCID iDs Nils Wittemeier  https://orcid.org/0000-0002-1437-8659 Pablo Ordejón  https://orcid.org/0000-0002-2353-2793 Zeila Zanolli  https://orcid.org/0000-0003-0860-600X References [1] Giustino F et al 2020 J. Phys. Mater. 3 042006 [2] Hasan M Z and Kane C L 2010 Rev. Mod. Phys. 82 3045–67 [3] Qi X L and Zhang S C 2011 Rev. Mod. Phys. 83 1057–110 [4] Moore J E 2010 Nature 464 194–8 [5] Koroteev Y M, Bihlmayer G, Chulkov E V and Blügel S 2008 Phys. Rev. B 77 045428 [6] Liu Z, Liu C X, Wu Y S, Duan W H, Liu F and Wu J 2011 Phys. Rev. Lett. 107 136805 [7] Ma Y, Li X, Kou L, Yan B, Niu C, Dai Y and Heine T 2015 Phys. Rev. B 91 235306 [8] Zhou M, Ming W, Liu Z, Wang Z, Yao Y and Liu F 2015 Sci. Rep. 4 7102 [9] Singh S, Zanolli Z, Amsler M, Belhadji B, Sofo J O, Verstraete M J and Romero A H 2019 J. Phys. Chem. Lett. 10 7324–32 [10] Reis F, Li G, Dudy L, Bauernfeind M, Glass S, Hanke W, Thomale R, Schäfer J and Claessen R 2017 Science 357 287–90 [11] Kawakami N, Lin C L, Kawai M, Arafune R and Takagi N 2015 Appl. Phys. Lett. 107 031602 [12] Novoselov K S, Geim A K, Morozov S V, Jiang D, Zhang Y, Dubonos S V, Grigorieva I V and Firsov A A 2004 Science 306 666–9 [13] Soler J M, Artacho E, Gale J D, García A, Junquera J, Ordejón P and Sánchez-Portal D 2002 J. Phys.: Condens. Matter 14 2745–79 [14] García A et al 2020 J. Chem. Phys. 152 204108 [15] Cuadrado R and Cerdá J I 2012 J. Phys.: Condens. Matter 24 086005 [16] Cuadrado R, Robles R, García A, Pruneda M, Ordejón P, Ferrer J and Cerdá J I 2021 Phys. Rev. B 104 195104 [17] Hamann D R 2013 Phys. Rev. B 88 085117 [18] García A, Verstraete M J, Pouillon Y and Junquera J 2018 Comput. Phys. Commun. 227 51–71 [19] van Setten M, Giantomassi M, Bousquet E, Verstraete M, Hamann D, Gonze X and Rignanese G M 2018 Comput. Phys. Commun. 226 39–54 [20] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett. 77 3865–8 [21] Dion M, Rydberg H, Schröder E, Langreth D C and Lundqvist B I 2004 Phys. Rev. Lett. 92 246401 [22] Lazić P, Atodiresei N, Alaei M, Caciuc V, Blügel S and Brako R 2010 Comput. Phys. Commun. 181 371–9 [23] Gresch D, Autès G, Yazyev O V, Troyer M, Vanderbilt D, Bernevig B A and Soluyanov A A 2017 Phys. Rev. B 95 075146 [24] Boys S and Bernardi F 1970 Mol. Phys. 19 553–66 [25] Quereda J, Castellanos-Gomez A, Agrait N and Rubio-Bollinger G 2014 Appl. Phys. Lett. 105 053111 [26] Wu R et al 2016 Phys. Rev. X 6 021017 [27] Bieniek M, Woźniak T and Potasz P 2017 J. Phys.: Condens. Matter 29 155501 9