Heterostructures of graphene and topological insulators Bi$_2$Se$_3$, Bi$_2$Te$_3$, and Sb$_2$Te$_3$

Prototypical three-dimensional topological insulators of the Bi$_2$Se$_3$ family provide a beautiful example of the appearance of the surface states inside the bulk band gap caused by spin-orbit coupling-induced topology. The surface states are protected against back scattering by time reversal symmetry, and exhibit spin-momentum locking whereby the electron spin is polarized perpendicular to the momentum, typically in the plane of the surface. On the other hand, graphene is a prototypical two-dimensional material, with negligible spin-orbit coupling. When graphene is placed on the surface of a topological insulator, giant spin-orbit coupling is induced by the proximity effect, enabling interesting novel electronic properties of its Dirac electrons. We present a detailed theoretical study of the proximity effects of monolayer graphene and topological insulators Bi$_2$Se$_3$, Bi$_2$Te$_3$, and Sb$_2$Te$_3$, and elucidate the appearance of the qualitatively new spin-orbit splittings well described by a phenomenological Hamiltonian, by analyzing the orbital decomposition of the involved band structures. This should be useful for building microscopic models of the proximity effects between the surfaces of the topological insulators and graphene.

In contrast to the topological insulators, the Dirac states in graphene are not topologically protected. However, due to the two-dimensional nature of graphene, one can easily manipulate its electronic states via so called proximity effects [13]. Within van der Waals heterostructures [14,15,16] with other two dimensional materials, one can induce magnetism, as well as strong spin-orbit coupling (SOC) in graphene [17,18,19,20]. Combining proximity-induced exchange and SOC in graphene can, under the right conditions, also lead to topologically protected edge states [21,22], which could be important for novel spintronics applications.
There have already been numerous studies considering graphene/topological insulator bilayers [23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39], in which two different kinds of Dirac electrons are simultaneously present. More specifically, it is possible to grow highquality topological insulators such as Bi 2 Se 3 epitaxially layer-by-layer on graphene, with small defect density [29,30]. Such heterostructures can actually be used as broadband photodetectors [36]. Relevant to our work are spin properties of graphene/topological insulator slabs. It has been argued that these structures can still exhibit quantum spin Hall states [31,35], and that the type and magnitude of proximity SOC in graphene can be tuned by the twist angle [26,23]. Finally, spin transport experiments have demonstrated spin-to-charge conversion in graphene on (Bi 0.15 Sb 0.85 ) 2 Te 3 [40] and on Bi 2 Te 2 Se [38]. Similar to transition-metal dichalcogenides [17,18], topological insulators strongly enhance the negligible intrinsic SOC of the graphene Dirac states from 10 µeV [11,41,42], by two orders of magnitude to about 1 meV [23,25,43]. The proximity-induced SOC in graphene is giant, drastically reducing spin relaxation times, and of valley-Zeeman type, leading to giant spin relaxation anisotropies [23,44]. Such graphene/topological insulator bilayers are ideal for the interaction of topological surface states, with in-plane spin-momentum locking, and the proximity induced spin-orbit fields in graphene.
In this manuscript, we first review the properties of Bi 2 Te 3 as a representative example of the three dimensional topological insulator family. Our first-principles results, considering 8 layers of Bi 2 Te 3 where the Dirac surface states have already formed, are consistent with literature. We include this background information to set the stage for the discussion of main results. Since we will be interested in the effects of an external (transverse) electric field, we also investigated the gate effect on the degenerate surface states in the Bi 2 Te 3 slab, which exhibit spinmomentum locking. Indeed, we show that these topological states can be efficiently separated in energy by an applied electric field. The splitting, ∆E, of the surface states, increases linearly with the slope of about 6.5 meV per mV/nm with the applied field. This DFT prediction agrees well with a simple estimate based on electrostatics.
In the second part, we consider graphene/topological insulator bilayers where we are interested in the proximityinduced SOC in graphene. We quantify the magnitude and type of induced SOC by fitting a symmetry-derived model Hamiltonian to the low energy bands of graphene. The proximity-induced SOC in graphene is similar, but with variations in magnitude (0.1-1 meV), for the considered topological insulators Bi 2 Se 3 , Bi 2 Te 2 Se, and Sb 2 Te 3 . Moreover, the charge transfer between materials and the resulting doping level of graphene can be significantly different for different topological insulators, ranging from 0 to 350 meV in terms of the Fermi energy. When the Dirac points of both materials are located near the Fermi level, as the case of Bi 2 Te 2 Se indicates, the simultaneous study of two very different spin-orbit fields is possible. Motivated by the recent spin-charge conversion experiments in graphene on (Bi 0.15 Sb 0.85 ) 2 Te 3 [40], we extensively discuss the case of graphene/Sb 2 Te 3 , including spin-orbit fields, and the gate tunability of proximity-induced SOC and the doping level. We find a giant electric field tunability of Rashba and intrinsic SOC, in magnitude and sign, which is important to interpret the above experimental data.
2 Monolayer graphene in proximity to Bi 2 Se 3 and Bi 2 Te 3 2.1 Topological band structure of Bi 2 Te 3 We begin by describing the topological band structure of Bi 2 Te 3 , in order to analyze the spin projection of the Dirac electrons as well as the orbital decomposition of the states.
For the calculation of the topological insulator Bi 2 Te 3 , we set up the atomic structure with the Atomic Simulation Environment (ASE) [45]. We consider 8 quintuple layers (QLs) of Bi 2 Te 3 , using the lattice constants [46] a = 4.386Å and c = 30.497Å, with the atomic parameters (u, v) = (0.4000, 0.2097). In Figure 1(a) we show the geometry of 8 QLs of Bi 2 Se 3 , visualized with VESTA [47]. The unit cell contains 40 atoms.
The electronic structure calculations are performed by density functional theory (DFT) [48] with Quantum ESPRESSO [49]. Self-consistent calculations are performed with the k-point sampling of 30 × 30 × 1. The energy cutoff for the charge density is 600 Ry, and the kinetic energy cutoff for the wavefunctions is 70 Ry We consider relativistic pseudopotentials with the projector augmented wave method [50] employing the Perdew-Burke-Ernzerhof exchange correlation functional [51]. Dipole and van der Waals corrections [52,53,54] are included to get correct band offsets and internal electric fields. In order to simulate the 8 QL slab of Bi 2 Te 3 , we add a vacuum layer of 30Å.
In Figs. 1(c-e), we show the calculated low energy band structure, projected on the three different parts (top QL, bottom QL, bulk) of the geometry, defined in Figure 1(a), along the k-path shown in Figure 1(b). We find that the Dirac states, that cross the Fermi level, are localized in the top and bottom QL layer of the Bi 2 Te 3 . Due to inversion symmetry of the 8 QL structure, the Dirac states of top and bottom QL are at the same energy, but with opposite spin. In Figs. 1(g-h), we show the same band structure where the color corresponds to the s x , s y , and s z spin expectation value, respectively. As the chosen k-path is along the k xdirection, the Dirac states only have a s y spin component.
In 1(f) we show the spin-orbit field of the top QL Dirac bands around the Γ point. As expected, the Dirac states show spin-momentum locking. Away from the center of the Brillouin Zone, the Dirac bands also show some trigonal warping. Still within the bulk gap the Dirac states acquire perpendicular (s z ) spin, which increases as the states get closer to the bulk bands.
Our calculated low energy band structure agrees very well with ARPES measurements [55] and earlier DFT results [56,1]. Especially the Dirac point of Bi 2 Te 3 is at about −150 meV below the Fermi level and located within the bulk bands, see Figure 1(c). In contrast, other topological insulators such as Bi 2 Se 3 have the Dirac point at the Fermi level under ideal defect free conditions [1]. However, from the experimental point of view, unintentional intrinsic doping is present for all members of the Bi 2 Se 3 topological insulator family, and the Dirac point is typically located below the Fermi level [55]. To compensate this effect and to bring the Dirac states to the Fermi level, the multicompositional topological insulator crystals Bi 2−x Sb x Te 3−y Se y are considered [57,58,25]. Depending on the numbers x and y, the defect doping can be counteracted and bulk transport can be suppressed. Of course, the formation of the Dirac states depends also on the number of QLs [59,8,9,10,60], because for thin samples the surface state wave functions still interact with each other through the bulk, such that gapless surface states are absent.
A possibility to break the aforementioned degeneracy of the Dirac states of top and bottom QL is by the application of a transverse electric field along the c-axis or through a substrate, breaking the inversion symmetry. Depending on the potential difference on the two sides, the Dirac states will be separated in energy [61]. The electric field that we apply is modeled by a sawtooth potential, and we can directly estimate the potential energy difference ∆V = e × A × d, from the electric field amplitude A, the thickness of the 8 QL structure d = 78.7Å (the distance between the outermost Te atoms), and e is the charge of the electron. In Figure 2(a), we show a zoom to the calculated Dirac surface states of Bi 2 Te 3 for a transverse electric field of 5 mV/nm. We find that the states originating from top and bottom QL are still intact and separated by ∆E in energy. The dipole of the structure, see Figure 2(b), grows linearly with the applied field. In Figure 2(c), we compare the energy splitting ∆E, extracted from the calculated band structures, with the estimated potential difference ∆V , as function of the applied electric field. Both depend linearly on the applied field, as expected, but the energy splitting ∆E is smaller than the estimated potential difference ∆V for all field values. This can be attributed to the fact that the surface states are localized within top and bottom QL and their spatial separation is not exactly equal to the thickness d, as we use in the estimation for ∆V . The splitting ∆E increases with a slope of roughly 6.5 meV per mV/nm of applied field.

