# Proximity band structure and spin textures on both sides of topological-insulator/ferromagnetic-metal interface and their transport probes

###### Abstract

The control of recently observed spintronic effects in topological-insulator/ferromagnetic-metal (TI/FM) heterostructures is thwarted by the lack of understanding of band structure and spin texture around their interfaces. Here we combine density functional theory with Green’s function techniques to obtain the spectral function at any plane passing through atoms of BiSe and Co 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 the spectral function on the surface of BiSe in contact with Co near the Fermi level , where circular and snowflake-like constant energy contours coexist around which spin locks to momentum. 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 for both BiSe/Co and BiSe/Cu interfaces while hosting distorted helical spin texture wounding around a single circle. These features explain recent 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 . Interestingly, three monolayers of Co adjacent to BiSe host spectral functions very different from the bulk metal, as well as in-plane spin textures signifying the spin-orbit proximity effect. We predict that out-of-plane tunneling anisotropic magnetoresistance in vertical heterostructure Cu/BiSe/Co, where current flowing perpendicular to its interfaces is modulated by rotating magnetization from parallel to orthogonal to current flow, can serve as a sensitive probe of spin texture residing at .

The recent experiments on spin-orbit torque (SOT) Mellnik2014 (); Fan2016 () and spin-to-charge conversion Shiomi2014 (); Kondou2016 () in topological-insulator/ferromagnetic-metal (TI/FM) heterostructures have ignited the field of topological spintronics. In these devices, giant non-equilibrium spin densities Pesin2012 (); Chang2015 (); Mahfouzi2014a (); Mahfouzi2016 () are expected to be generated due to strong spin-orbit coupling (SOC) on metallic surfaces of three-dimensional (3D) TIs and the corresponding (nearly Bansil2016 ()) helical spin-momentum locking along a single Fermi circle for Dirac electrons hosted by those surfaces Hasan2010 (). Such strong interfacial SOC-driven phenomena are also envisaged to underlie a plethora of novel spintronic technologies Soumyanarayanan2016 ().

These effects have been interpreted almost exclusively using simplistic models, such as the Dirac Hamiltonian for the TI surface with an additional Zeeman term describing coupling of magnetization of the FM layer to the surface state spins Hasan2010 (), , where is the momentum operator, is the vector of the Pauli matrices, is the magnetization unit vector and is the Fermi velocity. Thus, the only effect of FM layer captured by is proximity effect-induced exchange coupling which opens a gap in the Dirac cone energy-momentum dispersion Hasan2010 (), thereby making Dirac electrons massive. On the other hand, recent first-principles calculations Luo2013 (); Lee2014a () demonstrate that band structure of even TI/ferromagnetic-insulator (TI/FI) bilayers, where hybridization between TI and FI states is largely absent, cannot be captured by simplistic models like . The properties of TI/FM interfaces are far more complex due to injection of evanescent wave functions from the FM layer into the bulk gap of the TI layer, which can hybridize with surface state of TI and blur its Dirac cone (as already observed in tight-binding models of TI/FM interfaces Mahfouzi2014a (); Zhao2010 ()), as well as related charge transfer. Thus, the key issue for topological spintronics Mellnik2014 (); Fan2016 (); Shiomi2014 (); Kondou2016 (); Soumyanarayanan2016 () is to understand band structure and spin textures (including the fate of the Dirac cone and its spin-momentum locking) in hybridized TI with FM or normal metal (NM) Kondou2016 () layers at nanometer scale around the interface where they are brought into contact, where properties of both TI side and FM or NM side of the interface can be quite different from the properties of corresponding bulk materials.

