Single and bilayer graphene on the topological insulator Bi$_2$Se$_3$: Electronic and spin-orbit properties from first principles

TL;DR

This study uses first-principles calculations to analyze how graphene's electronic and spin-orbit properties are affected when placed on top of the topological insulator Bi$_2$Se$_3$, revealing tunable SOC effects and potential for spintronic applications.

Contribution

It provides a detailed first-principles analysis of proximity-induced spin-orbit coupling in graphene on Bi$_2$Se$_3$, including effects of electric fields and interlayer distance, which was not previously characterized.

Findings

Graphene becomes hole doped by 350 meV while maintaining linear dispersion.

Proximity-induced SOC in graphene is about 1 meV, valley-Zeeman type, weakly dependent on Bi$_2$Se$_3$ layers.

Electric fields can tune band offsets and SOC, enabling a spin-orbit valve in bilayer graphene.

Abstract

We present a detailed study of the electronic and spin-orbit properties of single and bilayer graphene in proximity to the topological insulator Bi$_2$Se$_3$. Our approach is based on first-principles calculations, combined with symmetry derived model Hamiltonians that capture the low-energy band properties. We consider single and bilayer graphene on 1--3 quintuple layers of Bi$_2$Se$_3$ and extract orbital and proximity induced spin-orbit coupling (SOC) parameters. We find that graphene gets significantly hole doped (350 meV), but the linear dispersion is preserved. The proximity induced SOC parameters are about 1 meV in magnitude, and are of valley-Zeeman type. The induced SOC depends weakly on the number of quintuple layers of Bi$_2$Se$_3$. We also study the effect of a transverse electric field, that is applied across heterostructures of single and bilayer graphene above 1 quintuple layer of Bi$_2$Se$_3$. Our results show that band offsets, as well as proximity induced SOC parameters can be tuned by the field. Most interesting is the case of bilayer graphene, in which the band gap, originating from the intrinsic dipole of the heterostructure, can be closed and reopened again, with inverted band character. The switching of the strong proximity SOC from the conduction to the valence band realizes a spin-orbit valve. Additonally, we find a giant increase of the proximity induced SOC of about 200%, when we decrease the interlayer distance between graphene and Bi$_2$Se$_3$ by only 10%. Finally, for a different substrate material Bi$_2$Te$_2$Se, band offsets are significantly different, with the graphene Dirac point located at the Fermi level, while the induced SOC strength stays similar in magnitude compared to the Bi$_2$Se$_3$ substrate.

Tables

Table 1. Table 1: Fit parameters of Hamiltonian ℋ SLG subscript ℋ SLG \mathcal{H}_{\textrm{SLG}} for the SLG/Bi 2 Se 3 stacks for different number of QLs. The Fermi velocity v F subscript 𝑣 F v_{\textrm{F}} , gap parameter Δ Δ \Delta , Rashba SOC parameter λ R subscript 𝜆 R \lambda_{\textrm{R}} , intrinsic SOC parameters λ I A superscript subscript 𝜆 I A \lambda_{\textrm{I}}^{\textrm{A}} and λ I B superscript subscript 𝜆 I B \lambda_{\textrm{I}}^{\textrm{B}} , and PIA SOC parameters λ PIA A superscript subscript 𝜆 PIA A \lambda_{\textrm{PIA}}^{\textrm{A}} and λ PIA B superscript subscript 𝜆 PIA B \lambda_{\textrm{PIA}}^{\textrm{B}} . The Dirac point energy E D subscript 𝐸 D E_{\textrm{D}} , as defined in Fig. 2 (a).

QLs	$v_{\text{F}}/{10}^5\left[\frac{\text{m}}{\text{s}}\right]$	$\Delta$ [$\mu$eV]	${\lambda }_{\text{R}}$ [meV]	${\lambda }_{\text{I}}^{\text{A}}$ [meV]	${\lambda }_{\text{I}}^{\text{B}}$ [meV]	${\lambda }_{\text{PIA}}^{\text{A}}$ [meV]	${\lambda }_{\text{PIA}}^{\text{B}}$ [meV]	$E_{\text{D}}$ [meV]
1	8.134	0.6	-0.771	1.142	-1.135	0.465	0.565	353.2
2	8.131	1.1	-0.691	1.221	-1.211	1.834	1.733	352.1
3	8.126	0.4	-0.827	1.343	-1.330	2.904	2.976	509.1
1111with additional hBN above SLG, parameter $\Delta$ in meV	8.110	21.32	-0.799	1.638	-1.517	2.708	2.433	341.2

Table 2. Table 2: Fit parameters of Hamiltonian ℋ BLG subscript ℋ BLG \mathcal{H}_{\textrm{BLG}} for the BLG/Bi 2 Se 3 stacks for different number of QLs. Intra- and interlayer hoppings γ 𝛾 \gamma of the BLG, when the lower (upper) GRP layer is placed in potential V 𝑉 V ( − V 𝑉 -V ). The asymmetry in the energy shift of the bonding and antibonding states Δ Δ \Delta . Intrinsic SOC parameters of the corresponding layer and sublattice λ I subscript 𝜆 I \lambda_{\textrm{I}} . Global and local space inversion symmetry breaking SOC parameters λ 0 subscript 𝜆 0 \lambda_{0} and λ R subscript 𝜆 R \lambda_{\textrm{R}} .

QLs	${\gamma }_0$	${\gamma }_1$	${\gamma }_3$	${\gamma }_4$	$V$	$\Delta$	${\lambda }_{\text{I}}^{\text{A1}}$	${\lambda }_{\text{I}}^{\text{B1}}$	${\lambda }_{\text{I}}^{\text{A2}}$	${\lambda }_{\text{I}}^{\text{B2}}$	${\lambda }_0$	${\lambda }_{\text{R}}$	$E_{\text{D}}$
	[eV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]	[meV]
