Quantitative comparison of electrically induced spin and orbital polarizations in heavy-metal/3d-metal bilayers

Electrical control of magnetization is of crucial importance for integrated spintronics devices. Spin-orbit torques (SOT) in heavy-metal/ferromagnetic heterostructures have emerged as promising tool to achieve efficiently current-induced magnetization reversal. However, the microscopic origin of the SOT is being debated,with the spin Hall effect (SHE) due to nonlocal spin currents and the spin Rashba-Edelstein effect (SREE) due to local spin polarization at the interface being the primary candidates. We investigate the electrically induced out-of-equilibrium spin and orbital polarizations in pure Pt films and in Pt/3d-metal (Co, Ni, Cu) bilayer films using ab initio electronic structure methods and linear-response theory. We compute atom-resolved response quantities that allow us to identify the induced spin-polarization contributions that lead to fieldlike SOTs, mostly associated with the SREE, and dampinglike (DL) SOTs, mostly associated with the SHE, and compare their relative magnitude, dependence on the magnetization direction, as well as their Pt-layer thickness dependence. We find that both the FL and DL components contribute to the resulting SOT at the Pt/Co and Pt/Ni interfaces, with the former contributions being larger at the Pt interface layer and the latter larger in the Co or Ni layers. Our calculations show that the electrically-induced transverse orbital polarization is exceedingly larger than the induced spin polarization and present even without spin-orbit coupling, in contrast to the spin polarization.


I. INTRODUCTION
Electrical control of magnetization has attracted considerable attention because of its potential for high-speed spin-based memories with low-power consumption. Following theoretical predictions [1,2] it was shown that the magnetization of a ferromagnetic layer in a multilayer stack can be switched with a spin-transfer torque (STT) exerted by a spin-polarized electric current flowing through the magnetic layer in perpendicular direction [3][4][5][6][7]. STT enabled the development of current-operated nonvolatile spin-logic devices, such as the STT-magnetic random access memory (STT-MRAM) [8]. While STTbased technology is a step forward, there are still shortcomings, such as unintended switching that can occur as the write and read currents flow in the same direction [9].
A different concept to electrical magnetization switching is the more recently discovered spin-orbit torque (SOT) [10][11][12][13]. SOT can be observed in a heavymetal/ferromagnetic bilayer film where the current flows dominantly through the heavy metal and parallel to the ferromagnetic layer. In this configuration, reversible magnetization switching can be achieved in a very energy efficient way and, moreover, have read and write currents flow in distinct directions through the device [14][15][16][17].
While it is evident from experiments that the SOT can be used to efficiently reverse the magnetization in the magnetic layer, its microscopic origin is still to be fully understood. Two candidates for driving the SOT have attracted much attention: the spin Hall effect (SHE) [18,19] and the spin Rashba-Edelstein effect (SREE) [20]. * leandro.salemi@physics.uu.se Both effects are caused by the spin-orbit interaction, either in the bulk of the material or at an interface, yet their microscopic appearance is drastically different. The SHE is a nonlocal effect wherein an electrical current generates the flow of a transverse spin current to the boundary of the conducting slab (see [21][22][23][24]) where it exerts a torque on the adjacent ferromagnetic layer. The SREE conversely is a local effect: as pointed out by Edelstein [20], a nonequilibrium spin polarization is generated at a symmetry-broken interface by an electric current in the presence of Rashba spin-orbit coupling (SOC) [25]. Both effects have been discussed in the context of SOT switching, in some cases the SHE was considered as the dominant effect [12,13] whereas in other cases the focus was on the SREE [10,11,14,26]. In heavy-metal/ferromagnetic bilayer structures both effects are expected to be present simultaneously and will contribute to the fieldlike (FL) SOT and dampinglike (DL) SOT [27][28][29][30][31][32][33], yet their relative contribution remains disputed and continues to be a topic of contemporary investigations [34][35][36] (see also [37] for a recent review).
The SHE and SREE are however not the only magnetic effects that can occur. It was discovered theoretically that, in addition to the spin polarization induced by a current through the SHE, also a nonequilibrium orbital polarization can be induced, which represents an orbital Hall effect (OHE) [41][42][43][44][45]. Similarly, the presence of spatial symmetry breaking in a material was recently shown to lead to a significant local orbital polarization, i.e., an orbital Rashba-Edelstein effect (OREE) [46]. Both the OHE and OREE are currently only poorly understood, in terms of their relative magnitudes as well as directions of arXiv:2004.11942v5 [cond-mat.mtrl-sci] 14 Jun 2021 the induced orbital torques. So far several first-principles calculations have been reported for the OHE [41][42][43][44][45]. A direct observation of the induced orbital polarization is yet to be achieved in experiments (see Refs. [47,48] for recent studies).
In this work, we employ relativistic density functional theory (DFT) and Kubo linear-response theory to compute the spin and orbital response to an external electric field for realistic metallic bilayer structures in which Pt is chosen as the heavy-metal material. Specifically, four different systems are investigated: a pure Pt system and three Pt/3d-metal bilayer systems, where the 3d element is Ni, Co or Cu. For these we compute the spin and orbital conductivity and magneto-electric (ME) tensors resolved for the individual atomic layers in the metallic heterostructures.
In the following, we first introduce the theoretical framework of linear response within DFT and subsequently apply our formalism to compute the spin and orbital responses for the considered bilayer systems, for various Pt thicknesses. We analyze the spatial symmetry of the spin response, which is embodied in the spin ME susceptibility tensor χ s . We show that the direction of the induced spin magnetization δS depends on the relative directions of the applied electric field E, the equilibrium magnetization direction M , and the system geometry. The tensors can be decomposed into odd-in-M and even-in-M components. Based on this parity with respect to M , we discuss their relationship to the DL and FL components of the SOT.
The relative importance of those tensor contributions strongly depends on the position of the atomic layer in the slab, and, to a lesser extent, to the thickness of the Pt slab. We investigate furthermore the magnetizationdirection dependence of the spin and orbital responses. Whereas the spin response vanishes when the spin-orbit interaction is turned off, off-diagonal components of the orbital ME tensor remain. We compute effective spinorbit torques on the magnetic Ni and Co layers and compare our results with previously reported values. Finally, we compare the relative importance of the FL and DL components of χ S and discuss their relationship to the SHE ad SREE. We find that the induced spin polarization on the Pt side of the Pt/Ni and Pt/Co interfaces mainly leads to a FL SOT while the induced spin polarization further into the FM layer contributes more to the DL torque component.