For example, computational searches Zhang2009 () for new materials realizing 3D TIs (or other topologically nontrivial electronic phases of matter like Weyl semi-metals Soluyanov2015 () and Chern insulators Sheng2016 ()) have crucially relied on first-principles calculations of spectral function on their boundaries and its confirmation by spin- and angle-resolved photoemission spectroscopy (spin-ARPES) Shoman2015 (). A standard density functional theory (DFT)-based framework developed for this purpose—where DFT band structure around the Fermi level is reconstructed using the Wannier tight-binding Hamiltonian Marzari2012 () used to obtain the retarded Green’s function (GF) of semi-infinite homogeneous crystal and the spectral function on its surface in contact with vacuum Zhang2009 (); Soluyanov2015 (); Sheng2016 ()—is difficult to apply to complicated inhomogeneous systems like TI/FM bilayers due to strongly entangled bands in the region of interest around . Also, spin-ARPES experiments cannot probe buried interfaces below too many monolayers (e.g., penetration depth of low-energy photons is 2–4 nm) of FM or NM deposited onto the TI surface Shoman2015 ().

An attempt Spataru2014 () to obtain the spectral function, , directly from DFT computed energy-momentum dispersion ( is the band index and is the crystal momentum) and site-projected character of the corresponding eigenfunctions for TI/FM supercells has produced ambiguous results. This is due to arbitrariness in broadening the delta function , as well as due to usage of atomic sites within the whole quintuple layer (QL) of BiSe (one QL consists of three Se layers strongly bonded to two Bi layers in between) which effectively averages the spectral function over all geometric planes within QL. Similar ambiguities (such as setting the amount of electron density which is localized on the surface or within the whole interfacial QL) plague interpretation of projected DFT band structure of TI/FI Lee2014a () and TI/FM bilayers Zhang2016 ().

Here we develop a framework which combines the noncollinear DFT Hamiltonian , represented in a basis of variationally optimized localized atomic orbitals Ozaki2003 (), with retarded GF calculations from which one can extract the spectral function and spin textures at an arbitrary geometric plane of interest within a junction combining TI, FM and NM layers. It also makes it possible to compute their spin and charge transport properties in the linear-response regime or at finite bias voltage. We apply this framework to two BiSe-based heterostructures whose supercells are depicted in Fig. 1, where we assume that those supercells are periodically repeated in the transverse -direction.

The heterostructure in Fig. 1(a) consists of BiSe, chosen as the prototypical 3D TI Bansil2016 (); Hasan2010 (), whose surface is covered by monolayers (MLs) of Co. The retarded GF of this heterostructure is computed as

(1) |

where is the transverse -vector, is the self-energy Velev2004 () describing the semi-infinite BiSe lead and is the Hamiltonian of the active region consisting of MLs of cobalt plus 6 QLs of BiSe to which the lead is attached. We choose 1–3 since ultrathin FM layers of thickness nm are typically employed in SOT experiments Kim2013 () in order to preserve perpendicular magnetic anisotropy (note that magnetocrystalline anisotropy does favor out-of-plane in BiSe/Co bilayers Zhang2016 ()). The spectral function (or local density of states) at an arbitrary plane at position within the active region is computed from

(2) |

where the diagonal matrix elements are obtained by transforming Eq. (1) from orbital to a real-space representation.

The heterostructure in Fig. 1(b) consists of semi-infinite Cu and Co leads sandwiching a BiSe layer of finite thickness, where we choose Cu as the NM layer similar to the very recent spin-to-charge conversion experiment of Ref. Kondou2016 (). Such a heterostructure is termed vertical or current-perpendicular-to-plane in spintronics terminology since applying bias voltage drives a current perpendicularly to the TI/FM interface. Its retarded GF is computed as

(3) |

where describes the active region consisting of 6 QLs of BiSe plus MLs of Cu and MLs of Cu. Its linear-response resistance is given by the Landauer formula

(4) |

where we assume temperature K in the Fermi-Dirac distribution function , and is the area of the two-dimensional (2D) Brillouin zone (BZ) within which vectors are sampled.