1	2.513	373.5	-274.9	-165.3	41.60	13.79	-1.056	1.170	0.012	0.012	-0.433	-0.273	183.3
2	2.512	373.5	-275.5	-165.6	43.12	13.79	-1.117	1.301	0.012	0.012	-0.167	-0.199	181.4
3	2.513	373.9	-264.9	-163.5	23.63	13.54	-1.030	1.168	0.012	0.012	-0.190	-0.188	19.95

Table 3. Table 3: Fit parameters of Hamiltonian ℋ SLG subscript ℋ SLG \mathcal{H}_{\textrm{SLG}} for the SLG/Bi 2 Te 2 Se stacks for different number of QLs. The Fermi velocity v F subscript 𝑣 F v_{\textrm{F}} , gap parameter Δ Δ \Delta , Rashba SOC parameter λ R subscript 𝜆 R \lambda_{\textrm{R}} , intrinsic SOC parameters λ I A superscript subscript 𝜆 I A \lambda_{\textrm{I}}^{\textrm{A}} and λ I B superscript subscript 𝜆 I B \lambda_{\textrm{I}}^{\textrm{B}} , and PIA SOC parameters λ PIA A superscript subscript 𝜆 PIA A \lambda_{\textrm{PIA}}^{\textrm{A}} and λ PIA B superscript subscript 𝜆 PIA B \lambda_{\textrm{PIA}}^{\textrm{B}} . The Dirac point energy E D subscript 𝐸 D E_{\textrm{D}} , as defined in Fig. 13 (a).

QLs	$v_{\text{F}}/{10}^5\left[\frac{\text{m}}{\text{s}}\right]$	$\Delta$ [$\mu$eV]	${\lambda }_{\text{R}}$ [meV]	${\lambda }_{\text{I}}^{\text{A}}$ [meV]	${\lambda }_{\text{I}}^{\text{B}}$ [meV]	${\lambda }_{\text{PIA}}^{\text{A}}$ [meV]	${\lambda }_{\text{PIA}}^{\text{B}}$ [meV]	$E_{\text{D}}$ [meV]
1	8.123	0.3	-0.669	1.353	-1.351	-1.091	-1.209	4.0
2	8.105	0.4	-0.487	1.446	-1.441	1.317	1.350	2.4
3	8.109	1.8	-0.521	1.466	-1.460	1.581	1.399	-0.7

Equations

$${\mathcal{H}}_{\text{SLG}}={\mathcal{H}}_0+{\mathcal{H}}_{\Delta }+{\mathcal{H}}_{\text{I}}+{\mathcal{H}}_{\text{R}}+{\mathcal{H}}_{\text{PIA}}+E_D,$$

$${\mathcal{H}}_0=\hslash v_{\text{F}}(\tau k_x{\sigma }_x-k_y{\sigma }_y)\otimes s_0,$$

$${\mathcal{H}}_{\Delta }=\Delta {\sigma }_z\otimes s_0,$$

$${\mathcal{H}}_{\text{I}}=\tau ({\lambda }_{\text{I}}^{\text{A}}{\sigma }_{+}+{\lambda }_{\text{I}}^{\text{B}}{\sigma }_{-})\otimes s_z,$$

$${\mathcal{H}}_{\text{R}}=-{\lambda }_{\text{R}}(\tau {\sigma }_x\otimes s_y+{\sigma }_y\otimes s_x),$$

$${\mathcal{H}}_{\text{PIA}}=a({\lambda }_{\text{PIA}}^{\text{A}}{\sigma }_{+}-{\lambda }_{\text{PIA}}^{\text{B}}{\sigma }_{-})\otimes (k_x{s}_y-k_y{s}_x).$$

$${\mathcal{H}}_{\text{BLG}}={\mathcal{H}}_{\text{orb}}+{\mathcal{H}}_{\text{soc}}+E_D,$$

$${\mathcal{H}}_{\text{orb}}=\left(\begin{array}{cccc}\Delta +V& {\gamma }_{0}f(\mathbf{k})& {\gamma }_4{f}^{\ast }(\mathbf{k})& {\gamma }_1\\ {\gamma }_0{f}^{\ast }(\mathbf{k})& V& {\gamma }_{3}f(\mathbf{k})& {\gamma }_4{f}^{\ast }(\mathbf{k})\\ {\gamma }_{4}f(\mathbf{k})& {\gamma }_3{f}^{\ast }(\mathbf{k})& -V& {\gamma }_{0}f(\mathbf{k})\\ {\gamma }_1& {\gamma }_{4}f(\mathbf{k})& {\gamma }_0{f}^{\ast }(\mathbf{k})& \Delta -V\end{array}\right)\otimes s_0,$$

$${\mathcal{H}}_{\text{soc}}=\left(\begin{array}{cccc}\tau {\lambda }_{\text{I}}^{\text{A1}}{s}_z& \text{i}({\lambda }_0+2{\lambda }_{\text{R}})s_{-}^{\tau }& 0& 0\\ -\text{i}({\lambda }_0+2{\lambda }_{\text{R}})s_{+}^{\tau }& -\tau {\lambda }_{\text{I}}^{\text{B1}}{s}_z& 0& 0\\ 0& 0& \tau {\lambda }_{\text{I}}^{\text{A2}}{s}_z& -\text{i}({\lambda }_0-2{\lambda }_{\text{R}})s_{-}^{\tau }\\ 0& 0& \text{i}({\lambda }_0-2{\lambda }_{\text{R}})s_{+}^{\tau }& -\tau {\lambda }_{\text{I}}^{\text{B2}}{s}_z\end{array}\right).$$

Full text

spintronics, graphene, heterostructures, proximity spin-orbit coupling, topological insulator

I Introduction