A. Linear response
The materials are modeled within DFT by the relativistic Kohn-Sham Hamiltonian as implemented in WIEN2k [49],Ĥ whereĤ 0 is the relativistic Kohn-Sham Hamiltonian, |nk the single-electron Kohn-Sham state for band n at wavevector k and nk the corresponding eigenenergy. Under the influence of an external perturbationV = −er·E where e is the electron charge, E the external electric field andr the position operator, the change δA in expectation value of a vectorial observable A associated to vector operatorÂ, can be expressed within the linear-response formalism [27,29,30,50] as The response χ A ij is expressed in terms of solutions of H 0 , with m e the mass of the electron, f nk the occupation of Kohn-Sham state |nk , Ω the Brillouin-zone volume, p j nmk thep j momentum-operator matrix element, A i mnk theÂ i -operator matrix element and ω nmk = nk − mk , the difference of Kohn-Sham eigenergies. As discussed below, we use forÂ the spin and orbital angular momentum operators,Ŝ andL, as well as the the spin and orbital current-density operators,Ĵ S andĴ L . The product J iS mnk p j nmk is proportional to the spin Berry connection [48]. The quantity τ inter (τ intra ) is the electronic lifetime for inter (intra) band transitions. In the main part of this work, τ inter and τ intra are set to τ −1 inter = 0.272 eV and τ −1 intra = 0.220 eV. Those values have been determined by comparing linear-response calculations to experimental conductivity data for Pt thin films [51]. The influence of the electronic lifetime on the computed nonequilibrium spin and orbital polarizations is discussed in Sec. III A 4.

B. Angular momentum and flow of angular momentum
The induced angular momentum is composed of a spin and orbital contribution. Let us first focus on the spin part.
The spin operatorŜ and spin-density current operator J S k can be defined aŝ whereσ x ,σ y , andσ z are the Pauli matrices, {. , .} denotes the anti-commutator, V is a reference volume and k (k = x, y, z) an index specifying the direction of the spin polarization carried by the spin-current density. In this work, V refers to the individual atomic spheres, allowing us to compute atom-projected quantities (see Appendix A for details).
Using the linear-response formalism, we can compute the out-of-equilibrium electrically induced spin angular momentum δS as well as the induced spin-current density J S k , using where χ S is the spin ME susceptibility tensor and σ S k the spin conductivity tensor. Both χ S and σ S k are real 2 nd -rank tensors, but note that, due to the spin component dependence, the spin conductivity tensor can be associated with a 3 rd -rank tensor, σ S . Analogous quantities can be straightforwardly defined for the orbital angular momentumL. Thus, we can define the orbital ME susceptibility tensor χ L and orbital conductivity tensor σ L k , where δL is the out-of-equilibrium electrically induced orbital angular momentum and J L k the induced orbital current density.

C. Computational methodology
The bilayer structures that are studied here consist of several Pt monoatomic layers that are covered with two monoatomic layers of the 3d elements Ni, Co or Cu (see Fig. 1). For comparison, we also study the pure Pt system, where the top two monolayers consist of Pt. The nomenclature used in this paper is the following: we denote our systems by nPt/2Y where n is the total number of Pt monolayers and Y is either Ni, Co, Cu or Pt. The minimum total number of Pt monolayers used in our calculations is 2 while the maximum is 18 (denoted as 16Pt/2Pt). The maximum thickness achieved is then ∼3.2 nm. The direction normal to the interfaces is taken as the z axis.
The monoatomic layers are labeled from z = 1 for the Pt monoatomic layer at the interface with vacuum (leftmost layer in Fig. 1) to z = n + 2 for the Y monoatomic layer at the interface with vacuum (rightmost layer in Fig. 1). Particular positions can be identified, like z = n for the Pt monoatomic layer at the Pt/Y interface and z = n + 1 for the Y monoatomic layer at the Pt/Y interface.
To compute the spin and orbital susceptibility and conductivity tensors, we use the following three-step procedure. As the DFT packages used employ full 3D periodic boundary conditions, all heterostructures contain 20Å of vacuum to avoid spurious interactions with neighboring simulation cells. More details on the computational recipe are given in Appendix A.

Symmetry of the ME tensors
The structures studied in this work are tetragonal, with nonmagnetic point group 4mm. The response tensor χ S can be expanded as where M k is the k-th component of the unit vector of magnetization, i.e., M = |M |M, with M the equilibrium magnetization. As detailed in Ref. [53], for M u z , χ S can be written as where χ S xy = −χ S yx and χ S xx = χ S yy = χ S zz . With M out-of-plane, the system exhibits an in-plane x/y spatial symmetry, which is fully recovered in our calculations.
The χ S tensor can be further decomposed into an oddin-M and even-in-M component, with specifically, A nonzero odd-in-M part can obviously not exist for nonmagnetic systems (nPt/2Pt and nPt/2Cu), which is as well recovered in our calculations. The spin response is highly dependent on the magnetization direction. Setting the magnetization in plane, M || u x , the χ S tensor can be written as being different from the M || u z case. Now the χ xy , χ yx elements are even-in-M and the χ xz , χ zx elements oddin-M . We can furthermore mention already that the orbital χ L tensor has the same nonzero elements with the same M parity. Note that the magnetization direction breaks the x/y symmetry such that the xy and yx components of χ S/L may be different. Such difference, if relevant, comes from terms with powers of M of order higher than two in Eq. (10). As we will see later, those higher order corrections are indeed relevant for the spin response. The symmetry of the system is naturally incorporated in the electronic density n(r) and magnetization density m(r) computed with DFT. As a result, our calculations, in which we do not impose any symmetry constraints, provide the full M -dependent response, i.e., the left hand-side of Eq. (10). Depending on the relative orientation of the induced spin polarization δS with respect to (1) the applied electric field E, (2) the normal direction u z and (3) the equilibrium magnetization vector M , the components of the ME susceptibility χ S can be classified according to three categories: • M -transverse components (M ⊥ ): • M -longitudinal component (M ):