The spectral function of the heterostructure in Fig. 1(a) computed at planes 1 and 2 within the BiSe layer is shown in Figs. 2(b)–(d) and 2(f)–(h), respectively, where plane 1 is passing through Se atoms on the BiSe surface in contact with Co layer and plane 2 is three QLs (or nm) away from plane 1. For comparison, we also show in Figs. 2(a) and 2(b) the spectral function at the same two planes within the semi-infinite BiSe layer in contact with vacuum, thereby reproducing the results from Ref. Zhang2009 () by our formalism. While the Dirac cone at the -point is still intact in Fig. 2(b) for ML of Co, its Dirac point (DP) is gradually pushed into the valence band of BiSe with increasing because of charge transfer from metal to TI. The charge transfer visualized in Figs. 6(c) and 6(d) is relatively small, but due to small density of states (DOS) at the DP it is easy to push it down until it merges with the larger DOS in the valence band of the TI. Adding more MLs of Co in Figs. 2(c) and 2(d) also introduces additional bands within the bulk gap of BiSe due to injection of evanescent wave functions which hybridize with the Dirac cone. The metallic surface states of BiSe itself penetrate into its bulk over a distance of 2 QLs Chang2015 (), so that in Fig. 2(e) the spectral function on plane 2 vanishes inside the gap of the semi-infinite BiSe layer in contact with vacuum, while the remaining states inside the gap in Figs. 2(f)–(h) can be attributed to the Co layer.

For infinitely many MLs of Co attached to 6 QLs of BiSe within the Cu/BiSe/Co heterostructure in Fig. 1(b), the remnant of the Dirac cone from the TI surface can be identified in Fig. 3(a) at around 0.5 eV below while it is pushed even further below in the case of Cu/BiSe interface in Fig. 3(e). The difference in work functions eV or eV and electron affinity eV determines Spataru2014 () the band alignment and the strength of hybridization, where -type doping [see also Figs. 6(c) and 6(d)] of the BiSe layer pins of the whole Cu/BiSe/Co heterostructure in the conduction band of the bulk BiSe. The remnant of the Dirac cone is quite different from the often assumed Hasan2010 () eigenspectrum of because of hybridization with the valence band of BiSe, as well as with states injected by the Co or Cu layers whose penetration into TI is visualized by plotting position- and energy-dependent spectral function in Fig. 5(a). On the other hand, the energy-momentum dispersion in the vicinity of and for an interval of vectors around the -point is surprisingly well-described by another simplistic model—ferromagnetic Rashba Hamiltonian Nagaosa2010 ().

In Figs. 3(b)—(d) and Figs. 3(f)—(h) we show constant energy contours of the spectral function at three selected energies denoted in Figs. 3(a) and 3(e) by dashed horizontal lines. Instead of a single circle as the constant energy contour for the eigenspectrum of , or single hexagon or snowflake-like contours (due to hexagonal warping Bansil2016 ()) sufficiently away from DP for the eigenspectrum of of the isolated BiSe layer, here we find multiple circular and snowflake-like contours close to the -point. The spin textures within the constant energy contours are computed from the spin-resolved spectral function. For energies near , the spin textures shown in Figs. 3(b) and 3c) are quite different from the helical ones in isolated BiSe layer Zhang2009 (). Nevertheless, Fig. 3(d) shows that the remnant Dirac cone still generates distorted helical spin texture wounding along a single circle but with out-of-plane component due to the presence of Co layer.

The envisaged applications of TIs in spintronics are based Mellnik2014 (); Fan2016 (); Shiomi2014 (); Kondou2016 (); Pesin2012 (); Mahfouzi2014a (); Mahfouzi2016 () on spin textures like the one in Fig. 3(d) since it maximizes Pesin2012 (); Chang2015 () generation of nonequilibrium spin density when current is passed parallel to the TI surface. However, utilizing spin texture in Fig. 3(d) in lateral TI/FM heterostructures would require to shift (by changing the composition of TI Kondou2016 () or by applying a gate voltage Fan2016 ()) by few tenths of eV below of the undoped heterostructures in Fig. 3(a). For example, extreme sensitivity of spin-to-charge conversion was recently observed Ref. Kondou2016 () on the surface of (BiSb)Te TI covered by a 8 nm thick Cu layer as of the TI layer was tuned, which is difficult to explain by assuming that the Dirac cone on the TI surface remains intact after the deposition of the Cu layer (e.g., Ref. Kondou2016 () had to invoke “instability of the helical spin structure”). On the other hand, it is easy to understand from Figs. 3(f)–(h) interface demonstrating how spin textures at BiSe/Cu interface change dramatically as one moves (even slightly) below or above . Comparing Figs. 3(a)–(d) with 3(e)–(h) makes it possible to understand the effect of the magnetization of the Co layer, which modifies Nagaosa2010 () Rashba dispersion around and the corresponding spin textures (particularly the out-of-plane component).

