{'fraction_non_alphanumeric': 0.04974361885300183, 'fraction_numerical': 0.015439632173468809, 'mean_word_length': 3.878796270654482, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum interference tuning of spin-orbit coupling in twisted van der Waals trilayers

We show that in van der Waals stacks of twisted hexagonal layers the proximity induced Rashba spin-orbit coupling can be affected by quantum interference. We calculate the quantum phase responsible for this effect in graphene--transition metal dichalcogenide bilayers as a function of interlayer twist angle. We show how this quantum phase affects the spin-polarization of the graphene bands and discuss its potential effect on spin-to-charge conversion measurements. In twisted trilayers symmetries can be broken as well as restored for certain twist angles. This can be used to deduce the effects of induced spin-orbit coupling on spin-lifetime anisotropy and magnetoconductance measurements.

Multilayer stacks of two-dimensional (2D) materials, commonly referred to as van der Waals (vdW) heterostructures, have become an important platform for exploring a wide range of exciting phenomena, because they offer electronic tunability and a large parameter space. Initially, bilayer structures served as a platform to study e.g., the effects of a moire potential in graphene/hBN [1][2][3][4] the induced spinorbit coupling (SOC) in graphene/transition metal dichalcogenide(TMDC) heterostructures [5][6][7], or strong electronic correlations in magic-angle twisted bilayer graphene (MATG) [8,9]. A natural next step is to add another layer, which introduces a second twist angle and/or a further layer of different properties and this can enhance the parameter space to tune the properties of the system in several ways. An exciting opportunity is to use MATG and induce SOC into the flatbands by adding a TMDC layer to the stack [10,11]. Further examples include the recent proposal of an engineered topological phase in WSe 2 /bilayer graphene (BLG)/WSe 2 system [12], the external gate tunability of SOC through changing layer polarization in TMDC/BLG structures [13][14][15][16] and imprinting double moire potential on graphene [17][18][19]. More generally, vdW multilayers can serve as quantum simulator platform for strongly correlated physics and topological materials [20].
The proximity induced SOC plays a pivotal role in many of the above proposals. A number of recent experimental works based on weak antilocalization (WAL) [5][6][7][21][22][23][24][25] spin-lifetime anisotropy measurements [26][27][28][29], spin-Hall [30][31][32][33] and Rashba-Edelstein effect [32,[34][35][36] proved that SOC is strongly enhanced in graphene/TMDC heterostructures, which also motivated theoretical work to understand these measurements [5,[37][38][39][40][41][42][43][44][45]. Recently, Refs. [46][47][48][49] have also discussed the interlayer twist angle dependence of the induced SOC in graphene/TMDC. In particular, it has been noted [47,48] that the most general form of the induced Rashba-SOC in twisted graphene/TMDC heterostructures that obeys time reversal T and three-fold rotation C 3 symmetries can be written as where s x,y,z (σ x,y,z ) are Pauli matrices acting in the spin (sublattice) space, and τ z is a Pauli matrix acting on the valley degree of freedom of graphene. Both ϑ R and λ R are functions of the interlayer twist angle θ, which we do not show explicitly in order to ease the notations. H R differs from the usual Rashba SOC term [50][51][52] H R = λ R 2 (τ z σ x s y − σ y s x ) by a rotation of angle ϑ R in spin-space. The terms containing ϑ R appear because for a general interlayer twist angle the symmetry of the heterostructure is lowered from C 3v to C 3 . Thus Eq. (1) is valid not only for graphene/TMDC heterostructures, but for a wide range of twisted heterostructures consisting of hexagonal layers, such as heterostructures of graphene with semiconductor [53,54], ferromagnetic [55], and topological insulator [56,57] structures. This general nature of the effect is one of the main motivations for our work. However, the physical significance of this spin-space rotation was not previously appreciated and the relation of ϑ R to the interlayer twist angle θ has not been discussed. As we will show, using graphene/TMDC twisted bilayers as an example, ϑ R can take finite values and it leads to quantum interference affecting the induced Rashba type SOC in twisted trilayers. While it has been known that λ R is tunable e.g., by pressure [58] in proximity structures, to our knowledge the possibility that quantum interference can affect its value has not been considered before.
Twisted bilayers In order to obtain the angle ϑ R we use the methodology recently developed to calculate the coupling strength λ R for graphene/TMDC bilayers [48]. ϑ R is a quantum phase which depends on the interlayer tunneling between the Bloch states of graphene and the TMDC layer and on certain off-diagonal matrix elements of the intrinsic SOC of the TMDC, for details we refer to Ref. [59]. One finds that λ R and ϑ R are the absolute value and the phase of the complex Rashba coefficient λ R = λ R e iϑR , which is given by the sum of the contributions from pairs of even (e) and odd (o) bands of the TMDC:λ R = q λ R,(e,o)q e iϑ (e,o)q = q λ R,q e iϑq . Here λ R,q and ϑ q are the magnitude and the phase of the contributions of the pairs of bands. In the calculations of λ R,q and ϑ q we have used the tight-binding model of Ref. [60]. For the initial steps of the calculations it is convenient to assume that the graphene layer is rotated with respect to the TMDC layer [48]. In the final steps, we change the representation of the Hamiltonian from the coordinate system fixed to the TMDC layer to the system fixed to the graphene layer with a transformation e −iτz σz 2 θ H R e iτz σz 2 θ , see Fig. 1(a). From the explicit form of the transformed Rashba Hamiltonian one finds that the non-zero matrix elements are ∝ λ R e ±i(ϑ R +θ) , i.e., the sum of the geometric angle θ and the quantum phase ϑ R plays an important role.
Regarding our numerical calculations, the first ingredient is the tight-binding (TB) Hamiltonian of Ref. [60] for TMDCs. This TB model itself is derived from density functional theory (DFT) calculations and we use it, among others, to calculate matrix elements of the spinorbit coupling Hamiltonian and the interlayer tunneling amplitude. The results also depend on i) the position of the Dirac point of graphene within the band gap E g of the TMDC, and ii) on the value of E g . Numerical DFT calculations are known to often underestimate E g and they seem to give [37] different results from experiments [61,62] for the energy alignment of graphene's Dirac point with the TMDC bands. We performed our calculations for two parameter sets to assess how sensitive are the results on the choice of these material parameters. Since Ref. [60] does not provide information on the band alignment of graphene and the TMDC layer, we use the DFT calculations of Ref. [37] for this purpose. The second parameter set for E g and the band alignment is extracted from ARPES measurements [61,62]. Since the coupling between the layers is weak at the Dirac point of graphene, the band alignment should mainly depend on the work function difference between the two materials. Therefore we assume that it does not depend on the interlayer twist angle. This is in agreement with the recent computational work of Ref. 63.
The results for the ϑ R + θ vs θ dependence for two selected TMDCs, MoS 2 and WSe 2 , are shown in Figs. 1(b) and (c), respectively. In the case of MoS 2 , using the DFT parameter set, one can see that ϑ R + θ remains in a limited range around π as θ varies from 0 to π/3 ( Fig. 1(b)). However, if parameters extracted from ARPES measurements [61] are used then ϑ R + θ covers the entire range [0, 2π]. For WSe 2 one finds that ϑ R + θ covers all of [0, 2π] ( Fig. 1(c)) and the results obtained from the two parameter sets qualitatively agree. The difference between the results for the two materials can be mainly traced back to the different energy alignment of the Dirac point in the TMDC band gap. One can also note in Figs. 1(b,c) that for θ = lπ/6, l = 0, 1, . . . one finds ϑ R + θ = nπ, where n is an integer. We give a detailed discussion of how this result for ϑ R + θ follows from our theoretical method in Ref.
[59], but already note at this point that for θ = lπ/6 the vertical mirror planes of the graphene and the TMDC lat-tice line up and the system, as a whole, has C 3v symmetry. In this case the Hamiltonian of the induced Rashba SOC reads H R = (−1) n+1 λ R (θ) 2 (τ z σ x s y − σ y s x ), hence H R simplifies to the form used previously in the literature [37,50]. We find that n can be, in general, both odd and even (see Figs. 1(b),(c)), which means that λ R (θ) can acquire a negative sign as θ is changed. Interestingly, when ϑ R + θ = nπ, the spin-orbit field s(q) = ( ŝ x , ŝ y , ŝ z ) T is not tangential to the Fermi surface as in the case of usual Rashba SOC (cf. Fig. 1(e) and Figs. 1(f),(g)). Instead, one can show that the inplane component ( ŝ x , ŝ y ) T is rotated by an angle ϑ R + θ with respect to the tangential direction. As we will discuss, this might have consequences on the interpretation of experimental results.
Twisted TMDC/graphene/TMDC trilayers Adding another TMDC layer, as shown in Fig. 1(d), introduces a second interlayer twist angle and the two twist angles θ (b) and θ (t) for the bottom and top TMDC layers allow an even broader control of the induced SOC in graphene.
Since the layers are only weakly coupled, the effective graphene Hamiltonian is H orb + H Here H vZ = λ vZ τ z s z is the Hamiltonian of the induced valley Zeeman SOC in the graphene layer. Note, that the contributions H R have a different sign. As a simple physical explanation, consider the case when the two TMDC layers are perfectly aligned, e.g., Then the graphene layer is horizontal mirror plane of the whole stack, which dictates that the Rashba SOC must vanish. (A more microscopic argument is given in Ref. [59].) One can define the complex Rashba coefficient for the trilayer system (tls) bȳ and its magnitude λ In terms of these quantities the induced Rashba type SOC can be written as H The importance of the phase e i(ϑ R +θ) discussed for bilayers becomes now more clear: it follows from Eq. (2) that the strength λ (tls) R of the induced Rashba SOC in trilayer stacks can be affected by quantum interference effects if ϑ R + θ (t) are non-zero. This can be interpreted as an interference of the virtual hopping processes to the two TMDC layers that give rise to the induced Rashba SOC. Calculations for the twist angle dependence of λ R have already been performed in Refs. [47,48], therefore we do not show these results here, see Ref.
where l is an odd integer, the trilayer stack is inversion symmetric. On the other hand, for an even l the stack has a horizontal mirror plane. Therefore the induced Rashba SOC must vanish for any integer l, as it can be seen in Figs. 2(a),(b). For WSe 2 encapsulation ( Fig. 2(a)) the maximum of λ (tls) R is found for |θ (b) − θ (t) | = (2l + 1)π/6. The maximum value of λ (tls) R is around 70% larger than in the graphene/WSe 2 case. Thus, double-sided encapsulation can significantly enhance the induced Rashba SOC. Surprisingly, when using MoS 2 for double encapsulation ( Fig. 2(b)), we find that λ (tls) R basically cannot be enhanced above the value obtained for one sided proximity effect. This can be understood by considering the explicit dependence of λ [59] for details. The origin of the extended regions where λ (tls) R is very small is due to the fact that for those twist angles λ R and ϑ R change slowly in both layers and therefore they approximately cancel in Eq. (2). In Figs. 2(c,d) we show the phase ϑ (tls) of Eq. (2), which determines the winding of the SOC field s(q) in the case of double encapsulation. The apparent diagonal lines correspond to λ (tls) R = 0 where ϑ (tls) is not defined. For completeness, we also discuss the twist angle dependence of the induced valley Zeeman SOC in TMDC/graphene/TMDC trilayers, see Figs. 2(e,f). The strength of the valley Zeeman type SOC is simply given by λ vZ (θ (t) ). If the two TMDC layers are (nearly) aligned, they can double the strength of the induced valley Zeeman SOC. When the whole stack has inversion symmetry, the effect of the two layers cancel and λ (tls) vZ = 0. The valley Zeeman SOC also vanishes along the lines θ (b) + θ (t) = π/3, π, this is a combined effect of time reversal and three-fold rotation symmetry of the TMDC layers. We also note that λ In the preceding discussions of graphene/TMDC and TMDC/graphene/TMDC heterostructures we have neglected a possible lateral shift of the graphene layer with respect to the TMDCs and moiré effects. Regarding the lateral shift, it does not affect our results [59]. Considering the moiré effects, in graphene/TMDC bilayers they are present only at energy scales larger than 2 eV[59, 61], i.e., they are negligible for energies close to the Dirac point of graphene. The situation might be different in TMDC/graphene/TMDC trilayers. Based on previous work of Refs. 65 and 66 on hBN/graphene/hBN trilayers, one may expect that supermoiré effects may become important when the TMDC layers are nearly aligned: [59] for further details). The discussion of supermoiré is beyond the scope of the present work.
Experimental predictions Several recent works [26,28,29] measured an anisotropy of the out-of-plane τ ⊥ and inplane τ spin lifetimes in graphene/TMDC heterostructures which can be interpreted as direct proof of induced SOC in graphene. Namely, according to the theoretical calculations of Ref [39], if there is a strong intervalley scattering in graphene, then the ratio of the spin lifetimes is given by τ ⊥ /τ = (λ vZ /λ R ) 2 (τ iv /τ p )+1/2, where τ iv (τ p ) is the intervalley (momentum) scattering time. (For ultraclean samples with SOC comparable to or larger than the disorder-induced quasiparticle broadening, a qualitatively different spin relaxation anisotropy is derived in Ref. [67].) Let us consider a WSe 2 /graphene/WSe 2 heterostructure and assume that θ (b) is kept fixed while θ (t) is changed. Note, that the ratio λ vZ /λ R can be tuned in a wider range in trilayer structures than in bilayers. For example, λ R is never zero for graphene/TMDC whereas one can choose θ (b) and θ (t) such that λ Using τ iv /τ p ≈ 5 as in Ref [39], we plot τ ⊥ /τ as a function of θ (t) for θ (b) = 0 • and θ (b) = 30 • in Fig. 3(a). When θ (b) = 0 and 1 • < θ (t) 15 • then λ (tls) R becomes small but λ (tls) vZ is finite, therefore τ ⊥ /τ strongly increases as a function of θ (t) . This enhancement happens before supermoiré ef-fects might become important for θ (t) 1 • . In contrast, if θ (b) = 30 • then τ ⊥ /τ remains finite for all θ (t) angles. This dramatic difference in the θ (t) dependence of τ ⊥ /τ is clearly a consequence of the wider tunability of a) Calculated spin lifetime anisotropy in a TMDC/graphene/TMDC heterostructure as a function of θ (t) for θ (b) = 0 • (blue) and θ (b) = 30 • (green). b) white color indicates the regions in the (θ (t) , θ (b) ) space where the REE is expected to vanish for a WSe2/graphene/WSe2 trilayer. In the calculations for a) and b) we used DFT parameters [64]. c) Conductance G (in units of 2e 2 /h) through an npn junction as a function of magnetic field Bz and Fermi energy EF for λvZ = 2.6meV. d) Cross section of c) at fixed Bz = 0.436T. λvZ = 0meV (blue) and λvZ = 2.6meV (orange).
One can expect that charge-to-spin conversion measurements should also be affected by the interlayer twist. Let us consider a graphene/TMDC bilayer and for simplicity assume that i) λ vZ λ R and ii) a dc electric field is applied along thex direction. In steady state theŷ component of the spin density is given by S y = 1 2 d 2 q 4π s y (q)δf q where δf q is the deviation of the quasiparticle distribution with respect to equilibrium. Using the semiclassical argument given in Ref. [41], for Fermi energies close to the Dirac point E D such that |E F − E D | < λ R (θ), i.e., when there is a single Fermi surface, one finds that S y ∝ cos(ϑ R + θ). This suggests that the Rashba-Edelstein effect (REE) may vanish (when ϑ R + θ = (2l + 1)π/2) even though λ R (θ) is not required to be zero by symmetry for any θ. S y can also change sign as a function of the interlayer twist angle, because ϑ R + θ ∈ [0, 2π], see Figs. 1(b),(c). Regarding the more realistic situation when there are two spinpolarized Fermi surfaces in graphene for |E F −E D | > λ R , looking at Figs. 1(g) one can understand that S y can be zero as a function of twist angle in this case as well. Furthermore, in twisted TMDC/graphene/TMDC trilayers, the REE can vanish either because λ (tls) R = 0 or because ϑ (tls) = (2l + 1)π/2, which can be expressed compactly as Re λ (tls) R = 0. As an example, the region where this condition is fulfilled in the (θ (t) , θ (b) ) space for WSe 2 /graphene/WSe 2 trilayers is shown in Fig. 3