Van der Waals (vdW) heterostructures Geim and Grigorieva (2013); Novoselov et al. (2016); Duong et al. (2017) and emerging proximity effects Žutić et al. (2019) are an ideal platform to induce tailored properties in two-dimensional (2D) materials. Prominent 2D material examples are semimetallic single layer graphene Castro Neto et al. (2009) (SLG), semiconducting transition-metal dichalcogenides Wang et al. (2012) (TMDCs), and insulating hexagonal boron-nitride Catellani et al. (1987) (hBN). Recently, also superconductors Frindt (1972) (NbSe2), and ferro- and antiferromagnets Dillon and Olson (1965); Huang et al. (2017); McGuire et al. (2015); Wiedenmann et al. (1981); Carteaux et al. (1995); Gong et al. (2017) (CrI3, Cr2Ge2Te6, MnPSe3) have been added to the list of 2D materials. Within this ever-expanding field of vdW structures, there already are subfields, such as valleytronics Langer et al. (2018); Zhong et al. (2017); Bussolotti et al. (2018); Schaibley et al. (2016), straintronics Lin et al. (2018); Roldán et al. (2015); Fang et al. (2018), twistronics David et al. (2019); Li and Koshino (2019), and spintronics Žutić et al. (2004); Han et al. (2014); Fabian et al. (2007), wherein several major achievements have been made, for example optical spin injection (Avsar et al., 2017) in SLG or tunable valley polarization in a TMDC Seyler et al. (2018), which are only possible due to vdW heterostructures and proximity effects.

Another important large class of materials are the three-dimensional (3D) topological insulators such as Zhang et al. (2009a) Bi2Se3, Bi2Te3, and Sb2Te3, which are also layered crystals, consisting of quintuple layers (QLs) of alternating chalcogen (Se, Te) and pnictogen (Bi, Sb) atoms, which are held together by vdW forces. However, the characteristic Dirac states with spin-momentum locking Hsieh et al. (2009) emerge only when the top and bottom surfaces of the topological insulator decouple, occurring at already 5–6 QLs, as demonstrated by angle resolved photoemission spectroscopy Zhang et al. (2010) and first-principles calculations Liu et al. (2010); Yazyev et al. (2010); Park et al. (2010). Since each QL is about 1 nm in thickness, these materials are in between the 2D and 3D regime, depending on how many QLs one investigates. Nevertheless, they are important for practical applications Tian et al. (2017), due to their topologically protected Hasan and Kane (2010); Zhang et al. (2009b) and well conducting surface states Koirala et al. (2015), and for proximity induced phenomena Song et al. (2018); Khokhriakov et al. (2018); Jafarpisheh et al. (2018); Zhang et al. (2014), since strong SOC is present. When the topological insulators act as a substrate, the 2D regime (1–2 QLs) is sufficient, as proximity effects are of short range nature.

Recently, the interface engineering of 2D materials has become an important topic Hesjedal and Chen (2017); Zhao et al. (2017). Experimentalists and theorists are searching for material combinations with novel properties. Graphene, due to its extremely high electron mobility Banszerus et al. (2015) and intrinsically small SOC Gmitra et al. (2009), is perfectly suited for spintronics. In addition, this monolayer carbon sheet can be efficiently manipulated by short range proximity effects.

One can induce strong SOC, as well as magnetism in SLG Gmitra and Fabian (2015); Gmitra et al. (2016); Song et al. (2018); Zollner et al. (2016); Jafarpisheh et al. (2018); Khokhriakov et al. (2018). Similar to a TMDC Gmitra and Fabian (2015); Gmitra et al. (2016), a topological insulator strongly enhances the rather weak intrinsic SOC of SLG from 10 $\mu$eV Gmitra et al. (2009), to about 1–2 meV Song et al. (2018); Jafarpisheh et al. (2018). Phase coherent transport measurements of SLG on Bi1.5Sb0.5Te1.7Se1.3 have shown Dyakonov-Perel type spin relaxation with proximity induced SOC of at least 2.5 meV Jafarpisheh et al. (2018). First-principles calculations of SLG on Bi2Se3 have found either pure intrinsic or valley-Zeeman type SOC in the meV range, depending on the twist angle Song et al. (2018). As a consequence of the large induced SOC in SLG, the spin lifetimes of electrons significantly decrease Cummings et al. (2017); Song et al. (2018), from nanoseconds down to the picosecond range. Giant spin lifetime anisotropies, the ratio of out-of-plane to in-plane spin lifetimes, can be achieved Cummings et al. (2017); Song et al. (2018).

Bilayer graphene (BLG) is even more interesting, since only the layer closest to the proximitizing material gets modified, allowing for highly efficient tuning of the proximity properties by gating and doping. Recent studies have shown short range proximity induced exchange or SOC in BLG on Cr2Ge2Te6 or WSe2 Zollner et al. (2018); Gmitra and Fabian (2017); Khoo et al. (2017). Due to the unique and tunable low energy band structure of BLG, all-electrical control of spin relaxation and polarization can be achieved in such heterostructures. The proposed spin-orbit and exchange valve effects Zollner et al. (2018); Gmitra and Fabian (2017) in proximitized BLG can lead to new opportunities for spintronics devices. So far, there are only few experimental studies Zalic et al. (2017); Song et al. (2010); Dang et al. (2010); Steinberg et al. (2015) of BLG on Bi2Se3. However, there is no detailed experimental and theoretical study of the electronic and spin-orbit properties of BLG interfaced with topological insulators.

The open questions we would like to address are as follows. In the case of SLG on Bi2Se3, when valley-Zeeman SOC is predicted to be present Song et al. (2018) for 1 QL, how does the presence of more QLs influence this interesting result? (The result is interesting, since in the topological substrate the spin-orbit fields are in-plane, while the induced valley-Zeeman fields in graphene are out-of-plane). Another important question is, what is the influence of the topological insulator on BLG? What are the band offsets, doping levels, and orbital and spin-orbit proximity effects? How does the interlayer distance between SLG and Bi2Se3 affect the magnitude of proximity SOC? Also, can an electric field tune SOC in SLG and BLG in proximity to the topological insulator? Strong hole doping of the SLG on Bi2Se3 is predicted Song et al. (2018). Can one find a different topological insulator with a better band alignment, such that the graphene Dirac point is near EF?

