Scattering-induced and highly tunable by gate damping-like spin-orbit torque in graphene doubly proximitized by two-dimensional magnet Cr2Ge2Te6 and WS2

graphene doubly proximitized by two-dimensional magnet Cr2Ge2Te6 and WS2 Klaus Zollner, ∗ Marko D. Petrović, 3 Kapildeb Dolui, Petr Plecháč, Branislav K. Nikolić, † and Jaroslav Fabian Institute for Theoretical Physics, University of Regensburg, 93040 Regensburg, Germany Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA Department of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA


I. INTRODUCTION
The spin-orbit (SO) torque [1] is a phenomenon in which unpolarized charge current injected parallel to the interface of a bilayer of ferromagnetic metal (FM) and SO-coupled material induces magnetization dynamics of the FM layer. For many possible applications of torque [2], such as nonvolatile magnetic random access memories [3] or artificial neural networks [4], it is crucial to switch magnetization from up to down along the direction perpendicular to the interface. This has led to intense search for optimal SO-coupled materials and their interfaces with FM layers which yield large SO torque while using as small as possible injected current. Thus far, minimal current density j for magnetization switching has been achieved using topological insulators (j ∼ 6 × 10 5 A/cm 2 [5]) and Weyl semimetals (j ∼ 3 × 10 5 A/cm 2 [6]), which is two orders of magnitude smaller than j required in early SO torque-operated devices employing heavy metals [1,7].
Further optimization could be achieved by using van der Waals (vdW) heterostructures [8] of very recently discovered two-dimensional (2D) ferromagnets [9,10] and 2D SO-coupled materials where current flows only through few monolayers and no Joule heat is wasted by its flow through the bulk. Furthermore, vdW heterostructures offer atomically flat and highly transparent interfaces, as well as possibility to use external manipulations-such as gating, straining and controlling coupling between 2D materials-in order to tune [11] the magnitude and ratio of field-like (FL) and damping-like Schematic view of Cr2Ge2Te6/graphene/WS2 vdW heterostructure attached to macroscopic left and right reservoirs with a small bias voltage V b between them injecting unpolarized charge current into graphene layer. A back gate voltage U bg and a top gate voltage Utg applied over a smaller "active region" are assumed to control the Fermi energy and the local on-site potential in graphene, respectively. The unit vectors of magnetic moments on Cr, Te and C atoms are denoted as mCr, mTe and mC, respectively, where only mC experiences SO torque-driven dynamics.
(DL) components of SO torque. The traditional labeling of torque components stems from how they affect dynamics of magnetization m C viewed as a classical vector-FL torque changes precession around an effective magnetic field (such as along the y-axis in the case of vdW heterostructure in Fig. 1), while DL torque competes with damping and plays a key role [12] in magnetization switching (such by directing magnetization m C towards the y-axis in Fig. 1). Controlling their ratio and relative sign is of great interest for applications since it makes it possible to tailor the switching probability [13].
The critical problem in the field of SO torquesunderstanding of competing microscopic mechanisms behind these components, especially the DL one, and control of their magnitude and ratio-remains unresolved. For example, experiments on FM/heavy-metal bilayers suggest [7] that DL SO torque is primarily generated by the spin Hall current [14] from the bulk of a heavy metal, while interfacial SO coupling (SOC) is detrimental because it generates spin memory loss [15][16][17] and, therefore, reduction of spin Hall current. Conversely, first-principles quantum transport calculations [18,19] on FM/heavy-metal bilayers find that interfacial SOC and spin Hall current can equally contribute to total DL SO torque. The purely interfacial contribution is easily identified in experiments [20] and computational studies [18,19] as the one that is independent on the thickness of heavy-metal layer, but its microscopic mechanisms remain poorly understood. Although reflection and transmission of electron spins traversing SO-coupled interface can lead to interface-generated spin currents inducing DL SO torque [21,22], both this mechanism and the spin Hall current are inoperative in 2D transport within vdW heterostructures or in FM/topological-insulator bilayers [5,23]. Thus far, it has been argued that spinindependent impurities cannot generate DL torque in purely 2D transport [24][25][26], as well as that the same impurities can completely cancel [25,26] the intrinsic Berry curvature mechanism [27] of DL SO torque as the only other possibility discussed for purely 2D transport.
In this study, we design a realistic (i.e., built atom-by-atom) vdW heterostructure monolayer-Cr 2 Ge 2 Te 6 /graphene/monolayer-WS 2 , as depicted in Fig. 1. Here Cr 2 Ge 2 Te 6 [28,29] is an example of the recently discovered 2D ferromagnets [9,10], and WS 2 is a transition-metal dichalcogenide (TMD) with strong SOC in monolayer form [30]. Unlike isolated graphene which is nonmagnetic and it hosts minuscule intrinsic SOC [31], doubly proximitized graphene within Cr 2 Ge 2 Te 6 /graphene/WS 2 trilayer offers a versatile "theoretical laboratory" to differentiate between competing mechanisms of SO torque and thereby learn how to control them.  [25,26,32] based on the Kubo formula [33] for nonequilibrium spin density ŝ(r) CD within graphene, which exerts torque ∝ ŝ(r) CD × m C (r) [8,18,[34][35][36] on the local magnetization m C (r) in graphene. Note that m C (r) is induced by magnetic proximity effect originating from Te atoms of Cr 2 Ge 2 Te 6 , so that m C m Te but these two magnetic moments are antiparallel to m Cr on Cr atoms, as illustrated in Fig. 1.
Alternatively, one can discretize continuous Hamiltonian in Eq. (1) to obtain its TB version, as provided in Appendix D, and combine it with computational quantum transport. We feed such first-principles TB Hamiltonian into the calculations [34][35][36] of current-driven (CD) part of nonequilibrium density matrixρ CD , which yields ŝ i CD = Tr spin [ρ CDŝ ] on site i of TB lattice for s = (ŝ x ,ŝ y ,ŝ z ) as the vector of the Pauli matrices and the trace being performed only in the spin space. The matrix ρ CD is computed via the nonequilibrium Green function (NEGF) formalism [37] for the Landauer setup where vdW heterostructure in Fig. 1 is divided into the left (L) semi-infinite-lead connected to "active region" which is connected to the right (R) semi-infinite-lead. Both semiinfinite leads and the active region in Fig. 1 are made of the same trilayer, assumed to be also infinite in the transverse y-direction. The leads terminate into the L and R macroscopic reservoirs kept at electrochemical potential difference µ L −µ R = eV b with small bias voltage V b in the linear-response regime. Such nonperturbative (i.e., numerically exact) approach can be used to investigate SO torque in the diffusive transport regime [18,36]; as well as in the ballistic [38,39] and quasiballistic regimes which are not accessible in the Kubo-formula-based calculations requiring nonzero electric field and momentum relaxation mechanisms. By using homogeneous potential barrier in the active region due to the top electrostatic gate [40,41] in Fig. 1, and/or short-ranged impurities generating spinindependent random potential [32,42], we demonstrate [ Fig. 4] our principal result-spin-independent scatterers can be a sole generator of nonzero DL SO torque in purely 2D transport, so they can be manipulated [32] to tune its magnitude.
The paper is organized as follows. Section II explains our density functional theory (DFT) calculations for Cr 2 Ge 2 Te 6 /graphene/WS 2 trilayer, with additional details and cases (such as Cr 2 X 2 Te 6 /graphene/WS 2 trilayers and Cr 2 X 2 Te 6 /graphene bilayers where X = {Si,Ge,Sn}) provided in Appendices A-C. Section II also provides a continuous effective Hamiltonian of graphene with parameters derived from DFT calculations, with its tight-binding (TB) version given in Appendix D. The quantum transport calculations of current-driven nonequilibrium spin densities and SO torques, based on the TB Hamiltonian (Appendix D) of doubly proximitized graphene with additional gate-or disorder-induced on-site potential, are presented in Sec. III. We conclude in Sec. IV. [ Fig. 2(b)] of Cr 2 Ge 2 Te 6 /graphene/WS 2 trilayer using Quantum ESPRESSO [43] and WIEN2k packages [44], with additional details provided in Sec. A. The band structure in Fig. 2(a) shows that the Dirac cone of graphene is preserved and located in the global band gap, so it can be probed by charge and spin transport. Figure 3 shows a zoom to the fine structure around K and K points with a fit to our continuous Hamiltonian This is derived from first-principles calculations, as well as additional symmetry arguments [45,46], and it is valid in the vicinity of both Dirac points, K and K . It reproduces the gap at the Dirac point and the spin-splitting of the bands due to proximity SOC [47,48], of valley-Zeeman and Rashba types [30,46,49], which originates from WS 2 ; as well as due to proximity magnetism [50] from Cr 2 Ge 2 Te 6 into graphene. These interactions are captured by relevant terms in our continuous Hamiltonian [Eq. (1)], with parameters provided in Table II. The vdW heterostructure we consider has broken timereversal symmetry and only C 3 symmetry. We denote v F as the Fermi velocity, and wavevector components k x and k y are measured from ±K. The valley index is τ = ±1 for ±K and pseudospin matrices areσ i , acting on sublattice space (C A , C B ), with i = {0, x, y, z} where i = 0 denotes a unit 2 × 2 matrix. For notational convenience, we usê σ ± = 1 2 (σ z ±σ 0 ). The staggered potential gap is ∆, and the parameters λ A I and λ B I denote the sublattice-resolved  intrinsic SOC. The parameter λ R is the Rashba SOC, and the proximity exchange interaction parameters are λ A ex and λ B ex . The parameter ξ describes valley exchange coupling resulting from an in-plane magnetization component [45]. The four basis states are |Ψ A , ↑ , |Ψ A , ↓ , |Ψ B , ↑ , and |Ψ B , ↓ . The continuous Hamiltonian is centered around the Fermi level at zero energy. Since firstprinciples results capture doping effects, we also introduce parameter E D (termed Dirac point energy) which shifts the global band structure.
The fitting parameters for Cr 2 Ge 2 Te 6 /graphene/WS 2 trilayer are summarized in Table II. They are in agreement with previous calculations for individual graphene/TMD [30,49] and graphene/Cr 2 Ge 2 Te 6 slabs [51]. As seen from Fig. 3, Hamiltonian in Eq. (1) with parameters from Tab. II perfectly reproduces firstprinciples results, including the spin ŝ z expectation values in both valleys. The valley degeneracy is broken, which is best manifested at the highest (lowest) spin-↑ (spin-↓) bands at K and K . This is due to the interplay of the proximity exchange and SO couplings which splits the bands in the two valleys differently. Furthermore, the Rashba SOC mixes the spin states and opens a global band gap, different for the two valleys. Appendix B also provides equilibrium spin textures within Cr 2 Ge 2 Te 6 /graphene/WS 2 trilayer. For comparison, band structure and fitting parameters for other Cr 2 X 2 Te 6 /graphene/WS 2 trilayers and Cr 2 X 2 Te 6 /graphene bilayers with X = {Si,Ge,Sn} are provided in Appendix C.

III. QUANTUM TRANSPORT CALCULATIONS OF SO TORQUE
The NEGF is computed for the active region of vdW heterostructure in Fig. 1. Since current flows only through graphene, the active region is modeled as graphene nanoribbon with armchair edges that is composed of N a = 6 dimer lines across the ribbon width W . The nanoribbon is described by the first-principles To remove the effect of the edges and mimic wide junctions along the y-axis often employed experimentally [40], sites on the lower and upper edge of the nanoribbon are connected [41] by a hopping with e ikyW phase, where k y ∈ [−π/W, π/W ). Thus, all transport quantities are then k y -points sampled to take into account an infinitely wide system [41].
The conventional spin-transfer torque [34] in spin valves and magnetic tunnel junctions is decomposed into the FL and DL components, T = T DL + T FL , each of which has relatively simple angular dependence. Conversely, the SO torque with its complex angular dependence [52] is naturally decomposed [18,36], T = T e +T o , into the odd T o and even T e components in the magnetization m C receiving the torque. They are obtained by averaging over N A atoms of the triangular sublattice A and N B atoms of the triangular sublattice B of graphene The total torque is the sum over the first Brillouin zone (BZ), T e,o = W 2π BZ T e,o (k y )dk y . The NEGF-based algorithm to split the density matrix,ρ CD =ρ e CD +ρ o CD , is given in Refs. [34,35], so that tracing Pauli matrices withρ e CD andρ o CD yields directly ŝ i e CD and ŝ i o CD and the corresponding SO torque components in Eq. (2). The magnitudes of odd and even SO torque components are plotted in Fig. 4(a),(b) in the units of eV b /A, where A is the area of a single hexagon of the honeycomb lattice. We assume that back gate and top gate can change the carrier density globally or locally, respectively, thereby making it possible to tune the Fermi energy E F of the whole device while concurrently creating a potential barrier within its smaller region, as demonstrated experimentally [40]. Without a potential barrier, U tg = 0, and in the ballistic regime with no impurities, the only nonzero component in Fig. 4(a) is T o = 0. Upon introducing U tg = 0 in still clean graphene, Fig. 4(b) reveals emergence of T e = 0. Furthermore, combined tuning of E F and potential barrier height in Fig. 4(b) modulates the ratio |T o |/|T e | by an order of magnitude, which is of great importance for applications [13]. We fix the height of the homogeneous barrier, within the gray rectangle of the active region region in Fig. 4(d),(e), at U tg = 0.1 eV.
To clarify how T e emerges from zero value in graphene without any scattering to nonzero value, Figs. 4(c) and 4(d),(e) plot spatial profiles of nonequilibrium spin densities ŝ i e,o CD [Eq. (2)] in the absence or presence of potential barrier, respectively. In Fig. 4(c), only ŝ i o CD = 0 is nonzero along the y-axis and slightly out of the plane in Fig. 5(a) for current injected along the x-axis. This is the well-known inverse spin-galvanic (or Edelstein) effect [53][54][55] in which current through 2D electron system with SOC, and thereby generated spin texture with spinmomentum locking, induces ŝ i CD = 0. But here this effect is operative in the ballistic transport regime [38,39], instead of originally considered diffusive regime [53,54]. Upon introducing impurities [leading to dashed red line in Fig. 4(c)], as random on-site potential to establish the diffusive transport regime characterized by the shot noise Fano factor F 1/3 [56], ŝ i CD acquires additional components in the other two spatial directions. The came occurs in backscattering from a barrier [leading to solid red line in Fig. 4(e)]. This means that components of ŝ i CD emerge along the x- Previously discussed mechanisms [24][25][26] for purely interfacial T e = 0 require spin-active impurities, such as magnetic or SO ones [57], that scatter spin-↑ and spin-↓ electrons differently. However, the potential barrier and/or short-ranged impurities we employ in Fig. 4 are spin-independent. Nonetheless, they can generate skewscattering [32,42] due to proximity SOC in the band structure and sublattice-staggered potentials in Eq. (1) which conspire to tilt [ Fig. 7] the equilibrium spin textures out of the plane. This occurs even when m C = 0, while the barrier also requires m C = 0. Both cases lead to generation of ŝ i CD components in all spatial directions. To connect T e and T o to traditionally discussed DL and FL torque components, we need to take into account typically complex angular dependence of SO torque observed experimentally [52] or in first-principles quantum transport studies [8,18,36]. For this purpose, we fit computational data (dots in Fig. 6) with an infinite series [8,36,52] for T e and T o vector fields on the unit sphere of orientations of m C . The non-negligible terms, compatible with the lattice symmetry, are given by The computed angular dependence of even, T e , and odd, T o , SO torque components-as the magnetization vector  m C = (sin θ cos ϕ, sin θ sin ϕ, cos θ) rotates within the xz-, xy-and yz-planes in the setup of Fig. 1 with nonzero potential barrier-is shown by dots in Fig. 6. Note that for θ = 0 • , the results in Fig. 6 correspond to the results in Fig. 4(b) at E F = −0.02 eV. The functions in Eq. (3), with fitting parameters listed in Table I, are then plotted as solid lines in Fig. 6. Hereŷ andẑ are the unit vectors along the y-and z-axis, respectively, for the coordinate system in Fig. 1. Note that prior to fitting, data sets in all three planes (xz, xy, and yz) are joined in a single continuous line for each case (with and without the onsite potential U tg ), so that they could all be fitted with a single set of fitting parameters. In Eq. (3) we recognize the lowest-order term τ o 0ŷ × m C as the FL SO torque, which is the only term (Table I) generated in the ballistic transport regime of Fig. 4(a),(c).
The usual DL SO torque [1] is τ e 0 m C × (ŷ × m C ). Higher order terms in Eq. (3) can have properties of both FL and DL torques [18]. Fitted functions in Eq. (3) are also required for using [8] computed SO torque in micromagnetic simulations [12] of the dynamics of m C (t).

IV. CONCLUSIONS
In conclusion, using first-principles calculations we design vdW heterostructure where graphene is sandwiched between semiconducting monolayers of ferromagnet Cr 2 Ge 2 Te 6 and transition-metal dichalcogenide WS 2 to acquire both SO, of valley-Zeeman and Rashba types, and exchange couplings via the respective proximity effects. Such doubly proximitized graphene offers a playground to investigate different microscopic mechanisms of SO torque. We provide first-principles continuous and tight-binding effective Hamiltonians of doubly proximitized graphene which can be used as a starting point of analytical diagrammatic or computational quantum transport calculations, respectively. In particular, using computational quantum transport, we demonstrate how DL component of SO torque can be generated solely by skew-scattering off spin-independent potential barrier or impurities in purely two-dimensional electronic transport due to the presence of proximity SOC and its spin texture tilted out-of-plane. This occurs without any spin Hall effect [7] or three-dimensional transport across the interface [21,22] based mechanisms which are naturally absent from vdW heterostructure. The electronic structure calculations and structural relaxation for graphene on Cr 2 Ge 2 Te 6 , as well as for graphene sandwiched between WS 2 and Cr 2 X 2 Te 6 with X = {Si, Ge, Te}, are performed by DFT implemented in Quantum ESPRESSO package [43]. The heterostructure supercell consists of 5 × 5 supercell of graphene whose bottom surface is covered by a √ 3 × √ 3 supercell of Cr 2 Ge 2 Te 6 , while its top surface is covered by a 4 × 4 supercell of WS 2 . We stretch the lattice constant of graphene by roughly 2%-from 2.46Å to 2.5Å-and stretch the lattice constant of Cr 2 Ge 2 Te 6 by roughly 6%-from 6.8275Å [29] to 7.2169Å. The WS 2 lattice constant is compressed by roughly 1%-from 3.153Å [58] to 3.125Å. The heterostructure supercell has a lattice constant of 12.5Å and it contains 128 atoms. The distance between WS 2 and graphene is 3.28Å, in agreement with previous calculations [30]. The distance between Cr 2 Ge 2 Te 6 and graphene is about 3.52Å, also in agreement with previous calculations [51].
Bulk vdW crystal Cr 2 Ge 2 Te 6 , which is composed of weakly bound monolayers, has Curie temperature T C 60 K and PMA. Each monolayer is formed by edgesharing CrTe 6 octahedra where Ge pairs are located in the hollow sites formed by the octahedra honeycomb. The layers are ABC-stacked, resulting in a rhombohedral R3 symmetry.
We use 24 × 24 × 1 (18 × 18 × 1 for the sandwich structure) k-point grid in self-consistent calculations to ensure converged results for the proximity exchange interaction and proximity SOC. We also perform open shell calculations that provide the spin-polarized ground state of Cr 2 Ge 2 Te 6 . The Hubbard U = 1 eV is used for Cr d-orbitals, being in the range of proposed U values for this compound [28]. The Perdew-Burke-Ernzerhof parametrization [59] of the generalized gradient approximation for the exchange-correlation functional is employed. We use an energy cutoff for charge density of 500 Ry, and the kinetic energy cutoff for wavefunctions is 60 Ry for the scalar relativistic pseudopotential within the projector augmented wave method [60] describing electron-core interactions. When SOC is included, we use the relativistic versions of the pseudopotentials.
For the relaxation of the heterostructures, we add vdW corrections [61,62] and use quasi-Newton algorithm based on trust radius procedure. In order to simulate quasi-two-dimensional systems, we add a vacuum of 20Å which avoids interactions between periodic images in the slab geometry. To determine the interlayer distances, the atoms of graphene and WS 2 are allowed to relax only along the z-axis (perpendicular to the layers), while the atoms of Cr 2 Ge 2 Te 6 are allowed to move in all directions, until all components of all forces are reduced below 10 −3 Ry/a 0 (a 0 is the Bohr radius). Table II summarizes fitting parameters in the continuous Hamiltonian in Eq. (1) for graphene within Cr 2 X 2 Te 6 /graphene/WS 2 trilayers with X = {Si, Ge, Sn}. The same parameters are also used in the TB version [Eq. (D1)] of the continuous Hamiltonian.  . The textures are plotted in the k x k y -plane in the vicinity of the K and K points at energies belonging to the conduction and valence bands, respectively. Thus, they complement the low-energy band structure shown in Fig. 3. The bands are strongly s z -polarized with some small in-plane contribution. The pattern of the in-plane contribution and the corresponding spin-momentum locking is characteristic of the Rashba SOC-induced spin texture [63].
Appendix C: Band structure of Cr2Ge2Te6/graphene bilayers The supercell of Cr 2 Ge 2 Te 6 /graphene bilayer is depicted in Fig. 8(d), where a 5 × 5 supercell of graphene is placed on a √ 3× √ 3 supercell of Cr 2 Ge 2 Te 6 . We keep the lattice constant of graphene unchanged at a = 2.46Å and stretch the lattice constant of Cr 2 Ge 2 Te 6 by roughly 4%  from 6.8275Å [29] to 7.1014Å. The resulting supercell has a lattice constant of 12.3Å and contains 80 atoms in the unit cell. The average distance between graphene and Cr 2 Ge 2 Te 6 is relaxed to 3.516Å, consistent with previous calculations [51]. In the case of Cr 2 Si 2 Te 6 /graphene and Cr 2 Sn 2 Te 6 /graphene we get interlayer distances of 3.568Å and 3.541Å, respectively.
It has been shown previously [51] that proximity induced exchange interaction in graphene due to monolayer Cr 2 Ge 2 Te 6 is an order of magnitude larger than proximity induced SOC by the same monolayer. Therefore, we first perform calculations with SOC turned off. Figure 8(a) shows DFT-calculated band structure for Cr 2 Ge 2 Te 6 /graphene bilayer without SOC. In general, the band structure resembles that of an isolated ferromagnetic semiconductor Cr 2 Ge 2 Te 6 [28], with the graphene Dirac cone located in the gap of monolayer Cr 2 Ge 2 Te 6 . The linear dispersion of graphene is nicely preserved, however the bands are spin-split due to proximity induced exchange interaction.
Without SOC, the band structure around the two valleys K and K are the same. When SOC is turned off the global band gap at E − E F = 0 is absent. Figure 8(c) shows how the band structure gets modified when SOC is turned on. As demonstrated in Fig. 8(b), the bands of the continuous Hamiltonian in Eq. (1), with parameters provided in Table III, fit perfectly the DFT-calculated low energy bands in a small energy window around the Fermi level.
In addition, Fig. 9 shows a zoom to the fine structure around the K and K points, including a fit with bands of the continuous Hamiltonian in Eq. (1) in the main text, when SOC is turned on. The fit can perfectly reproduce the DFT-computed band structure around both valleys. The fitting parameters are summarized in Table III for the calculations with and without SOC. We find that a tiny gap opens due to SOC, which is responsible for mixing the two spin components. The proximity exchange parameters are on the order of 5 meV, and proximity SOC is an order of magnitude smaller. As shown in Ta- ble III, DFT calculations with and without SOC show no significant change of the proximity exchange parameters λ A ex and λ B ex , the staggered potential induced gap ∆ and the Fermi velocity v F .
In addition to matching the bands, Fig. 9(b),(d) demonstrates that ŝ z expectation values of the valence and conduction band agree perfectly with those from the model Hamiltonian. The bands are fully s z -polarized due to the strong proximity exchange coupling induced by monolayer Cr 2 Ge 2 Te 6 . Where SOC mixes the spin-↑ and spin-↓ bands, a gap is opened and ŝ z expectation values change in sign. The ŝ x and ŝ y expectation values (not shown) are very small and show a signature of Rashba SOC in agreement with the small λ R parameter in Table III.
As a further consistency check, we recalculated the band structure with WIEN2k package [44]. The corresponding fitting parameters are also given in Table III. The only difference is that the proximity exchange parameters λ A ex and λ B ex are increased compared to the calculation with Quantum ESPRESSO package [43], which we can attribute to the different basis sets, number of kpoints, and cutoffs used. next-nearest-neighbor sites on sublattices A and B; and d ij is the vector connecting nearest neighbor sites i and j. The parameter ν ij is either +1 or −1 depending on the hopping direction. Other symbols have the same meaning as in Eq. (1), and the values of the fitting parameters are the same as in Table II. The only difference between the continuous Hamiltonian in Eq. (1) and the TB Hamiltonian in Eq. (D1) is neglect ofĤ ξ term in Eq. (D1) due to small value of ξ parameter in Table II. Another difference is in that m C as the unit vector along the proximity induced magnetization of graphene inĤ ex can point in any direction to make it possible to study the angular dependence of the SO torque in Fig. 6. In contrast, this term is fixed only along the z-axis in Eq. (1).
We note that quantum transport calculations using TB Hamiltonian in Eq. (D1) are equivalent to computationally much more expensive calculations using DFT Hamiltonian [8], as long as one remains in the linear-response transport regime driven by small applied bias voltage. The small bias voltage V b E F is indeed utilized in the calculations in Figs. 4, 5 and 6. This is due to the fact that self-consistent re-calculation of particle densities and the Hamiltonian due to current flow is nonlinear effect in the bias voltage. Trivially, one should also keep the Fermi energy E F in Figs. 4 and 6 within the energy window around the Dirac point where the continuous Hamiltonians in Eq. (1) or TB Hamiltonian in Eq. (D1) with parameters in Table II are applicable.