(b).
While the spin-lifetime anisotropy measurements require diffusive samples, in Figs. 3 (c)-(d) we show an example of how strong SOC can affect ballistic transport properties. Highly transparent np and npn junctions in graphene have recently been realized in several experiments [68][69][70][71][72][73] and Fabry-Perrot type interference measurements [72,73]. Assume now that a TMDC/graphene/TMDC trilayer is tuned with interlayer twist angles such that the Rashba SOC is switchedoff and simultaneously the valley Zeeman SOC enhanced. We calculate the conductance[59] through a smooth graphene npn junction [74] as a function of outof-plane magnetic field B z and Fermi energy E F for λ vZ = 2.6meV, see Fig. 3(c). One can observe that for |B z | > 0.3T the high-conductance ridges are split. This is apparent in Fig. 3(d), where we compare the cases when λ vZ = 0 (blue) and λ vZ = 2.6meV (orange) for a fixed magnetic field B z = 0.436T. Our calculations also indicate that the Rashba type SOC does not have a similar effect on the conductance ridges (not shown). An additional experimental probe of the twist angle dependent SOC may be the measurement of the reflection of electrons at a planar junction, which is briefly discussed in Ref. [59].
Summary We found that the induced Rashba type SOC in twisted hexagonal bilayers can be parameterized by the strength λ R and a spin-rotation angle ϑ R . This latter can lead to interference effects for λ R in trilayer heterostructures. We also calculated the valley Zeeman SOC in twisted TMDC/graphene/TMDC trilayers. Finally, we discussed how the interlayer twist angle dependence of the induced SOC can be deduced from spinlifetime anisotropy, charge-to-spin conversion and magnetrotransport measurements.
In this work we have neglected possible lattice relaxation in vdW heterostructures [75]. An important future direction would be to study its effect on the results presented here.
Acknowledgments We acknowledge helpful conversations with Péter Makk and Jaroslav Fabian. This research was supported by the Ministry of Innovation and Technology and the National Research, Development and Innovation Office (NKFIH) within the Quantum Information National Laboratory of Hungary. P.R. and A.K. were also supported by the ELTE Institutional Excellence Program (TKP2020-IKA-05), and the Hungarian Scientific Research Fund (OTKA) Grants No.K134437 and No.NN127903 (Topograph FlagERA project). A.K. acknowledges support from the Hungarian Academy of Sciences through the Bolyai János Stipendium (BO/00603/20/11) as well.
Note added During the preparation of this manuscript we have become aware of Ref. [63], where the authors discuss the twist angle dependence of the induced SOC in graphene/TMDC structures using DFT calculations.