In this article we investigate these questions using first-principles calculations of SLG and BLG on the topological insulator Bi2Se3. We study the proximity induced SOC in SLG and BLG, originating from the topological insulator, by varying the number of QLs of Bi2Se3 from 1–3. Symmetry-derived low energy tight-binding model Hamiltonians for SLG and BLG are fitted to the first-principles band structures, to extract orbital and spin-orbit parameters of the proximitized materials. Our results show, that the dispersion of SLG (BLG) is preserved, but strong hole doping appears, as the Dirac point is about 350 meV (200 meV) above the Fermi level. The proximity induced SOC is about 1 meV in magnitude, but with opposite sign for A and B sublattice, the so called valley-Zeeman type. We find that the intrinsic SOC parameters increase by about 10%, for every QL of Bi2Se3 that we add, up to 3QLs. As proximity effects are short ranged, we expect this increase to saturate.

Furthermore, we study the effect of a transverse electric field on the low energy band parameters, for 1QL of Bi2Se3 proximitizing SLG and BLG. The electric field can tune the SOC parameters, which can have significant impact for tuning spin lifetimes in SLG and BLG. In addition, the surface states of the topological insulator, as well as the Dirac points of SLG and BLG, can be tuned with respect to the Fermi level, by the field. Most interesting is the BLG case, in which only the low energy conduction band is strongly spin-orbit split for zero field, as a consequence of short range proximity effects, and atom and layer localized energy states. The tuning of the orbital gap of BLG, by gradually increasing the field, leads to a gap closing and subsequent reopening, now with a strongly spin-orbit split valence band. Consequently, a spin-orbit valve effect can be realized, similar to BLG on a TMDC Gmitra and Fabian (2017). An interlayer distance study, between SLG and the Bi2Se3 substrate, shows a giant increase of the proximity induced SOC of about 200%, when we decrease the interlayer distance by only 10%. Furthermore, an atomically modified substrate Bi2Te2Se, leads to a significantly different band alignment and enhanced proximity SOC in SLG. The Dirac point of SLG, as well as the surface states of the topological insulator, are now located at the Fermi level. Especially this system holds promise, for the simultaneous study of two very different spin-orbit fields: in-plane spin-momentum locking from the topological insulator and out-of-plane proximity induced spin-orbit field from SLG.

II Geometry & Computational Details

For the calculation of SLG and BLG on the topological insulator Bi2Se3, we consider $5\times 5$ supercells of SLG and BLG (in Bernal stacking) on top of $3\times 3$ supercells of Bi2Se3. Initial atomic structures are set up with the Atomic Simulation Environment (ASE) Bahn and Jacobsen (2002). We marginally stretch the lattice constant of graphene Castro Neto et al. (2009) to $a=2.486$ Å and leave the Bi2Se3 lattice constants Nakajima (1963) unchanged with $a=4.143$ Å and $c=28.636$ Å, using the atomic parameters $(u,v)=(0.4008,0.2117)$. We consider only geometries without relaxation, using interlayer distances of $3.5$ Å between the lowest graphene layer and the topmost QL of Bi2Se3, in agreement with recent studies Song et al. (2018), and an interlayer distance of $3.3$ Å for the BLG, in agreement with experiment Baskin and Meyer (1955). In Fig. 1 we show the geometry of SLG or BLG on top of 1–3 QLs of Bi2Se3, visualized with VESTA Momma and Izumi (2011). Compared to Ref. Song et al. (2018), we only study what they call the ’large unit cell’, where no twist angle between the materials is present.

The electronic structure calculations are performed by density functional theory (DFT) Hohenberg and Kohn (1964) with Quantum ESPRESSO Giannozzi and et al. (2009). Self-consistent calculations are performed with the $k$-point sampling of $9\times 9\times 1$. Only for the largest heterostructures, when 3QLs of Bi2Se3 are considered, a smaller $k$-point sampling of $6\times 6\times 1$ is used, due to computational limitations. We use an energy cutoff for charge density of $500$ Ry, and the kinetic energy cutoff for wavefunctions is $60$ Ry for the relativistic pseudopotentials with the projector augmented wave method Kresse and Joubert (1999) with the Perdew-Burke-Ernzerhof exchange correlation functional Perdew et al. (1996). We also add vdW corrections Grimme (2006); Barone et al. (2009), and Dipole corrections Bengtsson (1999) are included to get correct band offsets and internal electric fields. In order to simulate quasi-2D systems, we add a vacuum of at least $24$ Å, to avoid interactions between periodic images in our slab geometry.

In contrast to the GW method, the generalized-gradient-approximiation (GGA) used here is certainly not the most accurate available choice to describe 3D topological insulators Aguilera et al. (2013); Nechaev et al. (2013); Yazyev et al. (2012). Consequently, the charge transfer between SLG (BLG) and the topological insulator, as well as doping and proximity SOC can be different in more sophisticated calculations. Still, the GGA does show the important band structure features and it allows to make predictions for proximity effects in SLG/Bi2Se3 heterostructures. Furthermore, the systems we consider are too large, making GW calculations computationally inaccessible. In our GGA-based study, we believe that proximity effects can be captured quite nicely on a qualitative and semi-quantitative level. For example, calculations Song et al. (2018), also using GGA have been used to analyze and interpret recent experiments of graphene/topological insulator structures Jafarpisheh et al. (2018).

III Model Hamiltonian

It has been shown that a topological insulator induces strong proximity SOC in SLG Song et al. (2018); Jafarpisheh et al. (2018), on the order of 1 meV. Depending on the exact geometry (twist angle), either pure intrinsic or valley-Zeeman type SOC can be realized Song et al. (2018). The valley-Zeeman SOC has also been observed in SLG/TMDC heterostructures Gmitra and Fabian (2015); Gmitra et al. (2016). For BLG on a TMDC, even a spin-valve effect is proposed to be present Gmitra and Fabian (2017). We want to analyze in detail the influence of the topological insulator Bi2Se3 on the low energy bands of SLG and BLG, in the previously mentioned (non-twisted) supercell configuration, where valley-Zeeman SOC has been found.

