Proximity band structure and spin textures on both sides of topological-insulator/ferromagnetic-metal interface and their transport probes J. M. Marmolejo-Tejada,1,2 K. Dolui,1 P. Lazi´c,3 P.-H. Chang,4 S. Smidstrup,5 D. Stradi,5 K. Stokbro,5 and B. K. Nikoli´c1,∗ 1Department of Physics and Astronomy, University of Delaware, Newark, DE 19716-2570, USA 2School of Electrical and Electronics Engineering, Universidad del Valle, Cali, AA 25360, Colombia 3Rudjer Boˇskovi´c Institute, PO Box 180, Bijeniˇcka c. 54, 10 002 Zagreb, Croatia 4Department of Physics and Astronomy, University of Nebraska Lincoln, Lincoln, Nebraska 68588, USA 7 5QuantumWise A/S, Fruebjergvej 3, Box 4, DK-2100 Copenhagen, Denmark 1 0 The control of recently observed spintronic effects in topological-insulator/ferromagnetic-metal 2 (TI/FM) heterostructures is thwarted by the lack of understanding of band structure and spin r texture around their interfaces. Here we combine density functional theory with Green’s function a techniques to obtain the spectral function at any plane passing through atoms of Bi Se and Co 2 3 M or Cu layers comprising the interface. In contrast to widely assumed but thinly tested Dirac cone gapped by the proximity exchange field, we find that the Rashba ferromagnetic model describes 8 the spectral function on the surface of Bi Se in contact with Co near the Fermi level E0, where 2 3 F 2 circularandsnowflake-likeconstantenergycontourscoexistaroundwhichspinlockstomomentum. The remnant of the Dirac cone is hybridized with evanescent wave functions injected by metallic ] layers and pushed, due to charge transfer from Co or Cu layers, few tenths of eV below E0 for l F l both Bi Se /Co and Bi Se /Cu interfaces while hosting distorted helical spin texture wounding a 2 3 2 3 h arounda singlecircle. Thesefeatures explainrecent observation[K.Kondou et al., Nat. Phys. 12, - 1027 (2016)] of sensitivity of spin-to-charge conversion signal at TI/Cu interface to tuning of EF0. s Interestingly, three monolayers of Co adjacent to Bi Se host spectral functions very different from 2 3 e thebulkmetal,aswellasin-planespintexturessignifyingthespin-orbitproximityeffect. Wepredict m thatout-of-planetunnelinganisotropicmagnetoresistanceinverticalheterostructureCu/Bi Se /Co, 2 3 . where current flowing perpendicular to its interfaces is modulated by rotating magnetization from t a paralleltoorthogonaltocurrentflow,canserveasasensitiveprobeofspintextureresidingatE0. F m - d Therecentexperimentsonspin-orbittorque(SOT)[1, (a) Bi2Se3 Lead Bi2Se3 (6 QL) Co (3 ML) n 2] and spin-to-charge conversion [3, 4] in topological- o Vacuum insulator/ferromagnetic-metal(TI/FM)heterostructures c 2 1 [ have ignited the field of topological spintronics. In these (b) x 2 devices, giant non-equilibrium spin densities [5–8] are Cu Lead Bi2Se3 (6 QL) Co Lead mϕ m expected to be generated due to strong spin-orbit cou- θ z v 2 pling (SOC) on metallic surfaces of three-dimensional y 2 1 3 4 6 (3D) TIs and the corresponding (nearly [9]) helical spin- V 4 momentum locking along a single Fermi circle for Dirac b 0 electrons hosted by those surfaces [10]. Such strong in- FIG. 1. Schematic view of TI-based heterostructures where: 0 terfacial SOC-driven phenomena are also envisaged to (a) semi-infinite Bi Se layer is attached to n monolayers of . 2 3 1 underlie a plethora of novel spintronic technologies [11]. Co(0001);(b)6QLsofBi2Se3 aresandwichedbetweensemi- 0 infinite Cu(111) layer and semi-infinite Co(0001) layer. Both 7 These effects have been interpreted almost exclusively heterostructures are infinite in the transverse direction, so 1 using simplistic models, such as the Dirac Hamiltonian that the depicted supercells are periodically repeated within v: for the TI surface with an additional Zeeman term de- the xy-plane. The magnetization m of the Co layer is fixed i scribingcouplingofmagnetizationoftheFMlayertothe along the z-axis in (a), or rotated within the xy-plane or the X surface state spins [10], HˆDirac =v (σˆ ×pˆ) −∆m·σˆ, xz-plane in (b). Applying the bias voltage Vb to the vertical r F z heterostructureinpanel(b)leadstoachargecurrentflowing wherepˆ isthemomentumoperator,σˆ isthevectorofthe a perpendicularlytobothBi Se /CuandBi Se /Cointerfaces. 2 3 2 3 Pauli matrices, m is the magnetization unit vector and v istheFermivelocity. Thus,theonlyeffectofFMlayer F captured by HˆDirac is proximity effect-induced exchange coupling ∆ which opens a gap in the Dirac cone energy- cannotbecapturedbysimplisticmodelslikeHˆDirac. The momentum dispersion [10], thereby making Dirac elec- propertiesofTI/FMinterfacesarefarmorecomplexdue trons massive. On the other hand, recent first-principles to injection of evanescent wave functions from the FM calculations [12, 13] demonstrate that band structure of layer into the bulk gap of the TI layer, which can hy- even TI/ferromagnetic-insulator (TI/FI) bilayers, where bridizewithsurfacestateofTIandbluritsDiraccone(as hybridizationbetweenTIandFIstatesislargelyabsent, already observed in tight-binding models of TI/FM in- 2 (a) n = 0 (b) n = 1 (c) n = 2 (d) n = 3 1 1 1 1 high 0.4 0.2 ) V e ( 0 0F E - E-0.2 -0.4 low (e) (f) (g) (h) 2 2 2 2 high 0.4 0.2 ) V e ( 0 0F E - E-0.2 -0.4 low X Γ Y X Γ Y X Γ Y X Γ Y FIG. 2. Spectral function, defined in Eq. (2), at plane 1 for panels (b)–(d) or plane 2 for panels (f)–(h) within Bi Se /Co(n 2 3 ML) heterostructure in Fig. 1(a) with m(cid:107)zˆ. For comparison, panels (a) and (e) plot the spectral function at planes 1 (akin to Ref. [15]) and 2, respectively, within semi-infinite Bi Se crystal in contact with vacuum (i.e., n=0). From Γ to Y we plot 2 3 A(E;k =0,k ;z∈{1,2}), while from Γ to X we plot A(E;k ,k =0;z∈{1,2}). x y x y terfaces [7, 14]), as well as related charge transfer. Thus, the TI surface [18]. thekeyissuefortopologicalspintronics[1–4,11]istoun- An attempt [20] to obtain the spectral function, derstandbandstructureandspintextures(includingthe A (E;k) = (cid:80)i∈QLjwi δ(E −ε ), directly from DFT fate of the Dirac cone and its spin-momentum locking) j n,i nk nk computed energy-momentum dispersion ε (n is the in hybridized TI with FM or normal metal (NM) [4] lay- nk band index and k is the crystal momentum) and site- ers at nanometer scale around the interface where they projected character wi of the corresponding eigenfunc- arebroughtintocontact,wherepropertiesofbothTIside nk tions for TI/FM supercells has produced ambiguous re- andFMorNMsideoftheinterfacecanbequitedifferent sults. Thisisduetoarbitrarinessinbroadeningthedelta from the properties of corresponding bulk materials. functionδ(E−ε ),aswellasduetousageofatomicsites nk iwithinthewholej quintuplelayer(QL )ofBi Se (one For example, computational searches [15] for new ma- j 2 3 QL consists of three Se layers strongly bonded to two Bi terials realizing 3D TIs (or other topologically non- layersinbetween)whicheffectivelyaveragesthespectral trivial electronic phases of matter like Weyl semi- function over all geometric planes within QL . Similar metals [16] and Chern insulators [17]) have crucially j ambiguities (such as setting the amount of electron den- relied on first-principles calculations of spectral func- sity which is localized on the surface or within the whole tion on their boundaries and its confirmation by spin- interfacial QL) plague interpretation of projected DFT and angle-resolved photoemission spectroscopy (spin- band structure of TI/FI [13] and TI/FM bilayers [21]. ARPES) [18]. A standard density functional theory (DFT)-based framework developed for this purpose— Herewedevelopaframeworkwhichcombinesthenon- where DFT band structure around the Fermi level E0 collinearDFTHamiltonianHDFT,representedinabasis F is reconstructed using the Wannier tight-binding Hamil- of variationally optimized localized atomic orbitals [22], tonian [19] used to obtain the retarded Green’s function with retarded GF calculations from which one can ex- (GF) of semi-infinite homogeneous crystal and the spec- tract the spectral function and spin textures at an arbi- tral function on its surface in contact with vacuum [15– trary geometric plane of interest within a junction com- 17]—is difficult to apply to complicated inhomogeneous bining TI, FM and NM layers. It also makes it possible systems like TI/FM bilayers due to strongly entangled to compute their spin and charge transport properties in bands in the region of interest around E0. Also, spin- the linear-response regime or at finite bias voltage. We F ARPES experiments cannot probe buried interfaces be- apply this framework to two Bi Se -based heterostruc- 2 3 lowtoomanymonolayers(e.g.,penetrationdepthoflow- tures whose supercells are depicted in Fig. 1, where we energy photons is 2–4 nm) of FM or NM deposited onto assume that those supercells are periodically repeated in 3 (a) (b) E-E0= 0 eV (c) E-E0= 0.05 eV (d) E-E0= -0.35 eV 1 F F F high1x10-3 0.2 0.5 0.1 ) V e 0(F 0 kx 0 zS E - E −0.1 -0.5 −0.2 -1 low-1x10-3 X Γ Y −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 k k k y y y (e) (f) E-E0= 0 eV (g) E-E0= 0.05 eV (h) E-E0= -0.35 eV 1 0.2 F F F high1x10-9 0.5 0.1 ) V e 0(F 0 kx 0 zS E - E -0.5 −0.1 -1 −0.2 low-1x10-9 X Γ Y −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 k k k y y y FIG. 3. Spectral function at: (a)–(d) plane 1 in Fig. 1(b) which is passing through Se atoms at the Bi Se /Co interface 2 3 with m(cid:107)zˆ; and (e)–(h) plane 2 in Fig. 1(b) which is passing through Se atoms at Bi Se /Cu interface, where we re- 2 3 move Co layer to make Bi Se semi-infinite along the z-axis. In panels (a) and (e), we plot A(E;k = 0,k ;z ∈ {1,2}) 2 3 x y from Γ to Y and A(E;k ,k = 0;z ∈ {1,2}) from Γ to X. Panels (b)–(d) and (f)–(h) plot constant energy contours of x y A(E−E0 ∈{0.0 eV,0.05 eV,−0.35 eV};k ,k ;z∈{1,2}) at three energies marked by horizontal dashed lines in panels (a) F x y or (e), respectively, as well as the corresponding spin textures where the out-of-plane S component is indicted in color (red z for positive and blue for negative). The units for k and k are 2π/a where a is the lattice constant of a common supercell x y combining two unit cells of the two layers around the corresponding interface. the transverse xy-direction. The heterostructure in Fig. 1(b) consists of semi- The heterostructure in Fig. 1(a) consists of Bi Se , infinite Cu and Co leads sandwiching a Bi Se layer of 2 3 2 3 chosenastheprototypical3DTI[9,10],whosesurfaceis finite thickness, where we choose Cu as the NM layer covered by n monolayers (MLs) of Co. The retarded GF similar to the very recent spin-to-charge conversion ex- of this heterostructure is computed as periment of Ref. [4]. Such a heterostructure is termed vertical or current-perpendicular-to-plane in spintronics G (E)=[E−HDFT−ΣBi2Se3(E)]−1, (1) k(cid:107) k(cid:107) k(cid:107) terminology since applying bias voltage Vb drives a cur- rentperpendicularlytotheTI/FMinterface. Itsretarded where k = (k ,k ) is the transverse k-vector, (cid:107) x y GF is computed as ΣBi2Se3(E) is the self-energy [23] describing the semi- inkfi(cid:107)nite Bi2Se3 lead and HDk(cid:107)FT is the Hamiltonian of the Gk(cid:107)(E)=[E−HDk(cid:107)FT−ΣCk(cid:107)u(E)−ΣCk(cid:107)o(E)]−1, (3) active region consisting of n MLs of cobalt plus 6 QLs of where HDFT describes the active region consisting of 6 Bi2Se3 to which the lead is attached. We choose n=1–3 k(cid:107) QLs of Bi Se plus 4 MLs of Cu and 4 MLs of Cu. Its since ultrathin FM layers of thickness (cid:39)1 nm are typi- 2 3 linear-responseresistanceRisgivenbytheLandauerfor- cally employed in SOT experiments [24] in order to pre- mula serveperpendicularmagneticanisotropy(notethatmag- netocrystalline anisotropy does favor out-of-plane m in 1 e2 (cid:90) (cid:90) (cid:18) ∂f (cid:19) = dk dE − Tr[ΓCoG ΓCuG† ], Bi2Se3/Co bilayers [21]). The spectral function (or lo- R hΩBZ BZ (cid:107) ∂E k(cid:107) k(cid:107) k(cid:107) k(cid:107) cal density of states) at an arbitrary plane at position z (4) within the active region is computed from where we assume temperature T = 300 K in the Fermi- Dirac distribution function f(E), Γα =i(Σα −[Σα ]†) A(E;kx,ky,z)=−Im[Gk(cid:107)(E;z,z)]/π, (2) andΩ istheareaofthetwo-dimenks(cid:107)ional(2kD(cid:107))Brillko(cid:107)uin BZ where the diagonal matrix elements G (E;z,z) are ob- zone (BZ) within which k vectors are sampled. k(cid:107) (cid:107) tained by transforming Eq. (1) from orbital to a real- The spectral function of the heterostructure in space representation. Fig. 1(a) computed at planes 1 and 2 within the Bi Se 2 3 4 (a) (b) E-E0= 0 eV (c) E-E0= 0.05 eV (d) E-E0= -0.35 eV 1 F F F high1x10-1 0.2 0.5 0.1 ) V e 0(F 0 kx 0 zS E - E −0.1 -0.5 −0.2 -1 low-1x10-1 X Γ Y −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 k k k y y y (e) (f) E-E0= 0 eV (g) E-E0= 0.05 eV (h) E-E0= -0.35 eV 1 y F F F high1x10-2 0.2 0.5 0.1 ) V e 0(F 0 kx 0 zS E - E −0.1 -0.5 −0.2 -1 low-1x10-2 X Γ Y −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 −0.2 −0.1 0 0.1 0.2 k k k y y y FIG. 4. Spectral function at: (a)–(d) the surface of a semi-infinite Co layer in contact with vacuum where m is perpendicular to the surface; and (e)–(h) plane 3 in Fig. 1(b) which is passing through Co atoms at the Bi Se /Co interface with m(cid:107)zˆ. In 2 3 panel (e), we plot A(E;k = 0,k ;z ∈ 3) from Γ to Y and A(E;k ,k = 0;z ∈ 3) from Γ to X. Panels (b)–(d) and (f)–(h) x y x y plot constant energy contours of the spectral function at three energies marked by horizontal dashed lines in panels (a) or (e), respectively, as well as the corresponding spin textures where the magnitude of the out-of-plane S component is indicted in z color (red for positive and blue for negative). The units for k and k are 2π/a where a is the lattice constant of a common x y supercell combining two unit cells of the two layers around the corresponding interface. layer is shown in Figs. 2(b)–(d) and 2(f)–(h), respec- For infinitely many MLs of Co attached to 6 QLs tively, where plane 1 is passing through Se atoms on the of Bi Se within the Cu/Bi Se /Co heterostructure in 2 3 2 3 Bi Se surface in contact with Co layer and plane 2 is Fig. 1(b), the remnant of the Dirac cone from the TI 2 3 three QLs (or (cid:39)2.85 nm) away from plane 1. For com- surface can be identified in Fig. 3(a) at around 0.5 eV parison,wealsoshowinFigs.2(a)and 2(b)thespectral below E0 while it is pushed even further below in the F function at the same two planes within the semi-infinite case of Cu/Bi Se interface in Fig. 3(e). The differ- 2 3 Bi Se layer in contact with vacuum, thereby reproduc- ence in work functions Φ =5.0 eV or Φ =4.7 eV 2 3 Co Cu ing the results from Ref. [15] by our formalism. While andelectronaffinityχ =5.3 eVdetermines[20]the Bi2Se3 the Dirac cone at the Γ-point is still intact in Fig. 2(b) bandalignmentandthestrengthofhybridization, where for n = 1 ML of Co, its Dirac point (DP) is gradually n-type doping [see also Figs. 6(c) and 6(d)] of the pushedintothevalencebandofBi Se withincreasingn Bi Se layer pins E0 of the whole Cu/Bi Se /Co het- 2 3 2 3 F 2 3 because of charge transfer from metal to TI. The charge erostructure in the conduction band of the bulk Bi Se . 2 3 transfer visualized in Figs. 6(c) and 6(d) is relatively The remnant of the Dirac cone is quite different from small, but due to small density of states (DOS) at the the often assumed [10] eigenspectrum of HˆDirac be- DP it is easy to push it down until it merges with the cause of hybridization with the valence band of Bi Se , 2 3 larger DOS in the valence band of the TI. Adding more as well as with states injected by the Co or Cu lay- MLs of Co in Figs. 2(c) and 2(d) also introduces addi- ers whose penetration into TI is visualized by plot- tional bands within the bulk gap of Bi Se due to injec- ting position- and energy-dependent spectral function 2 3 tion of evanescent wave functions which hybridize with A(E;z) = 1 (cid:82)dk dk A(E;k ,k ;z) in Fig. 5(a). On ΩBZ x y x y the Dirac cone. The metallic surface states of Bi Se it- the other hand, the energy-momentum dispersion in the 2 3 selfpenetrateintoitsbulkoveradistanceof2QLs[6],so vicinity of E0 and for an interval of k vectors around F (cid:107) thatinFig.2(e)thespectralfunctiononplane2vanishes theΓ-pointissurprisinglywell-describedbyanothersim- insidethegapofthesemi-infiniteBi Se layerincontact plistic model—ferromagnetic Rashba Hamiltonian [25]. 2 3 with vacuum, while the remaining states inside the gap in Figs. 2(f)–(h) can be attributed to the Co layer. In Figs. 3(b)—(d) and Figs. 3(f)—(h) we show con- stant energy contours of the spectral function at three 5 selected energies E denoted in Figs. 3(a) and 3(e) by (a) 1Cu Bi2Se3 Cohigh (b) dashed horizontal lines. Instead of a single circle as 60 %)1 tHttˆhohDeehierecaixocga,negnsootrsanpsnaeitnlctgwerlnueaemrhpregioxnyfaggHc[oo9Dnn])FtooTsruusroffinfofctowiherflentatihlksyoeel-aalewitkiegaedeycnBofsrnpio2temSoceut3rrDuslPma(ydfeuoorerf, 0E−E(eV)F−00..055 TAMR (%)out23450000 TAMR (in-010 9EEE0AFFn---gEEEl000eFF1 8===ϕ0 (-00d0e.e0.gV325)75e0eVV360 here we find multiple circular and snowflake-like con- 10 F F tours close to the Γ-point. The spin textures within the −1 low 0 10 20 30 40 50 60 70 -90 -60 -30 0 30 60 90 constant energy contours are computed from the spin- z (Å) Angle θ (deg) resolvedspectralfunction. ForenergiesnearE =E0,the F spin textures shown in Figs. 3(b) and 3c) are quite dif- FIG.5. (a)Theposition-andenergy-dependentspectralfunc- tion A(E;z) = 1 (cid:82)dk dk A(E;k ,k ;z) from the left Cu ferent from the helical ones in isolated Bi2Se3 layer [15]. ΩBZ x y x y lead, across Bi Se tunnel barrier, toward the right Co lead Nevertheless, Fig. 3(d) shows that the remnant Dirac 2 3 for the heterostructure in Fig. 1(b). (b) The out-of-plane cone still generates distorted helical spin texture wound- TAMR (θ) ratio defined in Eq. (5) as function of angle θ out ing along a single circle but with out-of-plane S compo- z betweenthemagnetizationmandthedirectionofcurrentin- nent due to the presence of Co layer. jected along the z-axis in Fig. 1(b). Inset in panel (b) shows The envisaged applications of TIs in spintronics are angular dependence of the in-plane TAMR (φ) ratio. In or- in based[1–5,7,8]onspintexturesliketheoneinFig.3(d) dertoconvergetheintegrationoverthetransversewavevector k inEq.(4),weemployauniformgridof101×101k-points since it maximizes [5, 6] generation of nonequilibrium (cid:107) for TAMR (θ) and 251×251 k-points for TAMR (φ). spin density when current is passed parallel to the TI out in surface. However, utilizing spin texture in Fig. 3(d) in lateralTI/FMheterostructureswouldrequiretoshiftE F by the charge distribution at the metal/vacuum inter- (by changing the composition of TI [4] or by applying face and thereby confine wave functions into a Rashba a gate voltage [2]) by few tenths of eV below E0 of F spin-split quasi-2D electron gas [28]. Thus, Figs. 4(a)– the undoped heterostructures in Fig. 3(a). For example, (d) explains the origin of recently observed [29] SOT in extreme sensitivity of spin-to-charge conversion was re- theabsenceofanyadjacentheavymetalorTIlayersince centlyobservedRef.[4]onthesurfaceof(Bi Sb ) Te 1−x x 2 3 passingcurrentparalleltoMLsofCohostingnonzeroin- TIcoveredbya8nmthickCulayerasE oftheTIlayer F plane spin textures will generate a nonequilibrium spin was tuned, which is difficult to explain by assuming that density[6]S andspin-orbittorque∝S ×m[8,26]. theDiracconeontheTIsurfaceremainsintactafterthe neq neq Finally,weproposeapurelychargetransportmeasure- deposition of the Cu layer (e.g., Ref. [4] had to invoke ment that could detect which among the spin-textures “instability of the helical spin structure”). On the other shown in Figs. 3(b)–(d) resides at the Fermi level of hand, it is easy to understand from Figs. 3(f)–(h) in- TI/FMinterface. Ourschemerequirestofabricateverti- terface demonstrating how spin textures at Bi Se /Cu 2 3 calheterostructureinFig.1(b)andmeasureitstunneling interface change dramatically as one moves E (even F anisotropic magnetoresistance (TAMR). The TAMR is a slightly) below or above E0. Comparing Figs. 3(a)–(d) F phenomenon observed in magnetic tunnel junctions with with 3(e)–(h) makes it possible to understand the effect a single FM layer [7, 27, 30, 31], where SOC makes the of the magnetization of the Co layer, which modifies [25] band structure anisotropic so that the resistance of such RashbadispersionaroundE0 andthecorrespondingspin F junctions changes as the magnetization m is rotated by textures (particularly the out-of-plane S component). z angle θ or φ in Fig. 1(b). The resistance change is quan- The theoretical modeling of SOT in TI/FM [1, 8] or tified by the TAMR ratio defined as [27, 31] heavy-metal/FM [26] bilayers is usually conducted by starting from strictly 2D Hamiltonians, such as HˆDirac R(α)−R(0) TAMR (α)= . (5) or the Rashba ferromagnetic model [25], respectively, out(in) R(0) so that the FM layer is not considered explicitly. Fig- ures 4(e)–(h) show that this is not warranted since the Here α ≡ θ for TAMR where magnetization in out Bi Se layerinducesproximitySOCandthecorrespond- Fig.1(b)rotatesintheplaneperpendiculartotheTI/FM 2 3 ing in-plane spin texture on the first ML of Co, which interface,andα≡φforTAMR wheremagnetizationin in decays to zero only after reaching plane 4 in Fig. 1(b). Fig.1(b)rotateswithintheplaneoftheTI/FMinterface. In fact, we find non-trivial in-plane spin texture even In the case of TAMR (θ), R(0) is the resistance when out on the surface of Co in contact with vacuum, as shown m (cid:107) zˆ in Fig. 1; and in the case of TAMR (φ), R(0) is in in Figs. 4(b)–(d), which is nevertheless quite different the resistance when m(cid:107)xˆ in Fig. 1. Thus, TAMR (θ) out from those in Figs. 4(f)–(h). The in-plane spin texture changes due to the different orientations of the magneti- in Figs. 4(b)–(d) is a consequence of the Rashba SOC zation with respect to the direction of the current flow, enabled by inversion asymmetry due to Co surface [27] while the situation becomes more subtle for TAMR (φ) in where an electrostatic potential gradient can be created where the magnetization remains always perpendicular 6 (a) (b) (c) (d) (e/Å3) FIG. 6. (Color online). Top and side view of common unit cells for (a) Bi Se /Cu(111) and (b) Bi Se /Co(0001) bilayers. 2 3 2 3 Panels (c) and (d) show charge rearrangement around the interface of bilayers in panels (a) and (b), respectively. to the current flow. Figure 5(b) demonstrates that the In order to avoid interaction with periodic images of the largest TAMR (θ = ±90◦) is obtained by tuning the bilayer, 18 ˚A of vacuum was added in the z–direction. out Fermi level to E −E0 =−0.35 eV so that nearly heli- F F ForthecaseofBi Se onCo(0001),themostfavorable cal spin texture in Fig. 3(d) resides at the Fermi level. 2 3 positionyieldsabindingenergyof460meVperCoatom. Another signature of its presence is rapid increase of Both ML of Co and QL of Bi Se in direct contact gain TAMR (θ)whentiltingmbysmallanglesθ awayfrom 2 3 out somecorrugation, roughlyaround(cid:39)0.1 ˚A, whiletheav- the current direction. The in-plane TAMR (φ) shown in erage z–distance between them is 2.15 ˚A. The average in the inset of Fig. 5(b) is much smaller (and difficult to distance between the ML of Cu and QL of Bi Se in di- converge in the number of transverse k-points) quantity 2 3 rect contact is around 2.26 ˚A with smaller corrugation whichdoesnotdifferentiatebetweenspintexturesshown than in the case of Co(0001), while the binding energy in Figs. 3(b)–(d). is 294 meV per Cu atom. For other relative positions of Bi Se layerwithrespecttoCu(111)andCo(0001)layers 2 3 the difference in binding energy is very small. Binding METHODS energies in both cases are rather small, thereby signal- ingthedominantvdWforces. Nevertheless,somecharge We employed the interface builder in the VNL [39] and rearrangement does occur at the interface due to push CellMatch[32]packagestoconstructacommonunitcells back/pillow effect [37], as shown in Figs. 6(c) and 6(d) for: (a) Bi Se /Cu(111) bilayer, where the common unit 2 3 where charge rearrangement is more pronounced for the cell is 5 × 5 in size compared to the smallest possible case of Bi Se /Cu(111) interface. 2 3 Cu(111) slab cell and copper is under compressive strain ThecalculationoftheretardedGFinEqs.(1)and (3) of 0.9 % while Bi Se lattice constant is unchanged; (b) 2 3 requires HDFT represented in the linear combination of Bi2Se3/Co(0001) bilayer where Co(0001) has the same k(cid:107) lattice constant as Bi Se , so the same unit cell as for atomic orbitals (LCAO) basis set which makes it possi- 2 3 Cu(111) is used without any strain on Co(0001). These bletospatiallyseparatesystemintotheactiveregionat- two unit cells are illustrated in Figs. 6(a) and 6(b), re- tached to one or two semi-infinite leads, as illustrated in spectively. In order to determine the best stacking of Figs. 1(a) and 1(b), respectively. We employ ATK pack- atomic layers and the distance of Bi Se atoms with re- age [38] for pseudopotential-based LCAO noncollinear 2 3 spect to surfaces of Cu(111) and Co(0001), we use DFT DFT calculations yielding HDFT, from which we obtain k(cid:107) calculations as implemented in the VASP package [33]. retarded GFs and the corresponding spectral functions, The electron core interactions are described by the pro- as well as the resistance in Eq. (4). In ATK calculations, jector augmented wave (PAW) method [34], and vdW- we use Perdew-Burke-Ernzerhof (PBE) parametrization DF [35] with optB88 is used as density functional [36] of generalized gradient approximation for the exchange- inordertodescribevanderWaals(vdW)forcesbetween correlationfunctional;norm-conservingpseudopotentials QLsofBi Se orbetweenBi Se andmetalliclayers. The for describing electron-core interactions; and LCAO ba- 2 3 2 3 cutoffenergyfortheplanewavebasissetis520eVforall sis set generated by the OpenMX package [22, 40] which calculations, while k-points were sampled at 3×3 surface consists of s2p2d1 orbitals on Co, Cu and Se atoms, and mesh. We use Cu and Co layers consisting of 5 MLs, s2p2d2 on Bi atoms. These pseudoatomic orbitals were where 3 bottom MLs are fixed at bulk positions while generated by a confinement scheme [22] with the cutoff the top two metallic MLs closest to Bi Se are allowed radius7.0and8.0a.u. forSeandBiatoms, respectively, 2 3 tofullyrelaxuntilforcesonatomsdropbelow1 meV/˚A. and 6.0 a.u. for Co and Cu atoms. The energy mesh 7 cutoff for the real-space grid is chosen as 75.0 Hartree. 1057 (2011). [11] A. Soumyanarayanan, N. Reyren, A. Fert, and C. Panagopoulos, Nature 539, 509 (2016). ACKNOWLEDGMENTS [12] W. Luo and X.-L. Qi, Phys. Rev. B 87, 085431 (2013). [13] A. T. Lee, M. J. Han, and K. Park, Phys. Rev. B 90, 155103 (2014). We thank Takafumi Sato and Jia Zhang for insight- [14] E. Zhao, C. Zhang, and M. Lababidi, Phys. Rev. B 82, ful discussions. J. M. M.-T., K. D. and B. K. N. were 205331 (2010). supported by NSF Grant No. ECCS 1509094. J. M. [15] H. Zhang et al., Nat. Phys. 5, 438 (2009). M.-T. also acknowledges support from Colciencias (De- [16] A. A. Soluyanov et al., Nature 527, 495 (2015). partamento Administrativo de Ciencia, Tecnologia e In- [17] X.-L.ShengandB.K.Nikoli´c,arXiv:1610.02719(2016). [18] T. Shoman et al., Nat. Commun. 6, 6547 (2015). novacion) of Colombia. P. L. was supported by the [19] N. Marzari et al., Rev. Mod. Phys. 84, 1419 (2012). Unity Through Knowledge Fund, Contract No. 22/15 [20] C.D.SpataruandF.L´eonard,Phys.Rev.B90,085115 and H2020 CSA Twinning Project No. 692194, RBI- (2014). T-WINNING. S. Smidstrup, D. Stradi and K. Stokbro [21] J.Zhang,J.P.Velev,X.Dang,andE.Y.Tsymbal,Phys. acknowledge support from the European Commission Rev. B 94, 014435 (2016). Seventh Framework Programme Grant Agreement IIIV- [22] T. Ozaki, Phys. Rev. B 67, 155108 (2003). MOS,ProjectNo. 61932;andHorizon2020researchand [23] J. Velev and W. Butler, J. Phys.: Condens. Matter 16, innovation programme under grant agreement SPICE, R637 (2004); I. Rungger and S. Sanvito, Phys. Rev. B 78, 035407 (2008); H. H. B. Sørensen et al., Phys. Rev. project No 713481. The supercomputing time was pro- B 77, 155301 (2008). vided by XSEDE, which is supported by NSF Grant No. [24] J. Kim et al., Nat. Mater. 12, 240 (2013). ACI-1053575. [25] N. Nagaosa et al., Rev. Mod. Phys. 82, 1539 (2010). [26] P. M. Haney et al., Phys. Rev. B 87, 174411 (2013). [27] A. N. Chantis, K. D. Belashchenko, E. Y. Tsymbal, and M.vanSchilfgaarde,Phys.Rev.Lett.98,046601(2007). [28] M. S. Bahramy et al., Nat. Commun. 3, 1159 (2012). ∗ [email protected] [29] S. Emori et al., Phys. Rev. B 93, 180402 (2016). [1] A. R. Mellnik et al., Nature 511, 449 (2014); Y. Fan [30] B. G. Park et al., Phys. Rev. Lett. 100, 087204 (2008). et al.,Nat.Mater.13,699(2014);Y.Wanget al.,Phys. [31] A. Matos-Abiague and J. Fabian, Phys. Rev. B 79, Rev. Lett. 114, 257202 (2015). 155303 (2009). [2] Y. Fan et al., Nat Nano 11, 352 (2016). [32] P. Lazi´c, Comp. Phys. Commun. 197, 324 (2015). [3] Y. Shiomi et al., Phys. Rev. Lett. 113, 196601 (2014); [33] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993); P. Deorani et al., Phys. Rev. B 90, 094403 (2014); M. G. Kresse and J. Furthmu¨llerb, Comput. Mater. Sci. 6, Jamali et al., Nano Lett. 15, 7126 (2015); J.-C. Rojas- 15 (1996). Sa´nchez et al., Phys. Rev. Lett. 116, 096602 (2016); H. [34] P. E. Blo¨chl, Phys. Rev. B 50, 17953 (1994); G. Kresse Wang et al., Phys. Rev. Lett. 117, 076601 (2016). and D. Joubert, Phys. Rev. B 59, 1758 (1999). [4] K. Kondou et al., Nat. Phys. 12, 1027 (2016). [35] M. Dion et al., Phys. Rev. Lett. 92, 246401 (2004). [5] D. Pesin and A. H. MacDonald, Nature Mater. 11, 409 [36] F. Mittendorfer et al., Phys. Rev. B 84, 201401 (2011). (2012). [37] H. Va´zquez, Y. J. Dappe, J. Ortega, and F. Flores, J. [6] P.-H. Chang et al., Phys. Rev. B 92, 201406(R) (2015). Chem. Phys. 126, 144703 (2007). [7] F. Mahfouzi, N. Nagaosa, and B. K. Nikoli´c, Phys. Rev. [38] Atomistix ToolKit (ATK) 2016.3, http://www. B 90, 115432 (2014). quantumwise.com. [8] F. Mahfouzi, B. K. Nikoli´c, and N. Kioussis, Phys. Rev. [39] Virtual NanoLab (VNL) 2016.3, http://www. B 93, 115419 (2016). quantumwise.com. [9] A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, [40] OpenMX 3.8, http://www.openmx-square.org/. 021004 (2016). [10] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010); X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83,