Spin Proximity Effects in Graphene/Topological Insulator Heterostructures

Enhancing the spin-orbit interaction in graphene, via proximity effects with topological insulators, could create a novel 2D system that combines nontrivial spin textures with high electron mobility. In order to engineer practical spintronics applications with such graphene/topological insulator (Gr/TI) heterostructures, an understanding of the hybrid spin-dependent properties is essential. {However to date, despite the large number of experimental studies on Gr/TI heterostructures reporting a great variety of remarkable (spin) transport phenomena, little is known about the true nature of the spin texture of the interface states as well as their role on the measured properties. Here we use {\it ab initio} simulations and tight-binding models to determine the precise spin texture of electronic states in graphene interfaced with a Bi$_2$Se$_3$ topological insulator. Our calculations predict the emergence of a giant spin lifetime anisotropy in the graphene layer, which should be a measurable hallmark of spin transport in Gr/TI heterostructures, and suggest novel types of spin devices


Introduction
Following the discovery of graphene and the large family of van der Waals heterostructures based on two-dimensional materials, [1][2][3][4][5] the spectrum of practical applications harnessing the uniqueness of such materials has grown continuously. 6 In the area of spintronics, 7,8 the long room-temperature spin lifetime in graphene opens the possibility of large-scale integration of lateral spintronic devices and architectures. [9][10][11] Additionally, many recent reports indicate the benefit of using proximity effects to tune the spin properties inside the graphene layer and to engineer devices such as spin fieldeffect transistors. 12,13 This provides exciting opportunities in the search for innovative spin manipulation strategies and the development of non-charge-based information processing technologies. Proximity effects have been studied by combining graphene with magnetic insulators, [14][15][16][17] or by magnifying the spin-orbit coupling (SOC) in graphene through extrinsic chemical functionalization. [18][19][20] Another route recently proposed is to interface graphene with transition metal dichalcogenides (TMDCs) 21,22 such as WS 2 or WSe 2 , which leads to phenomena such as weak antilocalization (WAL) [23][24][25] or large nonlocal Hall signals. 26 Additionally, the fact that electronic states in Gr/TMDC sys-tems are spin-polarized primarily along the outof-plane direction 27 results in large spin lifetime anisotropy between in-plane and out-ofplane spin-polarized electrons, which is however weakly energy dependent and therefore not tunable. [28][29][30] Recently, a lot of attention has been paid to heterostructures of graphene and topological insulators (TIs), with reports of anomalous magnetotransport, giant Edelstein effect, and gate-tunable tunneling resistance, [31][32][33][34][35] as well as the possible existence of a quantum spin Hall phase. 36,37 On the more applied side, the fabrication of broadband photodetectors based on Gr/TI heterostructures has been realized, 38 as well as the injection of spin-polarized current from an ultrathin Bi 2 Te 2 Se nanoplatelet into graphene. 39 TI materials are distinguished by their strong intrinsic SOC, which leads to the formation of a bulk band gap and 2D surface states that host massless Dirac fermions with spin-momentum locking. [40][41][42][43][44][45] Proximity to a TI leads to a band gap opening and spinsplit bands in graphene, as discussed theoretically for the case of Bi 2 Se 3 , 46,47 or for graphene interfaced with Sb 2 Te 3 . 37 However, currently there is substantial variability in the literature concerning the precise spin characteristics of Gr/TI systems. Rajput and coworkers measured and calculated a spin splitting of ∼80 meV in graphene on Bi 2 Se 3 , 47 while Lee et al. calculated a band gap of up to 20 meV induced in graphene by Bi 2 Te 2 Se when all Dirac cones coincided. 48 Kou et al. predicted a SOC of ∼2 meV induced in graphene sandwiched between two layers of Sb 2 Te 3 . 49 They also pointed out, as did Lin et al., 50 the importance of the Kekulé distortion on the magnitude of the band gap in graphene. Jin and Jhi reported a TI thickness dependence of SOC induced in graphene by Sb 2 Te 3 , and they also hinted at unusual spin textures induced in the graphene bands. 51 Meanwhile, De Beule et al. 52 concluded that the spin texture imprinted on the graphene states should resemble the standard Rashba texture, as also found in Zhang et al. 53 Overall, these works indicate that TIs clearly induce strong proximity effects, resulting in gap opening and spin splitting of the bands, but the precise nature of the spin texture induced in the graphene layer, and the way to detect it experimentally are crucially lacking. Additionally, such information is not only essential for clarifying how proximity effects between graphene and TIs generate the measured properties, but could also enlarge the possibilities for tailoring spintronics applications.
In this article, we report fundamental spin transport properties of Gr/TI heterostructures, by performing ab initio calculations and fitting to tight-binding (TB) models that fully reproduce both the band structure and the spin texture in the graphene layer. Structures with different twist angles between the graphene and the TI are considered, but in all cases a giant spin lifetime anisotropy in the graphene layer, with in-plane spins relaxing much faster than out-of-plane spins is obtained. In the highly commensurate structure, with a twist angle of 30 • , the anisotropy is maximal near the graphene Dirac point, reaching values of tens to hundreds, and decays to 1/2 at higher energies. Meanwhile, in the larger unit cell, with a twist angle of 0 • , the anisotropy remains high at all Fermi energies and exhibits a strong electronhole asymmetry. The difference in these behaviors is driven by the dominating SOC terms in each structure, which depend on the specific interface symmetry. This contrasts with prior calculations of the spin texture in Gr/TI systems, 52,53 which predicted a purely Rashbalike spin texture with an energy-independent anisotropy of 1/2.
Our theoretical predictions could be experimentally confirmed by performing a gatedependent measurement of the anisotropy as recently achieved in Gr/TMDC samples; 29,30,54,55 heterostructures which however do not exhibit any gate-dependence. Differently, Gr/TI allows for a strong gate-dependent anisotropy enabling the fabrication of tunable spin filtering devices, while the in-plane spin-momentum locking could also make it possible to convert charge current to spin current and to control the spin orientation of the current. 39

Model and Methods
Ab initio calculations of the electronic structure of Gr/TI heterostructures were carried out using density functional theory (DFT), 56 implemented in the Vienna Ab initio Simulation Package (VASP), 57 with the wave functions expanded in a plane wave basis with an energy cutoff of 600 eV, using the projector augmented wave method. 58 The PBE form of the generalized gradient approximation 59 was used to compute the exchange-correlation energy, and a 24 × 24 × 1 (9 × 9 × 1) k-point mesh for the small (large) unit cell was used together with a convergence criterion of 10 −6 eV. The spin-orbit coupling was included through noncollinear calculations, while the Van der Waals force was accounted for based on the Tkatchenko-Scheffler method, 60 and all structures were fully relaxed until forces were smaller than 10 −2 eV/Å. Figure 1 shows the simulated Gr/TI heterostructures, which contain a Bi 2 Se 3 film and one monolayer of graphene in either a √ 3 × √ 3 or a 5 × 5 supercell, shown in Figs. 1(a) and (b), respectively. Since the graphene layer is attached to the TI substrate, the minimumenergy lattice constant of bulk Bi 2 Se 3 , 4.196 A, was adopted. The relaxed crystal structures exhibit a lattice mismatch of less than 3%, while the relaxed interlayer spacing between the graphene and TI is larger than 3.5Å. For the small unit cell we considered TI films of two different thicknesses, one and six quintuple layers (QLs), as indicated in Fig. 1(c), while for the large unit cell we only considered a thickness of 1QL. Since our results are qualitatively independent of the number of QLs, we restrict our discussion to the 1QL case, and show the 6QL results in the Supplemental Material. Figure  1(a) depicts the hollow configuration, with the top (green) Se atom in the center of a carbon ring. This is the most stable configuration, but when studying the spin texture in this unit cell we have also considered other alignments. To describe the electronic properties of graphene on a TI, we employ a TB Hamil- where c † is (c is ) is the creation (annihilation) operator of an electron at lattice site i with spin s, d ij (D ij ) is the unit vector pointing from site j to nearest (next-nearest) site i, s are the spin Pauli matrices, ν ij = +1(−1) for a clockwise (counterclockwise) hopping path from site j to i, ξ i = +1(−1) on sublattice A (B), and the single (double) brackets are sums over first (second) nearest neighbors. The first term in Eq. (1) describes the hopping between nearestneighbor carbon atoms. As depicted in Fig  1(a), this has two different strengths: t 0 is the hopping within the carbon ring surrounding the top (green) Se atom, and t p is the hopping between carbon rings. This describes a Kekulé distortion of the graphene lattice in the hollow configuration, and opens a band gap of 2|t 0 − t p | in the absence of SOC. In the large unit cell, Fig. 1(c), there is no Kekulé distortion and thus t 0 = t p . The second term describes intrinsic SOC in the graphene lattice, 61 λ I , and is assumed nonzero only for the carbon ring surrounding the top Se atom; this is highlighted by the solid triangle in Fig. 1(a). In the larger unit cell, this term is uniform. The third term is valley-Zeeman SOC, λ V Z , which couples spin and valley and arises when sublattice symmetry is broken in the graphene layer. 62 For this reason, it is only present in the larger unit cell. The fourth term is a uniform Rashba SOC, λ z R , induced by an electric field perpendicular to the graphene plane. The fifth term is a second Rashba SOC, λ ρ R , arising from a radial in-plane electric field. We have found the best qualitative fit to DFT by choosing a nonuniform in-plane field, with nonzero values of λ ρ R only along the green arrows in Fig. 1(a). In the larger unit cell this term does not exist, owing to the lack of radial symmetry. Finally, the last term is denoted PIA (pseudospin in-version asymmetry) SOC, λ P IA , which is akin to a second-order Rashba SOC and leads to a k-linear spin splitting of the bands. This particular term only arises in the presence of sublattice symmetry breaking plus a perpendicular electric field. 62

Results
We first focus on the highly commensurate structure, since this is the structure that has been exclusively studied in the literature up to now. Figure 2(a) shows the DFT band structure of this system, where the red symbols denote the projection onto the carbon atoms. Since the Bi 2 Se 3 layer is only 1QL thick, the surface states do not form, but they do appear in the 6QL structure (see Supplemental Fig. S1). In the Gr/Bi 2 Se 3 heterostructure, there is significant charge transfer between the graphene and the TI, which induces a strong p-doping of the graphene and pushes its Dirac point into the TI conduction band. An equivalent amount of charge transfer is also seen in the bridge and top configurations of this unit cell. Because of the √ 3 × √ 3 supercell, the graphene Dirac cones are folded from the K and K points onto the Γ point of the first Brillouin zone, resulting in eight nearly-degenerate bands. This degeneracy is broken by the Kekulé distortion and the SOC induced by the TI, resulting in a band gap opening and spin splitting of the conduction and valence bands. This can be seen in Fig. 2(b), which shows a closeup of the graphene bands near the Dirac point. Black symbols are the DFT results and the red lines are the fits using the TB model of Eq. (1). A band gap and spin splitting on the order of a few meV are observed. As mentioned above, each of these bands are nearly doubly degenerate due to the folding of K/K to Γ, but this degeneracy is broken by the Kekulé distortion and the SOC (see Supplemental Fig. S2). Figure 3 shows an overview of the spin texture of the graphene bands in the Gr/Bi 2 Se 3 heterostructure. The top, middle, and bottom rows show the projections of the spin of the highest conduction band along the x, y, Black, red, and green curves are energy contours corresponding to 34, 100, and 149 meV above the graphene Dirac point. and z axes respectively, plotted as a function of the angle θ around a constant energy contour. The first column contains the DFT results and the second column contains the TB fit. The black, red, and green curves indicate the energy dependence of the spin texture, at Fermi energies of 34, 100, and 149 meV above the graphene Dirac point. Several characteristic features of the spin texture can be seen. The first is that the x and y components exhibit an overall Rashba-like behavior, with S x ∼ − sin θ and S y ∼ cos θ. However, this overall trend is punctuated by sharp minima every 60 • . Second, the z component of the spin is generally nonzero and shows maxima at these same points. These peaks correspond to points of anticrossing between the K and K bands that were folded to the Γ point, which is enabled by the valley-mixing Kekulé distortion and the in-plane Rashba SOC (see Supplemental Fig. S2). Finally, the energy dependence shows that the in-plane components of the spin become weaker near the Dirac point, while the weight of the out-of-plane component increases. Table 1 shows the TB parameters that best reproduce the DFT band structure and spin texture of the Gr/TI heterostructures, with the left column showing the case for the smaller unit cell. The orbital gap induced by the Kekulé distortion is ∼6 meV, and the SOC strengths are on the order of a few meV. A notable result is the relative magnitude of the in-plane and out-of-plane Rashba terms. Recent work has found good fits to the DFT band structure when assuming λ ρ R λ z R , 51 but here to obtain the proper in-plane spin texture it is necessary to enforce λ ρ R λ z R . An overview of the dependence of the in-plane spin texture on λ z R and λ ρ R can be found in Fig. S3 of the Supplemental Material. The intrinsic SOC λ I is necessary for a proper fit to the DFT band structure, but has no impact on the spin texture. Meanwhile, the out-of-plane spin component S z depends crucially on the presence of the Kekulé distortion, which hybridizes the K and K bands; in its absence the magnitude of S z drops by three orders of magnitude (see Supplemental Fig. S2). As stated above, the valley-Zeeman and PIA terms do not appear in this system because sublattice symmetry is not broken. Here we note that this spin texture is quite different from that in prior works, which predicted a purely Rashba-like behavior. 52, 53 We attribute this difference to the choice of model used; the earlier works used a continuum model for the graphene and TI bands that does not account for trigonal warping, the Kekulé distortion, or the in-plane Rashba terms. We have found that both of these terms are crucial for properly capturing the DFT spin texture. Figure 4 shows the band structure and spin texture of the larger unit cell (see Fig. 1(c)), and the right column of Table 1 shows the TB fitting parameters. It is clear that, owing to the different interface symmetry of the larger unit cell, there are significant differences in the band structure, spin texture, and relevant SOC parameters compared to the smaller unit cell. In the band structure, the graphene Dirac cones remain separated at the K and K points of the Brillouin zone while the charge transfer between the graphene and the TI remains large (see Supplemental Fig. S6). Owing to this lack of band folding, hybridization between the valleys no longer occurs and the 60 • periodicity of the spin texture disappears. Instead, S z remains independent of the momentum direction, and its sign is valley-dependent. This behavior is driven by the presence of valley-Zeeman SOC, λ V Z , which is permitted by the sublattice symmetry breaking in the larger unit cell. The in-plane spin components follow the typical Rashba texture, and the PIA SOC determines their energy dependence. As mentioned above, the Kekulé distortion and in-plane Rashba SOC are not present in this system. Additionally, the intrinsic SOC is found to be vanishingly small. It it interesting to note that this TB model is identical to that used for Gr/TMDC systems, and the obtained fitting values are also quite similar. 27 A comparison of Figs. 3 and 4 indicates that the symmetry of the system can have a significant impact on the spin texture induced in graphene by proximity to a TI substrate. This then raises the question: which of these scenarios is most likely to be encountered in an experimental setup? For experiments that interface graphene with TIs via stacking of exfoliated layers, the situation seen in the larger unit cell seems more likely, as the alignment and relative orientation of the graphene and TI lattices remains largely uncontrolled. 63,64 However, through careful processing and device fabrication, more precise control over the interface may be achieved. 4,65 From the perspective of spin transport, a measurement of the spin lifetime anisotropy can be invaluable in determining the dominant SOC terms and the nature of spin relaxation in these systems. 54,55 Indeed, predictions of fast in-plane spin relaxation in Gr/TMDC heterostructures, driven by spinvalley locking and intervalley scattering, 28 have recently been confirmed by measurements of the spin lifetime anisotropy. 29,30 From our calculated spin textures we can predict the spin lifetime anisotropy of each system, defined as the ratio of out-of-plane to in-plane spin lifetime, ζ ≡ τ s,z /τ s,x . Assuming the D'yakonov-Perel' regime of spin relaxation, the lifetime of spins polarized along α is given by where Ω is the momentum-dependent effective magnetic field arising from SOC in units of spin precession frequency, τ * β is the time to randomize the βcomponent of Ω, with β ⊥ α, and the overline represents an average over the Fermi surface at a particular Fermi energy. 66 For a given energy band, the effective magnetic field can be decomposed as Ω = ωS, where ω = ∆E/ is the spin precession frequency associated with the spin splitting ∆E of the band, and S = ψ|s|ψ is the spin polarization of the eigenstates ψ associated with the band. The spin lifetime anisotropy arising from the spin-split band structure can then be written as where the sum over i includes each of the four conduction or valence bands in the Fermi surface average. In the small unit cell, because both Dirac cones are folded to the Γ-point we have τ * x = τ * z = τ p , with τ p the momentum relaxation time. In the large unit cell the graphene Dirac cones remain at K and K , and owing to the presence of λ V Z we have τ * x = τ p and τ * z = τ iv , where τ iv is the intervalley scattering time. 28 Figure 5 shows the spin lifetime anisotropy in the Gr/TI heterostructures, calculated from Eq. (2). Panels (a) and (b) are for the small and large unit cell respectively, the open circles show the DFT results, and the solid lines are from the TB fits. In the small unit cell, the anisotropy remains negligible away from the graphene Dirac point, on the order of 1/2, and reaches values in the hundreds at the lowest energies. This trend also holds within the TI bandgap, where the graphene bands are completely in-plane and the anisotropy is 1/2. Such behavior results from the increase in weight of S z near the Dirac point and the corresponding decrease of S x and S y . The anisotropy obtained from the TB fit is both qualitatively and quantitatively similar to the DFT results. It should be noted, however, that the TB model does not fully account for all aspects of the behavior of S z . In particular, the DFT results show that the lowest pair of conduction and valence bands exhibit no out-of-plane spin texture, while the upper pair of conduction and valence bands show S z similar to Fig. 3(c) (see Supplemental Fig. S4). Meanwhile, the TB results show identical magnitude of S z for all bands; this leads to an overstimation of ζ by approximately a factor of two. Figure 5(b) shows the anisotropy in the large unit cell, in units of τ iv /τ p . Both the TB fit and the analytical prediction derived for Gr/TMDC systems 28 show nice agreement with the DFT results. In this case, the anisotropy is characterized by a strong electron-hole asymmetry, which is driven by the relatively large value of λ P IA ; as shown in Ref. 28, the out-of-plane spin relaxation rate is proportional to (ak F λ P IA ± λ R ) 2 , where a is the graphene lattice constant, k F is the Fermi wave number, and the +(-) is for the conduction (valence) band. At sufficiently negative energies, when ak F λ P IA = λ R , this model predicts that the spin lifetime anisotropy will diverge. In reality, when τ ⊥ s becomes sufficiently long another source of spin relaxation, such as contact dephasing or magnetic impurities, would take over, placing an upper bound on ζ. In systems without a strong PIA SOC, the anisotropy would be independent of the Fermi energy. In the large unit cell, the anisotropy is driven by the SOC and the charge scattering through τ iv /τ p . In general, intervalley scattering is caused by structural defects such as dislocations, grain boundaries, vacancies, etc., as well as chemical adsorbates such as hydrogen, oxygen, or other hydrocarbons that could be deposited during device fabrication. 67 Bi 2 Se 3 is known to suffer from Se vacancies, which might also induce short-range Coulomb potentials and intervalley scattering in graphene. 67 Measuring τ p is straightforward, as it can be deduced from the mobility and charge density. For example, a typical carrier density of 2 × 10 12 cm −2 coupled with a mobility of 6000 cm 2 /V.s, as measured recently for a graphene/Bi 2 Se 3 system, 33 yields τ p ≈ 100 fs. Determining τ iv requires a measurement of weak localization (WL), but in Gr/TI or Gr/TMDC systems the strong SOC leads to weak antilocalization (WAL), making it difficult to extract τ iv . So far, the best that has been done for a Gr/TMDC system is to measure WL in a region of the device that is not covered by the TMDC, and to assume that value as an upper bound of τ iv in the Gr/TMDC region. 25 We are not aware of any estimates of τ iv in Gr/TI systems. Measurements of WL in graphene systems yield τ iv /τ p in the range of 3-20 depending on the sample quality and Fermi energy. 25,68,69 Very generally, τ p is in the range of tens of fs and τ iv is on the order of hundreds of fs to a few ps. Assuming τ iv ≈ 10τ p as a typical experimental situation, we would expect an anisotropy on the order of a few tens over the full range of gate voltage.
As mentioned in the Model and Methods section, we also used DFT to calculate the spin texture of the Gr/TI heterostructures in the top and bridge configurations, where the top Se atom sits below a carbon atom or a carboncarbon bond, respectively. In contrast to the hollow configuration, these configurations show a small anisotropy, ζ 1 (see Supplemental  Fig. S5). This arises from the lack of Kekulé distortion and/or in-plane Rashba SOC in these other lattice arrangements, which as discussed above, are both necessary to enhance the magnitude of S z and thus the anisotropy.

Conclusions
In summary, our study reveals the emergence of anisotropic spin transport in graphene in proximity with topological insulators, but the origin and energy dependence of this anisotropy vary significantly with the geometry of the interface. This arises from the very small lattice mismatch, which permits a highly commensurate unit cell at the appropriate twist angle. This is in contrast to the case of Gr/TMDC systems, which have a much larger lattice mismatch that precludes the formation of a small and highly commensurate unit cell. 23,27 Similarly to the case of Gr/TMDC, Gr/TI displays an almost energy-independent anisotropy for zero twist angle between the graphene and TI lattices. However, in the highly commensurate unit cell (with a twist angle of 30 • ), the spin anisotropy is connected to both a Kekulé distortion and an in-plane Rashba SOC induced in the graphene by the TI. As a result, the spin lifetime becomes highly anisotropic near the graphene Dirac point while vanishing at higher energies, suggesting a much stronger variability via electrostatic gating in experiments. Such spin anisotropy could be playing a role in the debated experimental results reported to date in Gr/TI heterostructures, 32,33,39 while simultanously suggesting new device engineering such as gate-tunable linear spin polarizers, which remove the in-plane component of a spin-polarized current but leave the out-ofplane component intact.
One useful observation is that, as shown in Fig. 2, the Fermi level initially lies in the Bi 2 Se 3 conduction band, which will generate parallel transport in the graphene and TI layers. However, given that the spin lifetime in the TI bulk should be exceptionally short (few femtoseconds), 70,71 any measured spin signal may still carry features of the spin transport in the graphene layer. To more optimally realize the conditions in which the TI surface states would play a role in the transport properties of Gr/TI heterostructures, ternary compounds, with the TI Fermi energy well within the TI bulk gap, 72 would be even more desirable.
Finally, it would be interesting to include defects and disorder in the ab initio simulations and TB models, since this could locally alter the strength and nature of the SOC parameters. Such analysis, beyond the scope of the present work, could be also extended by developing a full Gr/TI tight-binding model, using for instance the Fu-Kane-Mele model. 73,74 Acknowledgement ICN2 is funded by the CERCA Programme / Generalitat de Catalunya, and is supported by the Severo Ochoa program from Spanish MINECO (Grant No. SEV-2013-0295). The authors acknowledge funding from the Spanish Ministry of Economy and Competitiveness and the European Regional Development Fund (Project No.