III.1 Graphene

The band structure of proximitized SLG can be modeled by symmetry-derived Hamiltonians Kochan et al. (2017); Zollner et al. (2016, 2019); Gmitra and Fabian (2015); Gmitra et al. (2016); Di Sante et al. (2019); Kane and Mele (2005). For our heterostructures, the effective low energy Hamiltonian is

[TABLE]

Here $v_{\textrm{F}}$ is the Fermi velocity and the in-plane wave vector components $k_{x}$ and $k_{y}$ are measured from $\pm$K, corresponding to the valley index $\tau=\pm 1$. The Pauli spin matrices are $s_{i}$, acting on spin space ($\uparrow,\downarrow$), and $\sigma_{i}$ are pseudospin matrices, acting on sublattice space (C${}_{\textrm{A}}$, C${}_{\textrm{B}}$), with $i=\{0,x,y,z\}$. For shorter notation, we introduce $\sigma_{\pm}=\frac{1}{2}(\sigma_{z}\pm\sigma_{0})$. The lattice constant of pristine graphene is $a$ and the staggered potential gap is $\Delta$. The parameters $\lambda_{\textrm{I}}^{\textrm{A}}$ and $\lambda_{\textrm{I}}^{\textrm{B}}$ describe the sublattice resolved intrinsic SOC, $\lambda_{\textrm{R}}$ stands for the Rashba SOC, and $\lambda_{\textrm{PIA}}^{\textrm{A}}$ and $\lambda_{\textrm{PIA}}^{\textrm{B}}$ are the sublattice resolved pseudospin-inversion asymmetry (PIA) SOC parameters. The basis states are $|\Psi_{\textrm{A}},\uparrow\rangle$, $|\Psi_{\textrm{A}},\downarrow\rangle$, $|\Psi_{\textrm{B}},\uparrow\rangle$, and $|\Psi_{\textrm{B}},\downarrow\rangle$, resulting in four eigenvalues $\varepsilon_{1/2}^{\textrm{CB/VB}}$. Note that the model Hamiltonian describes the dispersion relative to the graphene Dirac point. For the first-principles results doping can occur, shifting the Fermi level off the Dirac point. Therefore we introduce another parameter $E_{\textrm{D}}$, which generates a shift of the global model band structure. We call it the Dirac point energy.

III.2 Bilayer graphene

We wish to describe the low energy band structure of the proximitized BLG in the vicinity of the K and K’ valleys. Therefore, we introduce the following Hamiltonian derived from symmetry Konschuh et al. (2012), where we keep only the most relevant terms

[TABLE]

We use the linearized version for the nearest-neighbor structural function $f(\bm{k})=-\frac{\sqrt{3}a}{2}(k_{x}-\textrm{i}k_{y})$, with the graphene lattice constant $a$ and the Cartesian components of the wave vector $k_{x}$ and $k_{y}$ measured from $\pm$K for the valley index $\tau=\pm 1$. Parameters $\gamma$ describe intra- and interlayer hoppings of the BLG, when the lower (upper) graphene layer is placed in potential $V$ ($-V$). The parameter $\Delta$ describes the asymmetry in the energy shift of the bonding and antibonding states. The parameters $\lambda_{\textrm{I}}$ describe the intrinsic SOC of the corresponding layer and sublattice. The combination of parameters $\lambda_{0}$ and $\lambda_{\textrm{R}}$ describe the global and local breaking of space inversion symmetry. For a more detailed description of the parameters, we refer the reader to Ref. Konschuh et al. (2012). The basis is $|\Psi_{\textrm{A1}},\uparrow\rangle$, $|\Psi_{\textrm{A1}},\downarrow\rangle$, $|\Psi_{\textrm{B1}},\uparrow\rangle$, $|\Psi_{\textrm{B1}},\downarrow\rangle$, $|\Psi_{\textrm{A2}},\uparrow\rangle$, $|\Psi_{\textrm{A2}},\downarrow\rangle$, $|\Psi_{\textrm{B2}},\uparrow\rangle$, and $|\Psi_{\textrm{B2}},\downarrow\rangle$. Similar to the SLG case, doping can occur, and again we denote the energy shift by the Dirac point energy $E_{\textrm{D}}$.

IV Band structure, fit results, and spin-orbit fields

Here, we analyze the dependence of the proximity SOC in SLG and BLG on the number of QLs. We show the full calculated band structures, as well as a zoom to the low energy bands originating from SLG or BLG, being proximitized by the topological insulator, and fit the individual model Hamiltonians. The orbital and spin-orbit fit parameters are summarized in a tabular form.

IV.1 Graphene

In Fig. 2 we show the full band structures of SLG above one, two, and three QLs of Bi2Se3. At the K point, one can recognize the Dirac bands originating from bare SLG Gmitra et al. (2009), while at the $\Gamma$ point, the surface states of the topological insulator form. By gradually increasing the number of QLs, the Dirac bands of Bi2Se3 start to form. We do not show a zoom to the surface bands of Bi2Se3, as we are mainly interested in proximity induced SOC in SLG and BLG. However, one would see a pair of linear bands at the $\Gamma$ point, each originating from one surface of the topological insulator, which still hybridize for 3QLs only, exhibiting a gap. At first glance the bands of SLG seem to be not affected by adding more QLs. We find that SLG gets hole doped, as the Dirac point is shifted roughly $350$ meV above the Fermi level. Only for 3QLs, the Dirac point is at about $500$ meV, which happens to appear due to the thicker topological insulator, consistent with recent calculations Song et al. (2018).