Proximitized graphene: effective Hamiltonian with spin-orbit coupling
In order to understand the proximity effect on the electronic band structure of graphene, we first introduce a generic phenomenological model describing Dirac states of graphene with reduced symmetry due to external effects [62,63,64,20,17,18,65,66]. The model Hamiltonian given in the basis |Ψ A , ↑ , |Ψ A , ↓ , |Ψ B , ↑ , and |Ψ B , ↓ reads: The first term H 0 describes a gapless linear dispersion near Dirac points K (K') with two-fold spin-degenerate bands. The parameter v F denotes the Fermi velocity, and k x and k y are the Cartesian components of the electron wave vector measured from ±K, corresponding to the valley index τ = ±1. The Pauli spin matrices are s i and σ i are pseudospin matrices, with i = {0, x, y, z}. We also define σ ± = 1 2 (σ z ± σ 0 ) for shorter notation. The pristine graphene lattice constant is a.
When graphene is situated above a substrate, the pseudospin symmetry of graphene gets broken and H ∆ describes a mass term, opening a gap in the spectrum [67,68]. The corresponding parameter ∆ is called staggered potential and models the size of the induced gap. Of course, the pseudospin symmetry breaking depends on the interlayer distance [64,20,32] and on the actual arrangement of graphene above the substrate's surface which can be tuned by twisting [23,26]. However, the sublattice potential asymmetry is not always responsible for the gap opening. Also a Kekulé lattice distortion [69,70,71,23,72], leading to a nearest neighbor hopping asymmetry [73], and SOC, e. g., from adatoms [74,75,76], can open the band gap in graphene.
The Dirac bands of freestanding graphene show an intrinsic SOC of about 12 µeV [11,42,41]. Due to a substrate, also the SOC in the effective graphene p z orbitals, forming the Dirac bands, can be modified. The term H I accounts for the modification of the intrinsic SOC due to proximity effects, where λ A I and λ B I are the sublattice resolved intrinsic SOC parameters.
The presence of a transverse electric field (vertical to the graphene layer) or a substrate breaks all symmetries, that would allow to flip the orientation of the transverse z axis (inversion with respect to z or mirror with respect to the xy-plane). Two additional terms arise due to this symmetry breaking, namely H R and H PIA . The first term is the Rashba SOC with parameter λ R , which describes the amount of space inversion asymmetry. The second term is the sublattice resolved pseudospin-inversion asymmetry (PIA) SOC Hamiltonian with parameters λ A PIA and λ B PIA , which describe the strength of the mirror plane asymmetry.
Finally, E D accounts for electron or hole doping of the Dirac bands due to external influences and we call it the Dirac point energy.

Graphene/topological insulator van der Waals bilayers
We now realize, using atomistic simulations, the effective Hamiltonian introduced in the previous section by combining graphene with a single quintuplet of Bi 2 Se 3 , Bi 2 Te 2 Se, and Sb 2 Te 3 topological insulators. The resulting structure is essentially a van der Waals bilayer, as we show in Figure 3.
For the calculation of the graphene/topological insulator bilayers we consider a 5 × 5 supercell of graphene on top of a 3 × 3 supercell of a topological insulator. Initial atomic structures are set up with ASE [45] and the heterostructure was visualized with VESTA [47], see Figure  3. For periodic DFT calculations, we need to marginally strain the constituent layers in order to form a commensurate unit cell. Therefore, we strain the graphene lattice constant [77] to a = 2.486Å and use the lattice structure of Bi 2 Se 3 , according to Ref. [46], extracting only 1QL of

Figure 3
Geometry of graphene above one QL of (Bi/Sb) 2 (Se/Te) 3 . Different colors correspond to different atomic species, as in Figure 1.
the topological insulator. For the other topological insulators, Bi 2 Te 2 Se and Sb 2 Te 3 , we simply replace the relevant atoms without changing the geometry. We consider only bilayers without relaxation, using interlayer distances of 3.5Å between the graphene layer and the QL of the topological insulator [23,20]. The first-principles calculations are performed in a similar way as for the 8QL Bi 2 Te 3 structure, discussed above. For the bilayer structures, we use a k-point sampling of 9 × 9 × 1, an energy cutoff for the charge density of 500 Ry, and a kinetic energy cutoff for wavefunctions of 60 Ry. Dipole and van der Waals corrections are also included [52,53,54]. Moreover, a vacuum layer of 24Å is added, to avoid interactions between periodic images in our slab geometry.
Note that for describing the electronic structure of a 3D topological insulator, the GW method is a more accurate choice, as compared to the generalized-gradientapproximation (GGA) employed in this work [78,60,79,80]. However, as we can see in Figure 1, the GGA also captures the main band structure features of the topological insulator and matches the ARPES measurements [55]. This makes sense, since GW is usually employed to faithfully describe the orbital gap in semiconductors, while the topological surface states are gapless (up to finite-size hybridization) and GGA is fully adequate to describe them and yield reliable predictions. Furthermore, GW calculations are computationally very demanding and inaccessible for such large heterostructure systems we consider here. In addition, recent GGA-based calculation results [23] have already been successfully used for the interpretation of ex-perimental data for graphene/topological insulator structures [25].
In Figure 4 we show the calculated band structures of graphene on one QL of Bi 2 Se 3 , Bi 2 Te 2 Se, and Sb 2 Te 3 . In the case of Bi 2 Se 3 , the Dirac point energy E D is well above the Fermi level indicating strong hole doping, similar to Refs. [23,81]. In contrast, for the other two topological insulators, the Dirac point of graphene is located at the Fermi level. Since the topological insulator thickness is just 1QL, the surface states have not yet developed [7,8,9,10,20]. However, we indicate the topological insulator surface states with the energy E TI in Figure 4(a). By fitting the Hamiltonian H from the previous section to the graphene Dirac bands we can extract several relevant orbital and spin-orbit parameters. In Table 1 we show the fit parameters for the different bilayers. The accuracy of the fit is shown in the next section, where we analyze the graphene/Sb 2 Te 3 case in detail.
The Fermi velocity is roughly independent of the topological insulator substrate. The sublattice symmetry breaking of graphene due to the topological insulator, described by the staggered mass parameter ∆, is tiny and almost negligible compared to the other parameters. Consequently, the potential asymmetry of the graphene sublattices in our investigated structure is small. However, by twisting the layers [23,26] or by decreasing the interlayer distance between graphene and the topological insulator surface [64,20,32], the gap in graphene's spectrum can be enhanced.
Interestingly, the intrinsic SOC parameters λ A I and λ B I are almost equal in magnitude but opposite in sign for all the studied bilayers. Such a valley-Zeeman type SOC, i. e. λ A I ≈ −λ B I , can lead to giant spin-relaxation anisotropies in graphene [44,18,23]. A more detailed analysis of the graphene/Bi 2 Se 3 and graphene/Bi 2 Te 2 Se cases is given in Refs. [20,23]. More precisely, depending on the twist angle and the exact interface of graphene and the topological insulator, a giant spin relaxation anisotropy can be present [23]. Additional QLs of the topological insulator are necessary for the surface states to form, but will have a minor extra impact on graphene's band structure, since proximity effects are short-ranged. In contrast, two very efficient tunability knobs for band offsets and proximity SOC are the interlayer distance and a transverse electric field [20].
Different to monolayer graphene, bilayer graphene shows a giant band gap, due to the intrinsic dipole present in heterostructures with Bi 2 Se 3 . Moreover, the resulting proximity band structure of bilayer graphene can be tuned by gating and a spin-orbit valve can be realized [20,82].

Sb 2 Te 3 substrate
The proximity effect in graphene due to Sb 2 Te 3 has not yet been systematically studied. Below we provide both DFT results and phenomenological descriptions for these bilayers.
In Figure 5(b) we show the calculated band structure for the graphene/Sb 2 Te 3 heterostructure, in the absence of a transverse electric field applied across the bilayer structure. We find that the Dirac point of graphene, as well as the band edge originating from the topological insulator is lo-cated at the Fermi level. The overall band structure is comparable to ARPES measurements of graphene on a thick Sb 2 Te 3 substrate, showing the coexistence of both Dirac cones near the bulk Sb 2 Te 3 valence band edge [83].
When a negative transverse electric field of −2 V/nm is applied across the bilayer, see Figure 5 do not shift in energy, while the bands of the topological insulator do. In Figure 5(c), we also label the doping energy of the topological insulator with E TI , as these bands would correspond to the topological surface states in few QL structures. In Figure 6, we show the low energy band properties of the graphene Dirac states, fitted to the model Hamiltonian, for zero electric field. The model agrees perfectly with the DFT calculated band structure, capturing also the spin expectation values and band splittings, using the parameters summarized in Table 1 for the Sb 2 Te 3 substrate.
In Figure 7 we summarize the evolution of the fit parameters as function of a transverse electric field, applied across the bilayer. Most interesting are the intrinsic SOC parameters, which can be tuned from positive to negative values, but always of valley-Zeeman type. The Rashba and PIA SOC parameters are also strongly changing with the applied field and can be even tuned to zero. The resulting spin-orbit fields of the Dirac bands are due to a competition of Rashba and PIA SOC favoring an in-plane spin texture, and the intrinsic SOCs favoring an out-of-plane texture. Due to tunability of these parameters with the electric field, we have a potential knob to tune the spin-orbit fields, as well as the magnitude of the proximity-induced SOC. The spin-orbit fields of the four Dirac bands, as la-  beled in Figure 6(e), are shown in Figure 8 for the zero Copyright line will be provided by the publisher field case. We can see that bands show very pronounced Rashba spin-orbit fields. For example, the first conduction band (CB 1 ) shows counter-clockwise, while the second conduction band (CB 2 ) shows a clockwise rotating spinorbit field, both also with a significant and opposite out-ofplane spin component.
Recently, a gate-tunable spin-galvanic effect has been shown experimentally in graphene/topological insulator bilayers [40]. More precisely, they demonstrate an efficient spin-charge conversion at room temperature in graphene/(Bi 0.15 Sb 0.85 ) 2 Te 3 heterostructures, which should be well comparable to our graphene/Sb 2 Te 3 bilayers. Especially the electric field results in Figure 7 can be used to explain their gate-tunability of the conversion efficiency, due to tunable proximity SOC.
Based on the above results, we can conclude that for device applications, only a thin (1-2 QLs) topological insulator is sufficient to fully exploit it's proximity effect on graphene. A thicker topological insulator is necessary for the Dirac surface states to form, allowing to simultaneously study two types of Dirac electrons, with very different spin-orbit fields. An electric field can be used to tune both, the surface states of the topological insulator and the proximity SOC in graphene. The magnitude of proximity SOC and band offsets depend on the topological insulator crystal. Consequently, a multicompositional material Bi 2−x Sb x Te 3−y Se y might be the best choice for applications, since proximity effects can be maximized with energetically aligned Dirac states. Especially the mentioned gate tunable spin-charge conversion is important for novel spin-orbit technology, without the need of ferromagnets.

Summary
We have reviewed the basic properties of the topological insulator Bi 2 Te 3 and find gate tunable energy splitting of Dirac states, which results from the potential difference in the surface states. The energy splitting increases linearly with a slope of about 6.5 meV per mV/nm of applied field, which can be experimentally verified. Additionally, we have reported original results for graphene/Sb 2 Te 3 bilayers in the context of related graphene heterostructures with Bi 2 Se 3 and Bi 2 Te 2 Se. We find that the position of the graphene Dirac point strongly depends on the substrate, when considering a single quintuplet of Bi 2 Se 3 , Bi 2 Te 2 Se, or Sb 2 Te 3 . We quantify the proximity SOC by fitting a symmetry-derived low energy graphene Hamiltonian to the DFT simulated band structure. The overall results are similar for all different topological insulators; we find a strongly enhanced SOC in graphene, which is of the valley-Zeeman type. The effective model and fitted parameters provide realistic foundations for phenomenological modeling of especially spin transport, and for interpreting future experiments on such structures.
From the detailed analysis of the graphene/Sb 2 Te 3 case, we find also a strongly gate tunable proximity SOC and doping level. We show that by tuning the gate field the graphene Dirac point can be well isolated from the valence band of the topological insulator, and the spin-orbit parameters can change sign as a function of the electric field. Remarkably, for all the investigated electric fields the intrinsic SOC induced in graphene remains of the valley Zeeman type, although the corresponding parameters change sign (simultaneously) at around the fields of about −2 V/nm. For this particular field value the Rashba coupling is predicted to dominate the spin properties. Our results regarding the electric field tunability of the proximity SOC strength is important to interpret recent gate-tunable spin-charge conversion experiments.
As outlook, it will be important to make a systematic investigation of twisted bilayers of graphene and topological insulator quintuplets, to demonstrate further tunability of the proximity induced phenomena in the two important materials.

Table of
ContentsProximity effects are a vital route to modify the electronic states of neighboring materials. Three-dimensional topological insulators, such as Bi 2 Te 2 Se, host topological surface states, due to their strong spin-orbit coupling. When graphene is placed on the surface of a topological insulator, giant spin-orbit coupling is induced by the proximity effect, enabling interesting novel electronic properties of its Dirac electrons.

Figure 9 Table of Contents Graphic
Copyright line will be provided by the publisher