I. SYMMETRY PROPERTIES OF THE PHASE ϑR
In order to be self-contained and to set the stage, in Sec. I A and I B we carefully re-derive previous results of Ref. [1] for the Hamiltonian H gr R , which describes the proximity induced Rashba SOC in graphene. We then present two short, but important discussions. First, in Sec. I C we argue that the sign of the SOC coupling λ R depends on the layer stacking order and this is important in trilayer stacks. Next, in Sec. I D we discuss the effects of the interlayer twist on H gr R . Our most important results in this section follow in Sec. I E. We prove that the results for the quantum phase ϑ R simplify for interlayer twist angles θ = 0, π/6, π/3, where the bilayer stack has C 3v symmetry. Finally, in Sec. I F we give a further discussion of the results obtained in Sec. I E.
A. Preliminaries

Lattice vectors
Both graphene and the TMDC layer have hexagonal lattices. We define the lattice vectors a 1 = Fig. 1(a)) and the corresponding primitive reciprocal lattice vec- Here a = a gr (a = a tmdc ) is the lattice constant for graphene (TMDC). We will distinguish the reciprocal lattice vectors b 1,2 of graphene from the TMDC ones by using the notation b 1,2 for the latter. The position of the A and B sublattice in the unit cell is given by The metal (chalcogen) atoms occupy a position corresponding to the A (B) sublattice in the unit cell of the TMDC.