In Fig. 3 we show the low energy band properties for SLG/Bi2Se3 for 1QL, with a fit to the model Hamiltonian. We find that the model dispersion, energy splittings, and spin expectation values agree very well with the first-principles data. In Tab. 1 we summarize the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{SLG}}$ for the SLG/Bi2Se3 stacks for different number of QLs. We find that the fit parameters are almost independent on the number of QLs, indicating that only the closest QL is mainly responsible for proximity SOC. However, the intrinsic SOC parameters gradually increase from about $1.1$ meV, for 1QL, to $1.3$ meV, for 3QLs. Such an increasing (and saturating) behavior of the induced SOC, with respect to the number of QLs, has already been reported for other interface configurations, yielding much larger SOC parameters Jin and Jhi (2013). Also the PIA SOC parameters increase with the number of QLs, while the Rashba SOC stays roughly the same. Similar to SLG/TMDC heterostructures Gmitra and Fabian (2015); Gmitra et al. (2016), we find staggered intrinsic SOC, i.e., $\lambda_{\textrm{I}}^{\textrm{A}}\approx-\lambda_{\textrm{I}}^{\textrm{B}}$. In analogy to the SLG/WSe2 heterostructure Gmitra et al. (2016); Frank et al. (2018), we find an inverted band structure, as the gap $\Delta$, being in the $\mu$eV range, is much smaller than the spin splittings of the bands. As a consequence, SLG proximitized by the topological insulator could host protected pseudohelical states Frank et al. (2018). All our results are in agreement with a recent study of SLG on 1QL of Bi2Se3 Song et al. (2018). In particular, the SOC parameters are of the same magnitude (1 meV) and are also of valley-Zeeman type, for their ’large unit cell’ calculation. In addition, we have calculated a heterostructure consisting of hBN/SLG/Bi2Se3 for 1QL, which is especially interesting for comparison to experimental studies where graphene is protected by the insulating hBN and proximitized by the topological insulator. The fit parameters are summarized in Tab. 1. In agreement with Ref. Zollner et al. (2019), the orbital gap parameter $\Delta$ of the proximitized SLG is about 20 meV, much larger than for the case without hBN. Also the intrinsic SOC is a bit larger, while the rest of the fit parameters are barely different.

To further analyze the low energy bands, we have calculated the spin-orbit fields, see Fig. 4, of the four low energy bands $\varepsilon_{1/2}^{\textrm{VB/CB}}$, corresponding to Fig. 3(e). The spin-orbit fields of the two outer (inner) bands rotate clockwise (counter-clockwise), being a clear signature of Rashba SOC. The spin-orbit fields are isotropic around the K point and show no signs of trigonal warping. Only the two inner bands change their $s_{z}$ expectation value, due to the inverted band structure. The induced spin-orbit fields point out-of-plane, while the topological insulator surface states have in-plane spin-orbit fields.

As a comparative case, we also show the calculated spin-orbit fields for graphene above 3QLs, see Fig. 5, where the PIA SOC parameters are significantly larger, see Tab. 1, than in the 1QL case. However, comparing the 1QL and 3QL cases, we do not notice any difference in the spin-orbit fields, as the magnitude (meV) of PIA SOC is still comparable in the two cases, and the other SOC parameters are almost unchanged. Only a much larger PIA SOC can significantly change the low energy band properties Kochan et al. (2017).

IV.2 Bilayer graphene

In Fig. 6 we show the full band structures of BLG above one, two, and three QLs of Bi2Se3. At the K point, one can recognize the parabolic bands originating from bare BLG Konschuh et al. (2012), while at the $\Gamma$ point, the surface states of the topological insulator form, just as for the SLG case. By gradually increasing the number of QLs, the Dirac bands of Bi2Se3 start to form. The BLG bands exhibit a sizable band gap. Additionally, the BLG gets hole doped by about 200 meV, for one and two QLs, whereas the Fermi level is in the band gap of BLG for three QLs. Experiments have reported about 380 meV of hole doping in BLG on Bi2Se3 Zalic et al. (2017), which may be attributed to the additional SiO2 substrate.

In Fig. 7 we show the low energy band properties for BLG/Bi2Se3 for 1QL, with a fit to the model Hamiltonian $\mathcal{H}_{\textrm{BLG}}$. As for the SLG case, we find very good agreement of our model and the first-principles dispersion, energy splittings, and spin expectation values for the proximitized BLG. In contrast to SLG, the BLG exhibits a large band gap of about 80 meV. As a consequence of the Bi2Se3 substrate, and the resulting induced dipole field, the two graphene layers are at a different potential $V$. Since the low energy bands are formed by the non-dimer carbon atoms of BLG, a band gap opens. In this case, the low energy conduction band is formed by the graphene layer closer to the Bi2Se3. The reason is that the conduction band is strongly spin split (2 meV) around the K point, due to proximity induced SOC. Since proximity induced phenomena are short range effects, only the closest graphene layer is affected. The band structure of BLG on Bi2Se3 is very similar to the one of BLG on WSe2 Gmitra and Fabian (2017), where a spin-orbit valve effect has been proposed.

In Tab. 2 we summarize the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{BLG}}$ for the BLG/Bi2Se3 stacks for different number of QLs. We find that the fit parameters are almost independent on the number of QLs, indicating that only the closest QL is mainly responsible for proximity SOC, similar to what we have found from the SLG case. For the fit, we assume for simplicity the pristine graphene intrinsic SOC parameters Gmitra et al. (2009) for the top layer ($\lambda_{\textrm{I}}^{\textrm{A2}}=\lambda_{\textrm{I}}^{\textrm{B2}}=12~{}\mu$eV). In agreement to the monolayer case, we find staggered intrinsic SOC, i.e., $\lambda_{\textrm{I}}^{\textrm{A1}}\approx-\lambda_{\textrm{I}}^{\textrm{B1}}$. Consequently, also proximitized BLG can exhibit topologically protected phases, as recently shown Alsharari et al. (2018). Surprisingly, the gap parameter $V$ diminishes by about 50%, for 3QLs and the Fermi level is now located within the band gap of the BLG.