Symmetry of the SOT
The SOT T felt by the equilibrium magnetization M can be written as where δB is the SOT effective magnetic field, which originates from the electrically-induced change in the exchange-correlation effective field B XC . We can approximate δB as where S is the equilibrium value of the expectation value ofŜ. Essentially, Eq. (19) assumes that the direction of the exchange-correlation effective field changes while its strength is unaffected. Note that because the spatial distribution of δS may differ for different orbital characters (e.g., s, p, d), Eq. (19) is an estimate of δB.
Using M = −2µ B S, where µ B is the Bohr magneton, the torque can be expressed as where we used that M = S/|S|. With the induced spin polarization [Eq. (6)] this can be written as Since χ S can be decomposed in odd-in-M and even-in-M components, we can always decompose the SOT in an odd (o) and an even (e) component, We can compute T o and T e , e.g. for an in-plane field, E = |E| u x . When M ||u z , using Eq. (13) we obtain while for M || u x this yields [see Eq. (14)] We can now easily establish the link between the symmetry of the components of χ S , their evenness or oddness in M , and the commonly discussed fieldlike torque T FL and dampinglike torque T DL , specifically, We recognize that T o , and therefore the E ⊥ components, correspond to T FL while T e , and therefore the M ⊥ components, correspond to T DL . A similar classification can be carried out for the cases where M || u x and M || u y . The classification and symmetry relations for these magnetization directions are summarized for convenience in Table I. The spatial symmetry of the induced spin polarization δS and whether it leads to a dampinglike or fieldlike spin-orbit torque T are also provided.

III. RESULTS
A. Spin response

Magnetization out-of-plane
We start with the case where M || u z . It is instructive to consider first the thickest heterostructures, i.e., 16Pt/2Ni, 16Pt/2Co, 16Pt/2Cu, and 16Pt/2Pt. In Figs In all cases, we observe that the response of the Pt atomic-layer at the vacuum interface (z = 1) is virtually independent on the type of Y atom used, suggesting that these systems are thick enough to isolate the Pt/3d-interface properties. The inclusion of the two monoatomic layers of 3d elements mainly impacts the χ S and σ S profiles close to their interface.
For the E-transverse components (Figs. 2(a) and (d)), both χ S and σ S are qualitatively barely impacted in the bulk Pt region by the replacement of the two last Pt atomic monolayers by two 3d atomic monolayers. The  profile of σ Sy zx is in all cases mostly defined by a plateau in the center of the Pt layer. This specific component of σ S is customarily identified with the SHE in Pt-bulk calculations. The spin-accumulation profile related to χ S yx across the Pt layer strongly resembles the type spin accumulation that is expected from transverse spin flow. Reversing the magnetization of the Ni and Co layers from +u z to −u z in the calculations, moreover, does not have a notable effect on the (M -even) spin-accumulation given by χ S yx . These results are in agreement with what is observed in SHE calculations for bulk Pt [51,54]. Note that the accumulated spin moment is given by δM = For the 16Pt/2Cu system, the discussion cannot be performed on the basis of torques (M is ill-defined). For the spin accumulation related to χ S yx , compared to the pure Pt case, we observe an increase of accumulated spin density at the Pt side of the Pt/Cu interface as well as a decrease of spin density within the Cu layer.
The M ⊥ components (Figs. 2(b) and (e)), are nonexistent for non-magnetic pure Pt and Pt/Cu system. Siz-able values are only obtained close to the interface with Co and Ni, both for χ S and σ S . Remarkably, while the E ⊥ and M ⊥ components are comparable in size close to the interface, their features differ greatly: (1) there is no bulk-like behavior for σ Sy zy /σ Sx zx , (2) those components are non-existent in bulk Pt, and (3) they are magnetization direction dependent. Although there is a spatial symmetry breaking at the Pt/Cu interface, no spin polarization is induced, which strongly suggests the importance of spin-split electronic states at the interface for χ S xx and σ Sx zx . The odd-in-M spin conductivity σ Sx zx was recently named magnetic SHE [60], to distinguish it from the conventional SHE σ Sy zx that is even-in-M [61].
As aforementioned, the M -longitudinal components (see Figs. 2(c) and (f)) are peculiar in the sense that they are not SOT-related (the usual SOT configuration does not involve out-of-plane electrical fields and, also, no torque is generated by a spin accumulation parallel to the static moment). Although such components can in principle be obtained via a symmetry analysis (cf. Ref. [53]), they haven't been, to the best of our knowledge, investigated so far. As will be clarified further below, this effect is due to the spin-orbit interaction. Here, an electric field applied parallel to the out-of-plane magnetization causes a sign-changing spin polarization along M in the ∼5 topmost monolayers. This is clearly a magnetic effect, as it does not exist for the nonmagnetic systems. The spin conductivity σ Sz zz shows a decaying behavior from z = 16 to z = 1, but this decay is slower than that of the equilibrium spin magnetization in the systems. Also, similar to the M -transverse component, no "bulk-like" behavior is observed. A possible way of observing this previously unidentified SOC-induced effect could be achieved by gating the ferromagnetic layer from the top and monitor a change of its magnetization.
So far, we focused on the components of the spin conductivity tensor giving rise to spin currents flowing along u z . While those components are the ones that should be of interest for understanding SOT in bilayer structures, other nonzero components can be observed, as well. For a magnetic system with the 4mm point group, with M u z , we find that σ S can be written as The components associated to the SHE, i.e., σ S k ij , where the indices are such that ijk = 0 ( ijk is the Levi-Civita symbol), are nonzero whether the system is magnetic or not. However, while in cubic systems like bulk Pt they are all equal in magnitude, here, because of the symmetry breaking with respect to the z axis, the tensor elements are not invariant under exchange of z and x or y indices. The other 7 components, that is σ Sx xz , σ Sx zx , σ Sy yz , σ Sy zy , σ Sz xx , σ Sz yy , and σ Sz zz only exists for magnetic systems. Some of those magnetic components have been discussed recently [61][62][63][64][65].