Formalism to calculate the induced SOC
As it was shown in Ref [1], an important contribution to the proximity induced Bychkov-Rashba type SOC in graphene comes from virtual interlayer tunneling processes to the TMDC layer. In third order perturbation theory it is given by Here X = {A, B} runs over sublattice indices of graphene, s, s are spin indices and b = b are band indices running over the bands of the (isolated) TMDC layer. Spin-flip process can be facilitated by certain offdiagonal elements of the intrinsic SOC matrix of the monolayer TMDC denoted by (H soc ) bs,b s The spin-flip off-diagonal matrix elements are allowed between pairs of bands b, b if one of the bands is symmetric (even) and the other one is antisymmetric (odd) with respect to reflection on the horizontal mirror plane of the TMDC (see, e.g., Ref. [2] for further discussion of the SOC in monolayer TMDCs). We assume that the graphene layer is rotated with respect to the TMDC layer by a twist angle θ. Therefore, as indicated in Fig. 1(b), the BZ of graphene is also rotated. For each interlayer twist angle θ the right-hand side of Eq. (1) needs to be evaluated for wavenumbers τ k j in the TMDC Brillouin zone (see Fig. 1 Here τ K θ , τ = ±1 denote the wavevector to the Dirac point of graphene and G θ and (G ) are reciprocal lattice vectors of graphene (TMDC). As indicated in Fig. 1(b), for each τ K θ there are three wavevectors τ k j and corresponding reciprocal lattice vectors G θ j and G j that satisfy the quasimomentum conservation condition. These τ k j vectors are related by 2π/3 rotations. Explicitly, one finds G θ (τ k j ) denote TMDC bands energies and E gr D is the energy of the Dirac point of graphene. (T τ k j ) Xs,bs is the tunneling matrix between Bloch states at τ K θ to states at τ k j . We assume spin-independent tunneling and therefore (T τ k j ) Xs,bs is diagonal in the spin space. (T τ k j ) Xs,bs can be written as Here H orb is the orbital part of the Hamiltonian of the system (without the SOC), Ψ (X) tmdc (k j , r) are Bloch wavefunctions of graphene and of band b of the TMDC layer, respectively. A possible lateral shift of the graphene lattice with respect to the TMDC lattice is denoted by d 0 . The phase φ X j is defined as φ X j = t X · G j . The tunneling amplitude t b (τ k j ) to band b of the TMDC can be written as where t X (t X ), X = X is the vector pointing to the chalcogen (metal) atom position in the unit cell of the TMDC layer, and t p b (τ k j ) (t d b ((τ k j )) describes the interlayer tunneling between graphene and the p orbitals of the chalcogen (d orbitals of the metal) atoms. The calculation of the amplitude t p,d b (τ k 1 ) in terms of materials dependent parameters will be discussed in Sec. II A.
Here λ soc characterizes the strength of the intrinsic SOC in the TMDC layer.Ĥ soc ∼L xŜx +L yŜy +L zŜz = 1 2 (L +Ŝ− +L −Ŝ+ ) +L zŜz ,L x,y,z are angular momentum operators,L ± =L x ± iL y ,Ŝ x,y,z = 2 s x,y,z are spin operators and the Pauli matrices s x,y,z act on the spin degree of freedom. The complex numbers α One can make use of the fact that the τ k j vectors are related by 2π/3 rotations, which is also a symmetry of the TMDC lattice. Following Ref. [3] one may write Ψ tmdc (τ k j , r) where R ±2π/3 denotes rotation by ±2π/3. Substituting this into Eq. (4) and taking into account that where α e,o (τ k 1 ).

B. Derivation of ϑR and λR
It is instructive to consider the contribution of a single pair of b = e and b = o bands in Eq. (1). Using Eq.(2), one finds As explained in Sec.
There are a number of possible choices for the vectors t A,B and for the primitive reciprocal lattice vectors b 1,2 of graphene, and the latter determine G j . For a given j the phase φ A,B j depends on these choices, but Hamiltonians. Importantly, U(τ ) does not depend on the geometric twist angle θ. Using Sec. I A 1, the M j matrices read ] mod π and, as it will be explained in more detail in Sec. I C, the index µ = 0, 1 gives the sign of the induced Rashba SOC which depends on the layer stacking. Note that the phase factors e iτ G θ j ·d0 , which encode the effect of a lateral shift between the layers in Eq. (2), have dropped out of the calculations. This means that our results for the induced Rashba SOC do not depend on such shift. The same conclusion holds for the valley Zeeman SOC as well, which depends on |t b (τ k j )| 2 , see Ref. 1.
Because of the C 3 symmetry of the system, t e,o (τ k j ) does not depend on the index j and from time reversal symmetry follows that t e,o (−k 1 ) = t * e,o (k 1 ). Therefore, by taking into account Eq. (6), where Λ e,o (τ Here Im[. . . ] denotes the imaginary part.θ e,o (k 1 ) is a quantum phase which depends on i) the phase η(k 1 ) of the interlayer tunneling amplitude between the Bloch states of graphene and the TMDC layer and ii) on off-diagonal matrix elements of the intrinsic SOC of the TMDC layer defined in Eq.(4). Due to C 3 symmetry, the denominator of Eq. (7) is the same for all j and because of time reversal symmetry it does not depend on τ . Therefore the denominator can be pulled out of the summation. Substituting Eqs. (8) and (9) into Eq. (7) and performing the summation over j one finds Here we defined and λ R,(e,o) (k 1 ) is given by [1] Using the τ z Pauli matrix, one may rewrite Eq. (11) in a more compact form: Note, that in the derivation of Eqs. (12)- (14) we only made use of time reversal symmetry and that the k j vectors are related by 2π/3 rotations. Therefore these results are valid for arbitrary interlayer twist angle θ. On the other hand, we considered the contribution of a single pair of e and o bands of the TMDC to the induced Rashba type SOC. As Eq. (1) shows, one needs to sum over all q pairs of (e, o) bands and one may define the complex Rashba SOC coefficient bȳ In terms of its magnitude λ R (k 1 ) = |λ R (k 1 )| and phase ϑ R (k 1 ) = Arg[λ R (k 1 )], the Hamiltonian of the induced Rashba SOC can be written as Note, that both λ R (k 1 ) and ϑ R (k 1 ) depend implicitly on θ through k 1 , as shown in Fig. 1. We will discuss the θ dependence of ϑ R in more detail in Sec. I E.
C. Staking order dependent sign of λR In our theory the sign of λ R depends on the stacking order of the graphene and the TMDC layers, i.e., whether is graphene below or above the TMDC layer. Such a sign has no physical significance in TMDC/graphene bilayer stacks, but it is important in TMDC/graphene/TMDC trilayers, because it means that the contributions of the two TMDC layers to the induced Rashba SOC add up with an opposite sign.
As explained below Eq. (8), the induced Rashba SOC is a second order process in interlayer tunneling. The tunneling amplitude is given by t e,o (τ k j ) = t e (τ k j )t * o (τ k j ), i.e., it depends on the product of the tunneling amplitudes to an e and to an o band of the TMDC. Ψ (e) tmdc (r) of an e band of the TMDC is even with respect to the mirror plane that contains the metal atoms of the TMDC layer. However, the graphene p z orbitals have a different sign above and below the graphene layer. Therefore one can expect that the matrix element given in Eq. (2) between the graphene Bloch wavefunction Ψ gr (r) and Ψ (e) tmdc (r) should have a global sign difference depending on the stacking order of the two layers, i.e., whether the graphene layer is above or below the TMDC layer. On the other hand, similar consideration suggests that the sign of the amplitude t o (τ k j ) between Ψ gr (r) and Ψ (o) tmdc (r) does not depend on the stacking order. Therefore, there will be a stacking order dependent sign factor in t e,o (τ k j ) = t e (τ k j )t * o (τ k j ), which is described by the index µ = 0, 1 introduced below Eq. (8). The stacking order dependence of t e (τ k j ) and t o (τ k j ) can be explicitly shown if they are calculated e.g., using the two-center Slater-Koster parametrization of the corresponding transfer integrals, see Refs. [1,4] for details.
Because of this stacking order dependent sign, the induced Rashba type SOC in TMDC/graphene/TMDC trilayers, which is given by the sum of the contributions from the two TMDC layer, reads H Taking into account the interlayer twist angle dependence of λ R (see Sec. II C) and of θ + ϑ R (shown in Fig.1 of the main text) one can then easily see that H (tls) R vanishes whenever the TMDC/graphene/TMDC stack is either inversion symmetric or has a horizontal mirror symmetry. σ y δk y ) for the orbital Hamiltonian of graphene in the x−y coordinate system, where δk = (δk x , δk y ) is a small wavevector measured from the τ K point. In the x − y system H gr orb is given by H gr,θ orb = e −iτz σz 2 θ H gr orb e iτz σz 2 θ . The total effective Hamiltonian of graphene is therefore where we have also taken into account the valley-Zeeman induced SOC, which is described by H gr vZ is invariant under rotations around theẑ axis and has the same form in both the x − y and x − y coordinate systems. One may perform a unitary transformation e iτz σz 2 θ H gr,θ eff e −iτz σz 2 θ . This transformation leaves H vZ unchanged and one finds H gr eff = H gr orb + e iτz σz 2 θ H gr R e −iτz σz 2 θ + H gr vZ .
It is convenient to perform another unitary transformation H gr eff → σ x H gr eff σ x . This changes the orbital Hamiltonian v F (τ z σ x δk x − σ y δk y ) → v F (τ z σ x δk x + σ y δk y ), leaves H gr vZ unchanged, and transforms the Rashba Hamiltonian to Eq. (20) shows more explicitly the effects of interlayer rotation on the graphene Rashba HamiltonianH gr R = λ R 2 (τ z σ x s y − σ y s x ) derived previously in Refs. [5,6]: i)H gr R is rotated in real space by θ, an ii) if ϑ R = 2mπ, m integer, then a rotation in spin-space appears.