The theoretical modeling of SOT in TI/FM Mellnik2014 (); Mahfouzi2016 () or heavy-metal/FM Haney2013 () bilayers is usually conducted by starting from strictly 2D Hamiltonians, such as or the Rashba ferromagnetic model Nagaosa2010 (), respectively, so that the FM layer is not considered explicitly. Figures 4(e)–(h) show that this is not warranted since the BiSe layer induces proximity SOC and the corresponding in-plane spin texture on the first ML of Co, which decays to zero only after reaching plane in Fig. 1(b). In fact, we find non-trivial in-plane spin texture even on the surface of Co in contact with vacuum, as shown in Figs. 4(b)–(d), which is nevertheless quite different from those in Figs. 4(f)–(h). The in-plane spin texture in Figs. 4(b)–(d) is a consequence of the Rashba SOC enabled by inversion asymmetry due to Co surface Chantis2007 () where an electrostatic potential gradient can be created by the charge distribution at the metal/vacuum interface and thereby confine wave functions into a Rashba spin-split quasi-2D electron gas Bahramy2012 (). Thus, Figs. 4(a)–(d) explains the origin of recently observed Emori2016 () SOT in the absence of any adjacent heavy metal or TI layer since passing current parallel to MLs of Co hosting nonzero in-plane spin textures will generate a nonequilibrium spin density Chang2015 () and spin-orbit torque Mahfouzi2016 (); Haney2013 ().

Finally, we propose a purely charge transport measurement that could detect which among the spin-textures shown in Figs. 3(b)–(d) resides at the Fermi level of TI/FM interface. Our scheme requires to fabricate vertical heterostructure in Fig. 1(b) and measure its tunneling anisotropic magnetoresistance (TAMR). The TAMR is a phenomenon observed in magnetic tunnel junctions with a single FM layer Mahfouzi2014a (); Park2008 (); Chantis2007 (); Matos-Abiague2009 (), where SOC makes the band structure anisotropic so that the resistance of such junctions changes as the magnetization is rotated by angle or in Fig. 1(b). The resistance change is quantified by the TAMR ratio defined as Chantis2007 (); Matos-Abiague2009 ()

(5) |

Here for where magnetization in Fig. 1(b) rotates in the plane perpendicular to the TI/FM interface, and for where magnetization in Fig. 1(b) rotates within the plane of the TI/FM interface. In the case of , is the resistance when in Fig. 1; and in the case of , is the resistance when in Fig. 1. Thus, changes due to the different orientations of the magnetization with respect to the direction of the current flow, while the situation becomes more subtle for where the magnetization remains always perpendicular to the current flow. Figure 5(b) demonstrates that the largest is obtained by tuning the Fermi level to eV so that nearly helical spin texture in Fig. 3(d) resides at the Fermi level. Another signature of its presence is rapid increase of when tilting by small angles away from the current direction. The in-plane shown in the inset of Fig. 5(b) is much smaller (and difficult to converge in the number of transverse -points) quantity which does not differentiate between spin textures shown in Figs. 3(b)–(d).

## I Methods