Pt-thickness dependence
Next, we investigate the Pt layer thickness dependence of the E ⊥ , M ⊥ , and M components of χ S . The number of Pt monolayers for our nPt/2Y systems is varied from n = 2 (Pt thickness ∼ 0.38 nm) to n = 16 (Pt thickness ∼ 3.08 nm). Figure 3 shows the computed Pt-thickness dependence where each column of the figure focuses on one particular atomic monolayer, with, from left to right, the Pt monolayer at the Pt/Y interface, the Y monolayer at the Pt/Y interface, and the Y atomic monolayer at the Y /vacuum interface. Each row focuses on one particular component, namely, from top to bottom, the E ⊥ , M ⊥ , and M components. The values of the tensor elements that give rise to the SOT, the E ⊥ and M ⊥ components, barely fluctuate beyond n = 8 (Pt thickness ≥ 1.54 nm). Thus, both T FL and T DL approach their maximum values already for relatively thin Pt layers, consistent with recent investigations [33,66].
For the E ⊥ components, in the case of pure Pt (nPt/2Pt), χ S yx tends to increase the closer we come to the last layer, which is typically what we would expect from a transport-generated spin accumulation profile. When the two last layers are replaced by a magnetic element (Y = Co or Ni) drastic changes occur. First, we observe that χ S yx is bigger for the 3d monolayer closer to the Pt layer than for the second Y layer. Second, at a fixed position, χ S yx is bigger for Y = Ni than Y = Co. For the M ⊥ components, associated to T DL , the χ S xx for the Pt monolayer at the Pt/Y interface ( Fig. 3(d)) is virtually identical for Y = Ni or Co and for all Pt thicknesses considered. In the first and second magnetic monolayer (Figs. 3(e)) and (f)) the χ S xx is, in both cases, bigger for Y = Ni than Y = Co. The bottom row, lastly, shows the M -longitudinal spin accumulation. Also here, we obtain that a larger magnitude of χ S zz is generated for Y = Ni.

Magnetization-direction dependence of the spin response
The direction of the magnetization vector M is often rotated in out-of-plane or in-plane directions in experiments [35,67]. The magnetization-direction unit vector M can be written as M = (sin θ cos φ, sin θ sin φ, cos θ), where φ and θ are the azimuthal and polar angles, respectively. To capture the M-dependence of χ S , we selfconsistently compute χ S at different (θ, φ) values. We choose the 12Pt/2Co system for these simulations and consider the ME spin response on the first Co atom at the interface with the Pt layer. In particular, we compute χ S (θ, φ) for M constrained in the x − y plane (Fig. 4) and M constrained in the x − z plane (Fig. 5).
In Fig. 4 we show the computed angular dependence of the tensor components when the magnetization is inplane and rotated from angle φ = 0 • (i.e., parallel to u x ) to 360 • . The panels reveal that spin ME tensor is, to great precision, given by the following parametric form, where α i (i = 1, 2, · · · , 5) are constants in units of 10 −3 nm V , whose values are α 1 = 12.54, α 2 = −20.83, α 3 = 11.50, α 4 = 77.18, and α 5 = −8.36. Note that this functional form only gives the φ-dependence, that is, θ is kept constant to 90 • . The computed angle dependence illustrates that the tensor elements, and thus the resulting SOTs, obey a distinct angular symmetry.
Terms whose angular dependency varies as 2φ are even under magnetization reversal, while those that vary as φ are odd. Let us consider for example the case where M is along u x , that is φ = 0. If we apply a field E along M , the only relevant nonzero elements, that is, the elements that give rise to a torque, are χ S yx and χ S zx . The former tensor element is even in the magnetization, giving consistently a '−α 2 − α 3 cos 2φ' angular dependence, while the latter element is odd in M , giving a cos φ dependence. Note that the expression (21) for the SOT contains the cross product M × (χ S E), which implies that an additional trigonometric function (sin φ, cos φ) appears in the SOT. Analogous calculations have been carried out for M being varied along the polar angle θ, with φ set to 0. Figure 5 shows the computed θ-angle dependence of the ME tensor χ S for the 12Pt/2Co system. It can be recognized that the angular variation of the matrix elements obeys the θ-dependence where β i and β 2 are constants. Here, we only give the generic trigonometric dependence of the tensor elements.
In this case there are nine constants needed to describe the angle dependence precisely as in Eq. (27). The angular dependence of the χ S tensor was computed for a particular atomic position in the bilayer. Performing further calculations, we find that the obtained angular dependence [Eqs. (27) and (28)] is valid at any atomic position in the layer, albeit with vanishing constants for the odd-in-M components when the atomic layer is close to the vacuum/Pt interface. For each atomic position there are different constants for the tensor components, however, a decisive property is, whether the tensor component is E-transverse, M -transverse, or Mlongitudinal. As illustrated in Table I, in this way equivalent tensor elements can be identified for different M directions.
In Fig. 6 we compare the atomic-layer resolved profiles of χ S for the 6Pt/2Ni system, for M || u x and M || u x . We compare the E ⊥ , M ⊥ and M components, which allows us to track the physically equivalent quantities with respect to the SOT symmetry (FL or DL). We use in Fig.  6 the superscript u z (u x ) to denote quantities computed with M || u z (M || u x ). For most of the atomic positions, the corresponding tensor elements for both M || u z and M || u x show a strikingly similar behavior, with the exception of the χ S xy element that differs most in the Ni layer. The similar behavior and approximately identical values illustrate that the tensor elements transform in an almost rigid manner when the magnetization vector is rotated. The larger difference in the χ S xy component suggest that the orientation of M has a stronger influence on the Pt/Ni interface electronic structure. Note that although χ S,ux xy and χ S,uz xy differ significantly at the Pt/Ni interface (see Fig. 6(a)), χ S,ux xy does not become effective, as it would give rise to an induced magnetization along the equilibrium magnetization, and can therefore not contribute to SOT.
For the M ⊥ components, shown in Fig. 6(c), very similar magnitudes are observed but there is an opposite sign in the χ S,uz xx and χ S,ux zx components. This sign reversal is easily recovered using our proposed classification. Indeed, in the case M || u z , one finds from Eq. (16) for the corresponding component while for M || u x we have, which perfectly captures the sign reversal. The Mlongitudinal χ S components, lastly, follow the approximate transformation behavior quite well, see Fig. 6(d). These calculations illustrate the usefulness of the classification of the tensor elements according E ⊥ , M ⊥ , and M || for identifying their resulting torque properties.