θ = 0 twist angle
In Fig. 3 we show the wavevectors k 1,2,3 that satisfy the quasimomentum conservation as well as the three reflection planes R (1,2,3) v that the system possess in addition to the C 3 symmetry. First, we will derive a useful relation for the tunneling amplitude t e,o (k 1 ) = t e (k 1 )t * o (k 1 ), defined below Eq. (8). Using the definition t e(o) (τ k j ) given in Eq.
• The results of Sec. I E 1-I E 3 were obtained for an arbitrary pair of e and o bands. This means, that using Eq. (15) one finds ϑ R (θ l ) = n(θ l )π − θ l , where n(θ l ) is an integer.
• According to Eq. (20), the non-zero matrix elements of H gr R are ∝ λ R (θ)e ±i(θ+ϑ R ) . For θ = θ l one finds that these matrix elements are λ R (θ l )(−1) n(θ l ) . This means that for interlayer twist angles where the C 3v symmetry of the stack is restored, the Hamiltonian of the induced Rashba SOC simplifies to which is, apart from the sign (−1) [n(θ l )+1] , the well-known result of Refs. [5,6]. This shows our results are in agreement with general expectations based on the symmetry of the system.
• For twist angles θ ∈ (θ l , θ l+1 ), i.e., when the stack has only C 3 symmetry, ϑ R is a continuous function of θ (through the wavenumbersk j ). Therefore the matrix elements λ R (θ)e ±i(θ+ϑ R ) of H gr R are complex numbers.

II. DETAILS OF THE NUMERICAL CALCULATIONS FOR ϑR, λ R AND λ vZ
In this section we discuss questions related to the choice of material and other parameters in our numerical calculations and compare some of our results to corresponding results in Ref. [1].