By looking at the calculated spin-orbit fields, see Fig. 8, we find that the two proximity spin-orbit split conduction bands are strongly $s_{z}$ polarized. In contrast, the two valence bands exhibit a typical Rashba type spin-orbit field. In the case of SLG on TMDCs, the induced strong staggered SOC (opposite sign on the two sublattices) seems to be a consequence of the $s_{z}$ polarized spin-valley locked states from the TMDC, getting imprinted on the SLG (or BLG). In our investigated case, the surface states of the topological insulator have spin-momentum locking, however with in-plane ($s_{x}$ and $s_{y}$) spin components. Still, the induced spin-orbit fields are out-of-plane, similar to the SLG/TMDC case.

So far, there is no complete microscopic understanding why in these two cases (SLG/Bi2Se3 and SLG/TMDC), the induced SOC is of valley-Zeeman type. However, in the case of SLG/TMDC structures, there are already first approaches discussing how the induced SOC depends on the twist angle and the position of the SLG Dirac point with respect to the TMDC band edges David et al. (2019).

V Transverse electric field

We have seen that the parameters and proximity SOC are marginally affected by adding more QLs. We now study the effect of a transverse electric field only for 1QL of Bi2Se3. Can we tune proximity SOC by the electric field? Is a spin-orbit valve effect present in BLG?

V.1 Graphene

In Fig. 9 we show the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{SLG}}$ for the SLG/Bi2Se3 stack as function of a transverse electric field. We find that the dipole of the structure grows linearly with applied electric field. The Fermi velocity $v_{\textrm{F}}$ is only slightly affected by the field, but shows a linear dependence. The gap parameter $\Delta$, reflecting the sublattice symmetry breaking, stays tiny in magnitude without any noticeable trend. The intrinsic as well as Rashba SOC parameters grow in magnitude, when changing the field from negative to positive amplitude. While the intrinsic SOC parameters, $\lambda_{\textrm{I}}^{\textrm{A}}$ and $\lambda_{\textrm{I}}^{\textrm{B}}$, can be changed from about 1 to 1.2 meV, the Rashba SOC parameter $\lambda_{\textrm{R}}$ changes from about 0.6 to 0.9 meV in magnitude, when tuning the electric field from $-2.5$ V/nm to 2.5 V/nm. A roughly linear trend can also be observed in the PIA SOC parameters, which can be tuned from 0 up to 0.8 meV. The Dirac point energy $E_{\textrm{D}}$ stays at the same position, with respect to the Fermi level, as the field changes. The doping energy of the topological insulator $E_{\textrm{TI}}$ decreases, when a negative electric field is applied. Such a field tunability of the SOC parameters, especially the Rashba one, can lead to a giant control of spin-relaxation times and anisotropies in SLG Zollner et al. (2019); Song et al. (2018).

V.2 Bilayer graphene

Most interesting is the case of BLG under the influence of an electric field. As already mentioned, the band structure is very similar to the one of BLG on WSe2 Gmitra and Fabian (2017), where a spin-orbit valve effect has been proposed. In Fig. 10 we show the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{BLG}}$ for the BLG/Bi2Se3 stack as a function of a transverse electric field. We find that the dipole of the structure grows linearly with applied electric field, similar to the SLG case. Most of the orbital and spin-orbit parameters stay more or less constant as a transverse electric field is applied. However, the field can for example tune the doping energy of the topological insulator $E_{\textrm{TI}}$, and also slightly tune the interlayer hopping amplitude $\gamma_{3}$.

Surprisingly, the parameter $\lambda_{0}$ drastically decreases in magnitude for positive fields. This means that the electric field, felt by one graphene layer due to the presence of the other, in this proximity effect set-up, diminishes with applied external transverse field. Important for the previously mentioned spin-orbit valve effect is the closing of the orbital gap, and the subsequent reopening with inverted band structure. Indeed, such a situation can be realized, for a transverse electric field around 2.1 V/nm, as the two graphene layers are at the same potential, when the paramter $V$ goes through zero. For an electric field of 2 V/nm, the band structure is shown in Fig. 11(a). The orbital gap of the BLG bands is about 10 meV, and the conduction band is strongly spin-orbit split. For a field of 2.3 V/nm, the band structure is inverted, see Fig. 11(b), and the valence band is now strongly spin-orbit split. The spin-orbit split bands are always localized on the graphene layer closest to the substrate. In Fig. 2(f), we can see that the band splittings near the K point for conduction and valence band differ by orders of magnitude. Important for the spin-orbit valve effect is, that spin relaxation depends quadratically on the magnitude of the band splittings. Consequently, the spin relaxation for electrons (holes) is large (small) for an electric field of about 2 V/nm, and vice versa for a field of about 2.3 V/nm.

VI Additional considerations

Our DFT calculations, together with the extracted fit parameters, can give an insight on the magnitude of the induced SOC in SLG, and are helpful to analyze and interpret experimental data in SLG/topological insulator heterostructures, as recently proven Song et al. (2018); Jafarpisheh et al. (2018). However, there are still some open questions, that we would like to address in the following. Is the chosen interlayer distance of $3.5$ Å reasonable for the studied heterostructures? How does proximity SOC depend on it? How does the proximity SOC depend on the composition of the topological insulator? Can we tune the Dirac states of SLG down to the Fermi level, by using a different topological insulator?

VI.1 Distance study

It has been shown that proximity SOC and exchange effects can be significantly enhanced when decreasing the interlayer distance between materials Frank et al. (2016); Zollner et al. (2019); Zhang et al. (2018, 2015). Here, we look into the proximity effects, when we modify the interlayer distance between SLG and Bi2Se3. Similar results can be expected for the BLG case.