Electronic lifetime dependence
In the above calculations we employed the intraband and interband electron lifetimes ( τ intra = 0.220 eV, τ inter = 0.272 eV) obtained from fitting the calculated linear-response conductivity to experimental electrical conductivities of Pt films [51]. In general, these electronic lifetimes depend on the sample purity and microstructure. To investigate the influence of these lifetimes on the nonequilibrium spin and orbital polarizations we have varied these parameters in our calculations, adopting the 12Pt/2Co and 14Pt films as typical model systems. We consider the case M || u z for the 12Pt/2Co film, for which the nonzero spin tensor elements are χ S yx , χ S xx , and χ S zz . For the 14Pt film only the tensor element χ S yx is nonzero, cf. the results shown in Fig. 2. At this point it is relevant to mention that the M -transverse tensor element χ S xx has only interband contributions, whereas the E-transverse spin polarization χ S yx has both intra-and inter -band contributions. The calculated dependence of the torque-related, spin ME responses on the electronic broadening δ (= /τ ) is shown in Fig. 7. The intraband part of χ S yx increases linearly for decreasing lifetime broadening, as expected from Eq.
(3). The interband contributions, conversely, vary with the lifetime broadening but become more broadening independent for large values of δ. For small δ the interband contributions converge to nonzero values (but the values for χ S yx become small). For realistic lifetime values ( /τ ∼ 0.25 eV) the intra-and interband induced spin polarizations have the same magnitude and need to be treated on equal footing. Note that semiclassical Boltzmann transport theory only captures the intraband spin ME component [68], an approximation that is viable for very pure crystals with long electron lifetimes.
In the limit of very small electronic broadenings, the intraband contribution to χ S yx will dominate completely the electrically induced spin polarization. This is the component that leads to a fieldlike SOT, see Table I. The intraband part of the spin response is commonly associated with the Rashba-Edelstein effect [53,60], consistent with the shape of the Rashba spin-orbit coupling Hamiltonian [25] that leads to a fieldlike torque in Landau-Lifshitz-Gilbert equations of spin dynamics. The spin Hall effect is conversely commonly associated with an interband contribution [24]. Such contribution is present both in χ S yx and χ S xx , i.e., in the E ⊥ and M ⊥ terms. The SHE will thus contribute to both a fieldlike and a dampinglike SOT. The SHE is considered to lead to a DL torque [12] when using the Slonczewski model [1]. The SHE was however also proposed to lead to a FL torque [28,69], consistent with our calculations. The latter identification was made within the microelectronic cir- cuit model, wherein the M ⊥ (E ⊥ ) component is mainly due to the SHE (SREE), when the spin mixing conductance is chiefly real [28,69]. To obtain such result, it is assumed that the transverse spin current and spin accumulation exists in the HM layer only, i.e., these quantities are zero in the FM layer. This might be a reasonable approximation for thicker layers, but in our case, where we compute atomistic quantities, we find a nonzero spin conductivity and spin accumulation (see Fig. 2) in the FM layer which is only two atoms thick. For a very small lifetime broadening the FL E ⊥ contribution will dominantly stem from the intraband SREE and the DL M ⊥ contribution from the SHE. We note however with respect to this discussion, that our DFT Hamiltonian contains the full form of the spin-orbit interaction and is thus different from the more elementary Bychkov-Rashba SOC [25], but it provides all materials' specific SOC effects.