A. Interlayer tunneling parameters
In Ref. [1] it was argued that for the description of the tunneling between the monolayer graphene (MLG) sheet and the MoS 2 layer it is sufficient to consider only the p orbitals of the chalcogen atoms and one can neglect the d orbitals of the metal atoms. Namely, the the d orbitals are farther away from the carbon atoms of MLG and therefore their overlap with the carbon p z orbitals should be small. We have extended the calculations of Ref. [1] by taking into account contributions from tunneling processes involving the d orbitals of the metal atoms as well. The results shown in the main text have been obtained in this way.
As already mentioned in Sec. I A, the tunneling matrix (T (τ k j )) b can be written as When only the p orbitals are considered in the interlayer tunneling, then the interlayer tunneling amplitude t p b (τ K θ , τ k j ) ≡ t p b (τ k j ) between the K θ point of graphene and electronic states at wavenumber k j in band b of the TMDC is given by Here c bx (τ k j ), c by (τ k j ) and c bz (τ k j ) are the complex amplitudes of the p x , p y and p z orbitals of the chalcogen atoms at the BZ point k j in band b of the TMDC. As it was shown in Ref. [1], t p b (τ k 2,3 ) = t p b (τ k 1 ). According to Eq. (35), two independent tunneling coefficients are needed: t p and t p ⊥ , and Ref. [1] obtained them from a fitting procedure using previous results [7] of DFT band structures calculations for MLG/TMDC heterostructures. The values of these parameters were estimated to be t p ≈ t p ⊥ ≈ 100 meV. Considering now the d orbitals of the TMDC metal atoms one can derive a similar expression for the tunneling amplitude t d b (τ k j ) between MLG and the d orbitals: Here c b,d (τ k 1 ) (d = {x 2 − y 2 , xy, xz, yz, 3z 2 − r 2 }) are the complex amplitudes of the d orbitals of the TMDC metal atoms at the BZ point k 1 in band b. The derivation of Eq(36) involves the same steps as that of Eq. (35) and therefore we do not give further details here. As one can see, t d b (τ k 1 ) involves three more overlap parameters: t d , t d || , t d ⊥ . Therefore there are altogether five t p,d parameters that describe the overlap between graphene's p z orbitals and the p, d orbitals of the TMDC layer. The approach used in Ref. [1] to estimate t p and t p ⊥ is not applicable to estimate all five t p,d parameters. Therefore we used the Slater-Koster method to re-calculate the tunneling coefficients both for the chalcogen p orbitals and the d orbitals of the metal atoms. The steps needed to calculate t p , t p ⊥ , t d , t d || , and t d ⊥ in terms of the set of parameters of the Slater-Koster method are explained in Refs. [1,4]. We only note that the parametrization of the transfer integrals is based on the Harrison's model [8] and the numerical values of the necessary parameters, with the exception of η C−Se , are listed in Ref. [4], see Eqs(4), (5) therein. All our calculations for MLG/MoS 2 are based on the parameters given in Ref. [4].
However, for the calculations for MLG/WSe 2 we also needed the Harrison's model parameter η C−Se . We have fixed its value in the following way. We took the values λ vZ and λ R of the induced valley Zeeman and Rashba SOC from Ref. [7], where they were obtained by fitting band structure calculations of graphene/WSe 2 heterostructures at θ = 0 twist angle. Using the methodology of Ref. [1] to calculate λ vZ and λ R , we adjusted the value of η C−Se such that we obtain a reasonable agreement with Ref. [7], for θ = 0. We found that η C−Se = 0.748 minimizes the deviation to ≈ ±20% between our results for λ vZ , λ R and those of Ref. [7].
To summarize, we used tunneling coefficients shown in Table I. Regarding the t p coefficients for MoS 2 , they are indeed not very different from those estimated in Ref. [1] (≈ 100 meV) using a different method,  but the sign of t p || is opposite to the one given in Ref. [1]. This can be traced back to the fact that the TB model of Ref. [9] uses a different orientation of the TMDC lattice with respect to Ref. [1]. This gave rise to an overlooked k → −k transformation in reciprocal space, which resulted in a sign change in the estimation of the tunneling parameter. One can also see that the values of the t d coefficients are indeed much smaller than t p and t p ⊥ . On the other hand, the complex amplitude factors c b,x 2 −y 2 , c b,xy , c b,xz , c b,yz , and c b,3z 2 −r 2 turn out to be significantly larger than c bx , c by and c bz in large regions of the BZ. Therefore t d b (τ k 1 ) is found to be small, but non-negligible with respect to t p b (τ k 1 ).

B. Further parameters
In addition to t b (τ k j ), there are a couple of further parameters that enter the numerical calculations. We will briefly discuss them here.
• The value of the band gap E g of the TMDC.
As briefly discussed in the main text, one can take its value either from experiments (when available) or from previous theoretical works. Regarding the latter, since we use the TMDC TB model of Ref. [9] to calculate the band structure, the spin-orbit coupling Hamiltonian matrix elements and the interlayer tunneling amplitude, we also used the E g of this model. (We implemented the TB model in the Kwant code [10]). Although the model of Ref. [9] itself is based on DFT calculations, the value of E g is different from what one can extract from the calculations of Ref. [7] that were performed for MLG/TMDC supercells at ϑ = 0. For example, in the case of MoS 2 the model of Ref. [9] has a band gap that is 17% larger than the corresponding E g given in Ref. [7]. See Table II for the E g values used in this work. When f G = 0, the Dirac point is aligned with the TMDC valence band edge, for f G = 1 the Dirac point has the same energy as the TMDC conduction band edge. When available, the value of f G can also be taken from experiments [11,12]. See Table II for the f g values used in this work.
• The number of bands in the model for the TMDC layer.
The explicit expression to calculate λ vZ involve a sum over contributions of individual bands. Similarly, in order to calculate ϑ R and λ R , one needs to sum over pairs of even (e) and odd (o) bands. Ref. [1] used the approximation that the conduction band (CB) and the valence band (VB) were taken into account for λ vZ , and three pairs of e-o bands for λ R . As it will be shown below, this already leads to qualitatively good results in most cases, see Sec. II C below. An exception is the value of λ R for θ ≈ π/6, where, due to band crossings and near-degeneracies of certain e bands, more pairs of e and o bands need to be taken into account. The TB model that we use involves 11 bands of the TMDC layer. Unless otherwise indicated, we use the contributions of all 11 bands to calculate λ vZ , and 30 pairs of e and o bands for the ϑ R and λ R calculations.
Eg (DFT) Eg (exp) fG (DFT) fG (exp) MoS2 1.807 eV [9] 2.0 eV [11] 0.974 [7] 0.55 [11] WSe2 1.638 eV [9] 1.95 eV [12] 0.161 [7] 0.426 [12]  • The TB model of the TMDC 13 We also note that the TB model of the TMDC appears to be less accurate for the bands above the band gap [9]. This can add further uncertainty to the results, especially in the case graphene/MoS 2 , where f G = 0.974, i.e., the relative contributions of e-o band pairs above E g is larger than the contribution coming from the valence bands.
These details, in addition to those already discussed in Sec. II A, also need to be considered when comparing our results to related previous works [1,4,7,13].