In Fig. 12 we show the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{SLG}}$ for the SLG/Bi2Se3 stack with 1QL as function of the interlayer distance between SLG and Bi2Se3. Our study shows, that the lowest energy is achieved for an interlayer distance of $3.4$ Å, very close to our chosen distance of $3.5$ Å according to Ref. Song et al. (2018). Note that the energetically most favorable interlayer distance depends on the specific DFT input and the chosen vdW corrections. The proximity-induced intrinsic and Rashba SOC parameters show the expected behavior and decrease in magnitude, as we increase the interlayer distance, see Fig. 12(d). The increase of the parameters is giant, about 200% when we decrease the interlayer distance by only about 10%.

Surprisingly, the PIA SOC parameters first increase with increasing the interlayer distance from $3.2~{}\textrm{\AA}$ to $3.4~{}\textrm{\AA}$. A possible explanation might be the position of the SLG Dirac point with respect to the bands of the topological insulator. Looking at the band structure Fig. 2(a), we see that there are bands of the topological insulator anti-crossing with the SLG Dirac bands, at around 0.6 eV. With increasing the distance, the Dirac point of SLG stays roughly at the same energy, while the bands of the topological insulator shift down in energy; compare the energies $E_{\textrm{D}}$ and $E_{\textrm{TI}}$ in Fig. 12(f). Consequently, the anti-crossing bands move closer towards the SLG Dirac point and we need larger PIA parameter to capture the low energy bands. For interlayer distances, larger than $3.4~{}\textrm{\AA}$, the two energies $E_{\textrm{D}}$ and $E_{\textrm{TI}}$ barely change anymore and the PIA SOC parameters decrease in magnitude as expected.

VI.2 Bi2Te2Se substrate

Experimentally, various atomic compositions are used to create the artificial topological insulator crystals Bi2-xSbxTe3-ySey, with some portions $x$ and $y$. The reason is that the unintentional intrinsic defect doping of the topological insulator can be compensated such that the Dirac surface states, with in-plane spin-momentum locking, are located near the Fermi level, and bulk transport can be suppressed Ren et al. (2011); Arakane et al. (2012). Recent transport measurements have shown giant proximity SOC in SLG on Bi1.5Sb0.5Te1.7Se1.3, of at least 2.5 meV Jafarpisheh et al. (2018). Here we show and discuss the first-principles calculated results for SLG on Bi2Te2Se, where we replace the outermost Se atoms of each QL by Te atoms in Fig. 1.

In Fig. 13, we show the full band structures of SLG above one, two, and three QLs of Bi2Te2Se. We find that the SLG Dirac point is now located near the Fermi level. Indeed, tuning the constituents of the topological insulator by $x$ and $y$, different band alignments can be formed. This case is also interesting for BLG, as the spin-orbit valve active bands could be shifted down to the Fermi level, by using a different topological insulator. In contrast to the case of Bi2Se3, the surface states of Bi2Te2Se are not yet gapless for 3QLs, but they are located much closer to the system Fermi level, compare Figs. 13 and 2.

In Tab. 3 we summarize the fit parameters of Hamiltonian $\mathcal{H}_{\textrm{SLG}}$ for the SLG/Bi2Te2Se stacks for different number of QLs. Again, we find that the fit parameters are almost independent on the number of QLs, indicating that only the closest QL is mainly responsible for proximity SOC. The proximity-induced intrinsic SOC parameters are about 20% larger, than for the Bi2Se3 substrate, and increase by about 10% for each QL that we add. The origin of the strong induced SOC in SLG is due to the nearest Se or Te atoms of the topological insulator. Because Te atoms have stronger atomic SOC than Se atoms, also the proximity induced SOC is enhanced.

VII Summary

In summary, we have shown that SLG and BLG on the topological insulator Bi2Se3 experience significant hole doping (350 meV for SLG, 200 meV for BLG) and strong proximity induced SOC of about 1-2 meV, giant compared to the small intrinsic SOC of the pristine SLG and BLG. Most surprising is, that the induced SOC is also of valley-Zeeman type, similar to TMDC substrates. The induced spin-orbit fields point mainly out-of-plane, even though the topological insulator hosts surface states with in-plane spin-momentum locking.

As we increase the number of QLs of Bi2Se3, below the proximitized SLG, the induced SOC can be increased by about 10%, each time a QL is added up to 3QLs. We expect the SOC to saturate as proximity induced phenomena are usually short range effects. In addition we show that an externally applied transverse electric field can tune band offsets and SOC parameters in SLG. This tunability of SOC by electric fields, can have significant impact on the spin relaxation properties. For the BLG, a spin-orbit valve effect can be realized, similar to the recently studied case of a TMDC substrate. In particular, we find that without applied electric field, the BLG states exhibit a band gap and a strongly spin-orbit split conduction band. For a moderate and experimentally accessible field of about 2.3 V/nm, the band structure can be inverted, now with strongly spin-orbit split valence band, offering the possibility to fully electrically control the magnitude of spin relaxation of electrons and holes in BLG.

Furthermore, in the case of SLG on 1QL of Bi2Se3, we have shown that a small decrease of the interlayer distance by only 10%, can strongly enhance proximity SOC by about 200%, allowing to tailor the magnitude of SOC by external pressure. Finally, we have extracted the orbital and SOC parameters for SLG on a different topological insulator Bi2Te2Se where the SLG is essentially undoped and still experiences strong valley-Zeeman SOC, even larger than in Bi2Se3. Experimentally, the tunability of the SLG doping level can be controlled by varying the constituents, $x$ and $y$, of the topological insulator Bi2-xSbxTe3-ySey. Tuning both Dirac states, from the SLG and the topological insulator, to the Fermi level allows to study the interplay of two very distinct spin-orbit fields at the same time.

Acknowledgements.

This work was supported by DFG SPP 1666, SFB 1277 (A09 and B07), the European Unions Horizon 2020 research and innovation program under Grant No. 785219.