We employed the interface builder in the VNL vnl () and CellMatch Lazic2015 () packages to construct a common unit cells for: (a) BiSe/Cu(111) bilayer, where the common unit cell is 5 5 in size compared to the smallest possible Cu(111) slab cell and copper is under compressive strain of 0.9 % while BiSe lattice constant is unchanged; (b) BiSe/Co(0001) bilayer where Co(0001) has the same lattice constant as BiSe, so the same unit cell as for Cu(111) is used without any strain on Co(0001). These two unit cells are illustrated in Figs. 6(a) and 6(b), respectively. In order to determine the best stacking of atomic layers and the distance of BiSe atoms with respect to surfaces of Cu(111) and Co(0001), we use DFT calculations as implemented in the VASP package Kresse1993 (). The electron core interactions are described by the projector augmented wave (PAW) method Blochl1994 (), and vdW-DF Dion2004 () with optB88 is used as density functional Mittendorfer2011 () in order to describe van der Waals (vdW) forces between QLs of BiSe or between BiSe and metallic layers. The cutoff energy for the plane wave basis set is 520 eV for all calculations, while -points were sampled at 33 surface mesh. We use Cu and Co layers consisting of 5 MLs, where 3 bottom MLs are fixed at bulk positions while the top two metallic MLs closest to BiSe are allowed to fully relax until forces on atoms drop below meV/Å. In order to avoid interaction with periodic images of the bilayer, Å of vacuum was added in the –direction.

For the case of BiSe on Co(0001), the most favorable position yields a binding energy of 460 meV per Co atom. Both ML of Co and QL of BiSe in direct contact gain some corrugation, roughly around Å, while the average –distance between them is 2.15 Å. The average distance between the ML of Cu and QL of BiSe in direct contact is around 2.26 Å with smaller corrugation than in the case of Co(0001), while the binding energy is 294 meV per Cu atom. For other relative positions of BiSe layer with respect to Cu(111) and Co(0001) layers the difference in binding energy is very small. Binding energies in both cases are rather small, thereby signaling the dominant vdW forces. Nevertheless, some charge rearrangement does occur at the interface due to push back/pillow effect Vazquez2007 (), as shown in Figs. 6(c) and 6(d) where charge rearrangement is more pronounced for the case of BiSe/Cu(111) interface.

The calculation of the retarded GF in Eqs. (1) and (3) requires represented in the linear combination of atomic orbitals (LCAO) basis set which makes it possible to spatially separate system into the active region attached to one or two semi-infinite leads, as illustrated in Figs. 1(a) and 1(b), respectively. We employ ATK package atk () for pseudopotential-based LCAO noncollinear DFT calculations yielding , from which we obtain retarded GFs and the corresponding spectral functions, as well as the resistance in Eq. (4). In ATK calculations, we use Perdew-Burke-Ernzerhof (PBE) parametrization of generalized gradient approximation for the exchange-correlation functional; norm-conserving pseudopotentials for describing electron-core interactions; and LCAO basis set generated by the OpenMX package Ozaki2003 (); openmx () which consists of s2p2d1 orbitals on Co, Cu and Se atoms, and s2p2d2 on Bi atoms. These pseudoatomic orbitals were generated by a confinement scheme Ozaki2003 () with the cutoff radius 7.0 and 8.0 a.u. for Se and Bi atoms, respectively, and 6.0 a.u. for Co and Cu atoms. The energy mesh cutoff for the real-space grid is chosen as 75.0 Hartree.

## Ii Acknowledgments

We thank Takafumi Sato and Jia Zhang for insightful discussions. J. M. M.-T., K. D. and B. K. N. were supported by NSF Grant No. ECCS 1509094. J. M. M.-T. also acknowledges support from Colciencias (Departamento Administrativo de Ciencia, Tecnologia e Innovacion) of Colombia. P. L. was supported by the Unity Through Knowledge Fund, Contract No. 22/15 and H2020 CSA Twinning Project No. 692194, RBI-T-WINNING. S. Smidstrup, D. Stradi and K. Stokbro acknowledge support from the European Commission Seventh Framework Programme Grant Agreement IIIâV-MOS, Project No. 61932; and Horizon 2020 research and innovation programme under grant agreement SPICE, project No 713481. The supercomputing time was provided by XSEDE, which is supported by NSF Grant No. ACI-1053575.

## References