C. Numerical calculations of the induced SOC in MLG/TMDC twisted bilayers
We now show our results for the twist angle dependence of λ vZ and λ R for graphene/MoS 2 and graphene/WSe 2 . In the former case we also compare our results to those of Ref. [1]. For the numerical calculations of λ R we evaluated Eqs. (13) and (15). To calculate λ vZ we used the result of Ref. [1]: where t b (k 1 ) is defined in Eq. (3) and the summation runs over all bands b of the tight-binding model.

Comparison to the results of Ref. [1]
First, we briefly discuss how the tunneling parameters shown in Table I modify the results for tunneling amplitude t b compared to the corresponding results in Ref. [1]. In Fig. 6 we plot the tunneling strength |t c | 2 and |t v | 2 to the CB and VB of MoS 2 as a function of the twist angle θ using the parameters from Table I  (solid lines). We also show the corresponding results of Ref. [1], where only t p ≈ t p ⊥ were non-zero. As one can see in Fig. 6, the present calculations tend to yield larger tunneling strengths. However, the twist angle dependence of |t c | 2 and |t v | 2 remain qualitatively the same.
In Fig. 7(a) we show calculations for λ vZ vs θ for graphene/MoS 2 using the DFT parameters for E g and f G , see Table II. The dashed curve in Fig. 7(a) again corresponds to the result of Ref. [1], which uses the contributions of the CB and the VB only. The thin solid line is obtained using the tunneling amplitudes from Table I and only tunneling to the VB and the CB is taken into account. One can see a pronounced increase of λ VZ around θ = 18 • with respect to Ref. [1], which is mainly due to the larger tunneling strengths, as shown in Fig. 6. Nevertheless, the curves remain qualitatively very similar. The thick solid line shows the result that we obtain by using the tunneling amplitudes from Table I Table II. Different lines correspond to different tunneling amplitudes and number of bands in the calculations, see text for details. The red marks at θ = 0 indicate the results of Ref. [7].
In Fig. 7(b) we show similar calculations for λ R vs θ. The larger tunneling strength t b used in this work leads again to an enhancement of λ R (thin solid line) compared to Ref. [1] (dashed line), but the θ dependence of the two results otherwise agree, including the deep minimum at θ = 30 • . However, if all possible e-o pairs of bands are included in the calculation (thick solid line), then this minimum becomes a small dip and the value of λ R is significantly enhanced. As already mentioned, this is because of band crossings and neardegeneracies of certain e bands which means that more pairs of e and o bands need to be taken into account than it was done in Ref. [1].
For comparison, in Fig. 7 we also indicate by red marks the results of the DFT calculations of Ref. [7] for θ = 0. As one can see, all our curves take values close to these reference values. As we explained in Sec. II B, the exact results depend on a number of details, including potentially fine-tuned Slater-Koster parameters. This leaves us with much freedom to adjust our parameters and fit previously published results, but this goes beyond the scope of this paper.
We found a qualitatively similar λ vZ and λ R vs θ dependence using the experimental parameters of E g and f G from Table II, but the maximum values of the induced SOC are significantly smaller.