Atomic layer-resolved orbital response
A similar analysis can be performed for the orbital response, both in terms of χ L and σ L . While similarities are observed, unique characteristic can be observed, too. To start with, we show in Fig. 8 the calculated layer-resolved orbital ME susceptibilities χ L and orbital conductivities σ L for the 16Pt/2Y systems, for M || u z , similar to the spin counterparts shown in Fig. 2. For the sake of completeness, we provide in Appendix B analogous plots to Fig. 3 for the Pt-thickness dependence and to Figs. 4, 5, and 6 for the angular dependence of the χ L tensors and their transformation properties under rotation of the magnetization direction. Also the dependence of the χ L tensor elements on the lifetime broadening is given in Appendix B.
The layer-resolved results, shown in Fig. 8(a), reveal that, just like for the spin, the E-transverse component in the pure Pt system resembles strongly the transport-induced accumulation of orbital angular momentum. The transverse conductivity σ Ly zx in Fig. 8(d) is the orbital counterpart of the SHE, i.e., the OHE conductivity. Notwithstanding the analogy to the spin response, the overall shapes of χ L yx and σ Ly zx show distinct features when compared to their spin counterparts. The shape of the χ L yx profile is considerably less smooth and the flat area of σ Ly zx in the interior of the Pt layer is far more extended for all systems. In this area the orbital susceptibility, and thus the local accumulated orbital polarization, vanishes. Notably, considering the values obtained, we obtain a huge orbital response χ L yx , roughly one order of magnitude larger than the spin counterpart. Also the OHE conductivity (Fig. 8(d)) is larger than the SHE conductivity. This finding is consistent with previous calculations of the OHE in bulk metals, which obtained an intrinsic OHE that is much larger than the SHE [42,43,70]. The huge induced E ⊥ -orbital suscepti- bility is consistent with previous calculations for noncentrosymmetric antiferromagnets that obtained an OREE that was much larger than the SREE [46].
When it comes to the relative magnitude of the different configurations, E ⊥ , M ⊥ , and M || , striking differences compared to the spin responses can be observed. Here, the orbital response at the interface is dominated by the E-transverse component. We associate the Etransverse χ L yx component, as before, to the orbital accumulation caused by both the OREE and the OHE. As we will see below, in particular the OREE (intraband) contribution to χ L yx is gigantic. The M -transverse and Mlongitudinal orbital ME susceptibilities (Figs. 8(b) and (c)), are an order of magnitude smaller. Again, it is evident that the latter two orbital susceptibilities have a purely magnetic origin (M -odd) as they vanish for the nonmagnetic systems and are furthermore caused by the breaking of inversion symmetry. Similar to the case of the spin angular momentum, we identify the M -transverse component χ L xx therefore as being an orbital polarization due to the magnetic OHE.
A further significant difference between the spin and orbital ME susceptibilities is the rapid variation of the orbital ME susceptibilities in the last few layers of the Pt/Y interface. While the χ S yx component has positive values for the atomic monolayers in the vicinity of the in-terface ( Fig. 2(a)), the orbital counterpart exhibits a sign change for the two topmost layers. A similar behavior can be observed for the M -transverse components, χ S xx and χ L xx . The unusual M -longitudinal components exist, too, for the orbital ME susceptibility and conductivity, Figs. 8(c) and (f), but these quantities are, interestingly, much smaller than their spin counterparts.
The dependence of the orbital responses on the Ptlayer thickness is shown in Fig. 11 in Appendix A. Ptlayer thicknesses of about 8 monolayers provide stable values, for both the E ⊥ and M ⊥ components of the the orbital ME susceptibilities.
Orbital transport and orbital polarization at interfaces are currently only poorly understood, and first measurements are being made [71][72][73] as well as theory developed [45,74]. Our calculations show that dependence of the orbital response χ L on the magnetization direction exhibits similarities with the spin response χ S , as the nonzero components are the same for both cases. However, while the pair χ S,uz xy /χ S,ux xy differs close to the Pt/Ni interface, we find that χ L,uz xy /χ L,ux xy are virtually identical. This suggests a different, much smaller, dependence of orbital polarization on the magnetization direction at an interface.
Next, we have investigated the dependence of the induced orbital polarization on the magnetization direc-tion. In the Appendix B, in Figs. 12 and 13, we show the computed dependence of the components of χ L on the magnetization angles φ and θ, respectively. The ab initio computed angle dependence of χ L bears several similarities to the angle dependence of the spin susceptibility χ S , shown in Figs. 4 and 5. The components of the tensor χ L follow a trigonometric dependence on the angles θ and φ, but not identical to the ones given by Eqs. (27) and (28). However, for angles where a χ S ij component is zero, the corresponding χ L ij component is also zero. The symmetry of the χ L and χ S ij tensor components with respect to M , i.e., odd or even in M , is also identical. A most significant difference is the dominance of the E-transverse components χ L xy ≈ −χ L yx that are an order of magnitude larger than the other components.
Further insight is obtained from considering the dependence of the orbital ME susceptibilities on the electronic lifetime. To exemplify the broadening effect we consider the 12Pt/2Co and 14Pt systems. In the Appendix B (Fig. 15) we show the computed dependence of the nonzero orbital ME tensor elements on the lifetime broadening. Similar to the spin ME tensor elements, the intraband contributions increase linearly for decreasing broadening, whereas the interband contributions approach nonzero values (except for χ L yx ). The orbital ME susceptibilities reach however considerably larger values; for example, the intraband contribution to χ L yx is a factor of 150 larger than its spin counter part for the last Co atom of the 12Pt/2Co system. This implies that in the limit of very pure crystals the E-transverse interband χ L yx element, associated with the OREE, will dominate completely the response of this system. This conversely implies that the electrically-induced orbital polarization has then always a pure Rashba symmetry, i.e., the induced orbital moment is perpendicular to the electric field direction. This is in accordance with previous calculations for noncentrosymmetric antiferromagnets that showed that the OREE has perfect Rashba symmetry and is considerably larger than the SREE that does not give a perfectly E-orthogonal spin response [46].

Dependence on spin-orbit coupling
To investigate the dependence of the spin and orbital ME susceptibilities and conductivities we can vary the strength of the spin-orbit coupling in the calculations. To do this, we artificially introduce a SOC scaling parameter α in the DFT calculations such thatĤ 0 can be written asĤ 0 =Ĥ sc + αĤ soc whereĤ sc is the scalar-relativistic part of the Hamiltonian andĤ soc the SOC part. We find that without theĤ soc term the whole spin ME susceptibility χ S and spin conductivity σ S k ij , with indices such that ijk = 0, vanish. Obviously, as the electric field E couples only to the electron's position operatorr, relevant to the electron's orbital motion, SOC is necessary to provide a coupling to the spin. Thus, these spin quantities are completely induced by the SOC. For χ L , the story is quite different. When α is set to zero, χ L xy and χ L xy , as well as σ Lx xz and σ Lx yz , are nonzero and are actually not really affected by the modified SOC strength, a feature of the that has been noted before for the OHE conductivities [42,43] and for the orbital polarizations induced by the OREE [46].
In Fig. 9 we show comprehensive results for the layerresolved profile of χ S xy and χ L xy for 6Pt/2Ni, computed for α = 0, 0.1, 0.5, and 1, with α = 1 corresponding to the intrinsic SOC strength. It is evident from Fig. 9(a) that spin ME susceptibility is a pure SOC effect that scales linearly with the SOC. The situation is different for the orbital ME susceptibility, which exhibits practically no dependence on the SOC strength, see Fig. 9(b). Clearly, the induced E-transverse orbital polarization represents the nonrelativistic response of the electron system to the electric field. That the induced orbital moment is nonrelativistic and perpendicular to the electric field E can be recognized from the one-electron operator expression that does not require SOC. For all other spin and orbital susceptibility components, as well as for the spin conductivity tensors elements, we find that these scale linearly with the size of the SOC, i.e., these are quantities induced by the SOC.