- (1) A. R. Mellnik et al., Nature 511, 449 (2014); Y. Fan et al., Nat. Mater. 13, 699 (2014); Y. Wang et al., Phys. Rev. Lett. 114, 257202 (2015).
- (2) Y. Fan et al., Nat Nano 11, 352 (2016).
- (3) Y. Shiomi et al., Phys. Rev. Lett. 113, 196601 (2014); P. Deorani et al., Phys. Rev. B 90, 094403 (2014); M. Jamali et al., Nano Lett. 15, 7126 (2015); J.-C. Rojas-Sánchez et al., Phys. Rev. Lett. 116, 096602 (2016); H. Wang et al., Phys. Rev. Lett. 117, 076601 (2016).
- (4) K. Kondou et al., Nat. Phys. 12, 1027 (2016).
- (5) D. Pesin and A. H. MacDonald, Nature Mater. 11, 409 (2012).
- (6) P.-H. Chang et al., Phys. Rev. B 92, 201406(R) (2015).
- (7) F. Mahfouzi, N. Nagaosa, and B. K. Nikolić, Phys. Rev. B 90, 115432 (2014).
- (8) F. Mahfouzi, B. K. Nikolić, and N. Kioussis, Phys. Rev. B 93, 115419 (2016).
- (9) A. Bansil, H. Lin, and T. Das, Rev. Mod. Phys. 88, 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, 1057 (2011).
- (11) A. Soumyanarayanan, N. Reyren, A. Fert, and C. Panagopoulos, Nature 539, 509 (2016).
- (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).
- (14) E. Zhao, C. Zhang, and M. Lababidi, Phys. Rev. B 82, 205331 (2010).
- (15) H. Zhang et al., Nat. Phys. 5, 438 (2009).
- (16) A. A. Soluyanov et al., Nature 527, 495 (2015).
- (17) X.-L. Sheng and B. K. Nikolić, arXiv:1610.02719 (2016).
- (18) T. Shoman et al., Nat. Commun. 6, 6547 (2015).
- (19) N. Marzari et al., Rev. Mod. Phys. 84, 1419 (2012).
- (20) C. D. Spataru and F. Léonard, Phys. Rev. B 90, 085115 (2014).
- (21) J. Zhang, J. P. Velev, X. Dang, and E. Y. Tsymbal, Phys. Rev. B 94, 014435 (2016).
- (22) T. Ozaki, Phys. Rev. B 67, 155108 (2003).
- (23) J. Velev and W. Butler, J. Phys.: Condens. Matter 16, R637 (2004); I. Rungger and S. Sanvito, Phys. Rev. B 78, 035407 (2008); H. H. B. Sørensen et al., Phys. Rev. B 77, 155301 (2008).
- (24) J. Kim et al., Nat. Mater. 12, 240 (2013).
- (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. van Schilfgaarde, Phys. Rev. Lett. 98, 046601 (2007).
- (28) M. S. Bahramy et al., Nat. Commun. 3, 1159 (2012).
- (29) S. Emori et al., Phys. Rev. B 93, 180402 (2016).
- (30) B. G. Park et al., Phys. Rev. Lett. 100, 087204 (2008).
- (31) A. Matos-Abiague and J. Fabian, Phys. Rev. B 79, 155303 (2009).
- (32) P. Lazić, Comp. Phys. Commun. 197, 324 (2015).
- (33) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993); G. Kresse and J. Furthmüllerb, Comput. Mater. Sci. 6, 15 (1996).
- (34) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994); G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
- (35) M. Dion et al., Phys. Rev. Lett. 92, 246401 (2004).
- (36) F. Mittendorfer et al., Phys. Rev. B 84, 201401 (2011).
- (37) H. Vázquez, Y. J. Dappe, J. Ortega, and F. Flores, J. Chem. Phys. 126, 144703 (2007).
- (38) Atomistix ToolKit (ATK) 2016.3, http://www.quantumwise.com.
- (39) Virtual NanoLab (VNL) 2016.3, http://www.quantumwise.com.
- (40) OpenMX 3.8, http://www.openmx-square.org/.