Calculations for MLG/WSe2
Our calculations for the twist angle dependence of λ vZ and λ R are shown in Fig. 8. We used the tunneling parameters given in Tables I and the DFT parameters in Table II. Contributions from all bands or pairs of bands of the TB model are taken into account. Red marks denote the results of Ref. [7] for θ = 0, they differ ≈ 20% from our results.
FIG. 8: a) λ vZ and b) λ R for MLG/WSe 2 as a function of interlayer twist angle θ using the DFT parameters given in Table II. The red marks at θ = 0 indicate the results of Ref. [7]. Importantly, one can see that both λ vZ (θ) and λ R (θ) are different from the MLG/MoS 2 heterostructures in Fig. 7. This is mainly due to the fact that the Dirac point of graphene is closer to the VB (we used f G = 0.161 in these calculations versus f G = 0.976 in Fig. 7, i.e., here the valence bands give larger contributions. Similarly to MLG/MoS 2 , λ vZ and λ R vs θ dependence using the experimental parameters for E g and f G from Table II are qualitatively similar, but the maximum values of the induced SOC are significantly smaller. Comparison of MLG/MoS 2 and MLG/WSe 2 The different θ dependence of the induced SOC for MLG/MoS 2 (Fig. 7) and MLG/WSe 2 (Fig. 8) explain the differences in the results for the corresponding twisted trilayer systems shown in Fig. 3 of the main text. We point out, in particular, that λ R has a single maximum for MLG/MoS 2 , while it has two maxima for MLG/WSe 2 . Together with the θ dependence of ϑ R given in Fig. 1  3. Comparison to the results of Refs. [14] and [13] During the preparation of this manuscript, a publication by Naimer et. al. [14] has appeared discussing the twist angle dependence of the proximity-induced SOC in graphene/TMDC bilayers. Their methodology is based on DFT calculations. The incommensurate graphene/TMDC heterostructures were approximated by periodic supercells, which required several approximations, e.g., straining the graphene sheet.
In agreement with our results, the calculations of Ref. [14] indicate that the strength of the proximity induced SOC in graphene/TMDC heterostructures is tunable by interlayer twist. The twist angle dependence of λ vZ for graphene/MoS 2 and graphene/WSe 2 are qualitatively similar: in the former case there is a maximum at θ ≈ 20 • , in the latter case λ vZ monotonically decreases as θ is changed from 0 • to 30 • . They also find that the phase factor ϑ R in H R is in general non-zero and depends on the twist angle θ. Regarding the θ dependence of λ R , for graphene/MoS 2 their result seems to be qualitatively similar to ours, while for graphene/WSe 2 the results show differences. On a more quantitative level, however, our results for λ vZ and λ R display a much stronger enhancement as a function of θ than in Ref. [14]. The reasons for this discrepancy remains to be investigated in future works.
Another relevant recent work is by Pezo et. al. [13], who performed supercell based DFT calculations for graphene/MoTe 2 bilayers for three different interlayer twist angles (θ = 0, π/12, π/6). It was noted in Ref. [13] that the position of graphene's Dirac point is closer to the valence band of MoTe 2 . Therefore one would expect that the θ dependence of the induced SOC should be similar to the one in graphene/WSe 2 , meaning, e.g., that magnitude of both λ vZ and λ R should decrease for θ = π/6 compared to their θ = 0 value. This seems to be in agreement with results of Ref. [13].
Finally, we note that both Refs. [14] and [13] emphasize the effect of the strain that is always present in their supercell based DFT calculations.

III. SUPERMOIRÉ EFFECTS
Let us first consider the moiré effects in a graphene/TMDC system. One can define the primitive moiré wavelength λ m for a hexagonal heterostructure [15] consisting of graphene and a substrate as where q = a s /a gr is the ratio of the lattice constants of the substrate (a s ) and of graphene (a gr ) and θ is the interlayer twist angle. In the case of graphene and e.g., monolayer MoS 2 this ratio is q ≈ 1.283. The maximum of λ m can be found for θ = 0 where λ m ≈ 1.1 nm (see Fig.9). As a comparison, the moiré length scale is λ m = 13.9 nm for a graphene/hBN bilayer at θ = 0. The moiré potential due to this periodic perturbation leads to, e.g., gap openings at the moiré Brillouin zone boundaries [15] at energies Using v F = 10 6 m/s, for graphene/MoS 2 this energy scale is E M ≈ ±2.2eV for θ = 0, which is far from the Dirac point. Indeed, in the ARPES experimental results of Ref. [11] on graphene/MoS 2 , superlattice effects were observed at binding energies E b = −3.55eV, from which the authors concluded that the interlayer twist angle was θ ≈ 6 • in their sample. In short, because of the relatively large difference between the graphene and TMDC lattice constants, no moiré effects are to be expected close to the Dirac point of graphene for any twist angle θ in a graphene/TMDC system.
The situation might be different in TMDC/graphene/TMDC trilayers. Namely, studying hBN/graphene/hBN trilayers, Refs. 16 and 17 found that band gaps can appear in the spectrum of these heterostructures at energies closer to the Dirac point of graphene than what was found in hBN/graphene bilayers. This can be understood as a consequence of an interference of the moiré patterns coming from the two hBN layers. Ref. [16] found that if one of the interlayer twist angles, e.g., θ (b) is kept fixed at θ (b) = 0 and θ (t) is changed, then the four longest wavelength component λ (i) sm of such a supermoiré potential in graphene are given by where we used the notation q −1 = δ. Eqs. (39) are valid also for TMDC/graphene/TMDC trilayers because the constituent layers have hexagonal lattices. We plot λ m as well as λ (i) sm as a function of θ (t) using the lattice constants of graphene and MoS 2 in Fig. 9. As one can see, λ sm , λ sm , λ sm are always quite short. However, λ sm does increase when the two TMDC layers are nearly aligned. For θ (t) = 1 • one finds λ (1) sm = 18nm corresponding to E M = 0.13eV, which is comparable to the ones found in hBN/Graphene/hBN. Based on these results one may expect an interesting interplay of induced SOC and supermoiré effects in TMDC/graphene/TMDC for |θ (b) − θ (t) | 1 • .

IV. MAGNETOTRANSPORT CALCULATIONS FOR GRAPHENE npn JUNCTIONS
In this section we briefly explain the theoretical model that was used in the magnetotransport calculations through a npn graphene junction with proximity induced SOC shown in Fig. 4 of the main text. We assume that the graphene flake is perfectly ballistic and due to external gates it is doped such that it hosts an npn junction, see Fig. 10. We numerically calculate the transmission through the junction as a function of perpendicular magnetic field magnetic field B z and Fermi energy E F using the in-house code EQuUs [18].
where t(q n , E F , B z ) is the transmission coefficient for mode n. When there are many transverse momenta the exact form of the boundary conditions is not important and therefore we used the infinite mass boundary condition [21] to obtain q n : q n = n + 1 2 π W , where n = 0, 1, 2, . . . and W is the width of the junction.

V. ELECTRON SCATTERING AT A PLANAR JUNCTION IN A GRAPHENE-WSE2 SYSTEM
We will discuss a further setup to investigate twist-angle dependent transport in a graphene-TMDC system. It is motivated by recent electron optics experiment of Ref. [22,23] where collimated electron beams were created in monolayer graphene. By making use of such collimated electron beams one can test the twist-angle dependence of the induced SOC through spin-dependent scattering.
To demonstrate this, we consider the TMDC/MLG/TMDC stack shown in Fig. 11. It consists of a MLG/WSe 2 bilayer of a fixed twist angle θ (b) (purple region). Part of this bilayer is covered by a top WSe 2 layer (orange region) of twist angle θ (t) . Thus the induced SOC in graphene changes at the bilayer-trilayer boundary and in the trilayer region it is assumed to be tunable by θ (t) . The left edge of the top layer defines a junction along the y axis, where ballistic electrons scatter into various forward-propagating and back-propagating modes, as indicated in Fig. 11. We calculate below the reflection coefficient of this junction for a collimated beam of electrons. vZ (θ (t) ) in H gr vZ . One can easily show that the reflection coefficient is the same for electrons in both valleys, for concreteness, we perform the calculation for the K valley. We express the electronic states as plane waves multiplied by spinors in the basis {|A ↑ , |B ↑ , |A ↓ , |B ↓ } and we assume that their energy is E = 20 meV. Due to translational symmetry alongŷ axis, the δk y component of the wavevector of the electronic states is conserved. The SOC leads to the splitting of the bands, therefore (at the E F given above) for each direction of propagation one can find two modes for the incoming electrons with slightly different wave vectors and opposite pointing spins. We assume a net zero spin polarization with an equal mixture of the two incoming modes, and for each incoming mode, one can find two reflected and two refracted solutions. In total this gives four back-scattered and four forward-scattered modes within one valley. Depending on the exact parameters and the angle of incidence, some of these modes can be decaying waves carrying no current. The particle current operator