A. Spin-orbit torque
Freimuth et al. [75] evaluated directly the SOT using a different approach to the perturbative DFT framework. While our computational method is distinct from theirs, we can evaluate the SOT T in a similar fashion. Using Eq. (21), we can write where V ↓ KS (V ↑ KS ) is the Kohn-Sham effective potential for minority (majority) spin electrons and S the equilibrium spin angular momentum. As mentioned in Sec. II D 2 this is an approximation of δB.
We define furthermore χ SOT as the SOT spin susceptibility tensor in units of TmV −1 . Since our computational approach involves quantities evaluated for each atomic site, we can access a layer-resolved B SOT .
For the thickest magnetic systems, 16Pt/2Co and 16Pt/2Ni, we find that the E ⊥ contribution to the SOT at the first (second) layer of Ni is 0.0032 (0.0020) mTcmV −1 and 0.0019 (0.0007) mTcmV −1 for Co. For the M ⊥ contribution, we find 0.0020 (0.0030) mTcmV −1 for the first (second) layer of Ni and 0.0019 (0.0020) Scaling behavior of (a) χ S yx , and (b) χ L yx as a function of the SOC scaling parameter α, calculated for the 6Pt/2Ni system (M || uz). The E ⊥ component of the spin ME susceptibility (a) scales linearly with α, and represents a SOCinduced quantity. The E ⊥ component of the orbital ME susceptibility (b) exists even without SOC. mTcmV −1 the first (second) layer of Co. These values are smaller than, but consistent with, those obtained by Freimuth et al. [27], because they used a much smaller broadening of electronic states.
The possible generation of large orbital torques has recently drawn attention [45]. However, although the orbital ME susceptibility is large, this does not automatically imply a large torque, because the induced orbital polarization can only couple to the static magnetic spin moment M via SOC. More precisely, in the commonly used formulation and implementation of DFT, the effective Kohn-Sham potential is not a functional of the orbital character, i.e., the exchange-correlation field only couples to spin angular momentum. The influence of the χ L on the SOT is nonetheless included in our calculations as the induced δL couples to S via SOC (and vice versa). Notwithstanding, theoretical efforts have recently been devoted to predicting the orbital torque [45,76] and experimental efforts are being undertaken to detecting the orbital polarization and torque and disentangling it from the spin torque [72,73,77,78].

B. Relative sizes of fieldlike and dampinglike torques
The layer-resolved SOTs are dominantly defined by the size of the spin ME susceptibility χ S . Any resulting torque can be decomposed, as customary is done, in an even-in-M FL component and and odd-in-M DL component. The linear dependence of the effective SOT magnetic field on χ S permits to compare the relative sizes of the FL and DL SOTs. Considering the case M || u z , the E ⊥ component χ S yx leads to the FL torque, while the M ⊥ component χ S xx corresponds to the DL torque. We can then quantify the relative importance of the two torques by computing the ratio A value of > 50% (< 50%) would then refer to a fieldlike (dampinglike) dominated torque. The square exponent accounts for the fact that we are comparing vectorial quantities. Note that |χ S xx | should be replaced by |χ S zx | for M || u x .
The calculated Pt-thickness dependence of this ratio for the nPt/2Ni and nPt/2Co systems is displayed in Fig. 10, for the last Pt monoatomic layer at the Pt/Y interface, as well as for the Y monolayer at the Pt interface and at the vacuum interface. It can be observed that there is virtually no change for the computed torque ratio for Pt layer thicknesses beyond eight Pt monolayers. For the Pt monolayer at the Pt/Y interface, the induced torque is to 90% composed of the FL torque component, see Fig. 10(a). For the Y monolayer at the Pt/Y interface, the torque consists for ∼ 75% of the FL component for Y = Ni and ∼ 50% for Y = Co ( Fig. 10(b)). For the Y monolayer at the Y /vacuum interface, the torque consist for ∼ 30% of the FL component for Y = Ni and ∼ 10% for Y = Co (Fig. 10(c)). This suggests that the Pt/Ni interface is more transparent to spin currents from the Pt than the Pt/Co interface, consistent with the better matching electronic structures of isoelectronic fcc Ni and Pt.
The torques resulting from the induced spin polarization on the two ferromagnetic Y monolayers will be the most important ones for the magnetization switching. The torque on the ferromagnetic layer at the vacuum interface is thus approximately both field-and dampinglike, whereas the torque at the ferromagnetic layer adjacent to the Pt layer has a larger FL contribution. As the relative contribution of the E ⊥ and M ⊥ components differs in both 3d monolayers, the direction of the total torque per monolayer will be different for each of the two Y monolayers. The calculated atomic-layer specific torques are ideally suited to investigate currentdriven magnetization switching dynamics using atomspecific Landau-Lifshitz-Gilbert spin-dynamics simulations (see e.g. [79][80][81]). Such simulations would provide layer-specific insight in how the magnetization of the ferromagnetic layers reverses under an applied electric field.

V. CONCLUSIONS
We have employed first-principles calculations to investigate the electric-field induced spin and orbital magnetoelectric susceptibility and the spin and orbital conductivity of heavy-metal/3d-metal bilayer structures. For each orientation of the 3d magnetization and the applied electric field we have shown that the susceptibility tensor and its associated conductivity tensor can be uniquely decomposed in components depending on the spatial symmetries, i.e., transverse electric E ⊥ , transverse magnetic M ⊥ , and longitudinal magnetic components M , as well as the magnetic symmetries (odd-in-M and even-in-M , respectively). Our atomic-layer specific calculations of the tensors show that all components are highly dependent on the position of the atomic layer in the considered heterostructure.
Analyzing the properties of the computed ME susceptibilities, we have identified the even-in-M , E ⊥components of χ S as spin accumulation associated with both the SHE and SREE, and the odd-in-M , M ⊥components with the magnetic SHE. We note however that our ab initio formulation uses a more general form of the SOC than the often used more elementary Rashba SOC. Extending the calculations to field-induced orbital polarizations, we have performed a similar analysis and decomposition for the orbital susceptibility tensor χ L and orbital conductivity, σ L . We have analyzed the relative importance of the different spin and orbital contributions as a function of Pt thickness. Both the out-ofequilibrium E ⊥ and M ⊥ spin responses lead to atomiclayer dependent SOTs that are of the same order of magnitude, but act in perpendicular directions. We find that the E-transverse spin accumulation is largest for the Pt layer at the Pt/3d-metal interface. The M -transverse spin accumulation, conversely, is larger at the 3d-vacuum interface. Our calculations show that both effects should be considered together when analyzing current-induced spin polarization in heavy-metal/ferromagnetic bilayer systems.
This perception is valid for electronic relaxation times that are realistic for metallic systems ( /τ ≈ 0.25 eV). For extremely pure materials, however, the intraband contribution to the E-transverse induced spin polarization will be much larger than other (interband) contributions, leading to a predominant fieldlike SOT.
Considering the electric-field induced orbital polarization, we find that the E-transverse orbital susceptibility and conductivity components are always much larger (∼ 10×) than their M -transverse orbital counterparts. In contrast to the spin counterparts, the E ⊥ -orbital ME susceptibility and OHE conductivity do not dependent on SOC. This exemplifies that the induced E-transverse orbital polarization is the primary response of the electron system to the electric field and that the other, both spin and orbital, induced polarizations are generated from the nonzero E ⊥ -orbital susceptibility by SOC. The nonrelativistic orbital ME susceptibilities are ten to a hundred times larger than the corresponding relativistic spin susceptibilities. However, although the electrically induced orbital polarization is huge, it can only couple to the equilibrium spin moment via SOC.
The computed induced spin and orbital polarizations follow trigonometric functional dependencies on the magnetization direction angles that are consistent with the M symmetry of the ME susceptibility components. Of particular interest are the large E ⊥ -orbital ME components that are practically angle independent and antisymmetric, i.e., χ L xy = −χ L yx . The induced E-transverse orbital polarization exhibits consequently a pure Rashba symmetry.
Our calculations show furthermore that there exists as well an electric-field induced spin and orbital polarization along the magnetization direction. This previously unobserved spin-orbit effect does not exert a torque on the static magnetization. We propose that it could be possible to observe this M -longitudinal effect in sensitive magneto-optical Kerr effect measurements (cf. [51]).
When the magnetization direction changes, the spin and orbital responses also change. We have shown that the magnetization direction does have a strong influence on the spin and orbital responses, but that it is possible to track the evolution of the individual components using simple, but robust, symmetry relations. This should aid the future investigation of SOT magnetization switching using atom-specific Landau-Lifshitz-Gilbert spin-dynamics simulations.
As mentioned in Sec. II C, the calculations are performed in 3 steps. First, the structures are fully relaxed with the DFT package SIESTA [52]. The cell parameters and atomic positions of the pure Pt films are relaxed until the pressure reaches values below 0.001 GPa and atomic forces on each atom are below 0.01 eV/Å. Then, the cell parameter is fixed and two monolayers of 3d elements (Ni, Co, or Cu) are added. The atomic positions are then relaxed using the same criterion as before. All SIESTA calculations are performed using a 15×15×1 Monkhorst-Pack grid [82] with an electronic temperature of 300 K. The double ζ with polarization pseudo-atomic basis set functions are used. The mesh-cutoff for real space integration is set to 250 Ry and we use the generalized gradient approximation (GGA) for the exchange-correlation functional in the PBEsol parametrization [83]. All structures contain 20Å of vacuum to avoid spurious interactions with neighboring simulation cells.
Second, once the structures are relaxed, the groundstate Kohn-Sham wavefunctions and energies are computed using the accurate full-potential, all-electron code WIEN2k [49], with spin-orbit interaction included [84]. The product between the smallest muffin-tin radius R M T and the largest reciprocal vector K max is set to R M T × K max = 8.5 and the self-consistent spin-polarized density is computed using a 30 × 30 × 1 k-points Monkhorst-Pack grid. The computed spin moments for the 16Pt/2Ni bilayer are 0.855 µ B and 0.760 µ B at Ni18 and Ni17, respectively. The spin moment on the Pt interface layer (Pt16) is 0.212 µ B . For the 16Pt/2Co bilayer the equivalent moments are 2.02 µ B , 1.942 µ B , and 0.251 µ B . The proximity induced moments in the Pt layer vanish within four layers.
Finally, the atom-resolved spin response tensors are then computed with a denser 200 × 200 × 1 k-mesh. As the WIEN2k code uses atom-centered wavefunctions, we exploit here this property to compute them in an atomprojected fashion. The simulation cell is divided into two subspaces: muffin-tin spheres around each atom, in which the wavefunction is expanded in terms of spherical harmonics, and the interstitial region in which the wavefunction is given in terms of plane waves, i.e., where the first right-hand term is the wavefunction about atom α and Ψ I (r) is the wavefunction in the interstitial. The atom-projected expected value of an operatorÔ is taken as where the integral is over the α th -muffin-tin volume and O can be replaced by any operator described in the Sec. II B.

Appendix B: Properties of the orbital responses
In this Appendix, we provide detailed calculated results for the orbital susceptibilities. We show in Fig. 11 the calculated dependence of the layer-resolved χ L tensor elements on the number of Pt monolayers n, similar to the results for the spin counterpart shown in Fig. 3. In Figs. 12 and 13 we provide the selfconsistently calculated dependence of the tensor elements on the magnetization angles φ and θ, respectively, for the 12Pt/2Co system. The angular dependencies are similar to the ones computed for the induced spin susceptibilities, shown in Figs. 4 and 5, and are even or odd-in-M , consistent with the symmetry classification. The slight asymmetry of some tensor elements, e.g., χ L xx and χ L yy in Fig. 12, can be ascribed to an angular dependence of the form sin 2φ/(a cos 2 φ+b sin 2 φ), with nonequal constants a and b.
In Fig. 14 we provide the computed magnetizationdirection dependence of the χ L for magnetization directions M || u z to M || u x . The similar shape of the curves and the very similar values of the χ L tensor elements illustrates the mapping according to the employed classification. In Fig. 15, lastly, we show the computed dependence of the intraband and interband orbital susceptibility elements on the electronic broadening, for the 12Pt/2Co and 14Pt systems.     FIG. 15. Dependence of the interband and intraband contributions of the orbital ME response χ L on the electronic broadening δ for the 12Pt/2Co and 14Pt systems. Note that the orbital ME response is roughly one to two orders of magnitude larger than the spin ME response.