Spin-charge coupled transport in van der Waals systems with random tunneling

We study the electron and spin transport in a van der Waals system formed by one layer with strong spin-orbit coupling and a second layer without spin-orbit coupling, in the regime when the interlayer tunneling is random. We find that in the layer without intrinsic spin-orbit coupling spin-charge coupled transport can be induced by two distinct mechanisms. First, the gapless diffusion modes of the two isolated layers hybridize in the presence of tunneling, which constitutes a source of spin-charge coupled transport in the second layer. Second, the random tunneling introduces spin-orbit coupling in the effective disorder-averaged single-particle Hamiltonian of the second layer. This results in non-trivial spin transport and, for sufficiently strong tunneling, in spin-charge coupling. As an example, we consider a van der Waals system formed by a two-dimensional electron gas (2DEG)--such as graphene--and the surface of a topological insulator (TI) and show that the proximity of the TI induces a coupling of the spin and charge transport in the 2DEG. In addition, we show that such coupling can be tuned by varying the doping of the TI's surface. We then obtain, for a simple geometry, the current-induced non-equilibrium spin accumulation (Edelstein effect) caused in the 2DEG by the coupling of charge and spin transport.


I. INTRODUCTION
In recent years experimentalists have been able to make very novel and high-quality heterostructures that allow the realization of new effects and states of great fundamental and technological interest [1]. Recently simple heterostructures formed by two graphene layers with a relative twist angle [2][3][4][5] have shown a phase diagram [6,7] that is remarkably reminiscent of the phase diagram of high-temperature superconductors. These are just some of the most striking examples that heterostructures can be used to realize novel effects that are not present in the single constituents. Applications of heterostructure engineering [8] can be found in tunnel junctions [9], plasmonic [10], photoresponsive [11], spintronics [12][13][14] and valleytronic [15] devices.
One of the essential elements to realize nontrivial topological states and spin-dependent transport phenomena is the presence of spin-orbit coupling (SOC). However, often the presence of spin-orbit coupling is not accompanied by other desirable properties such as high mobility, or superconducting pairing. For this reason heterostructures that combine one constituent with significant SOC and one constituent with no SOC but other distinct properties are very interesting both for fundamental reasons and for their potential for technological applications. So far, the theoretical studies of van der Waals heterostructures have focused on the regime when the tunneling is not random and a strong hybridization between the electronic states of the isolated systems can be achieved, see for example the case of graphene-topological-insulator systems [16][17][18][19][20][21][22][23][24][25]. However, in many situations we can expect the tunneling between the systems forming the heterostructure to be random, due for example to the incommensurate nature of the stacking configuration and/or the presence of surface roughness.
In this work we focus on this situation and study the electron and spin transport in a two-dimensional van der Waals system comprised of one component (layer) with strong SOC and one with no, or negligible, SOC, when the interlayer tunneling is random. Due to the random nature of the tunneling in most experimental situations the transport will be diffusive even in the absence of disorder. For this reason we consider only the diffusive regime, in which specific details of the system considered (like the value of the mean free path) do not affect the general expression of the transport equations that, therefore, have a somewhat universal character. We find that in general, if the diffusive transport in the layer with SOC exhibits spin-charge coupling [26][27][28] such coupling will be present also in the layer without SOC, i.e., in the most common experimental situation. To exemplify this general result we consider the case of a van der Waals system formed by a two-dimensional electron gas (2DEG) placed on the surface of a strong three-dimensional topological insulator (TI) [29][30][31]. Graphene and the surface of TIs in the tetradymite family such as Bi 2 Se 3 have almost commensurate lattices and as a consequence in many graphene-TI heterostructures the K, K points of the graphene's Brillouin zone are folded close to the TI's point [19]. This fact, combined with the random and finite-range nature of the interlayer tunneling, implies that the results that we obtain for a 2DEG-TI van der Waals system are directly relevant to graphene-TI heterostructures. Besides graphene, other 2D systems with low energy states around the point could be candidates for hosting the 2DEG [32][33][34]. We obtain the diffusive transport equations in the 2DEG layer and show that they describe a transport in which the charge and the spin degrees of freedom are coupled. Finally, we show how the diffusive equations give rise to spindependent transport effects, analogous to the ones obtained for a 2DEG with Rashba SOC [26] and an isolated TI's surface [27], that are tunable by simply varying the doping of the TI, and that can be used for possible spintronics applications.
The rest of this work has the following organization. In Sec. II, we define the heterostructure Hamiltonian and the disorder and tunneling potentials, and we calculate the relevant response function. In Sec. III, we derive the spin-charge coupled transport equations. In Sec. IV, we solve the diffusion equations for a particular experimental setup and compare the results with those for a 2DEG with spin-orbit coupling and a TI alone. Finally, in Sec. V we present our conclusions.

II. MICROSCOPIC APPROACH
In this section, we introduce the heterostructure Hamiltonian along with the impurity and random tunneling potentials. Then, we calculate the self-energies and relevant response function, including vertex corrections. Finally, from the response function, we extract the inverse diffuson.
The HamiltonianĤ for the heterostructure can be written asĤ where l is the layer index,Ĥ l is the Hamiltonian for layer l in the clean limit,V l is the term due to disorder located in layer l, andT is the term describing interlayer tunneling. For the 2DEG layer we haveĤ l =Ĥ 2d (k) = kss ψ † 2d,ks H 2dss (k)ψ 2d,ks , whereψ † 2d,ks (ψ 2d,ks ) is the creation (annihilation) operator for an electron with momentum k and spin s. Without loss of generality we can linearize the 2DEG dispersion around the Fermi surface and assume H 2d (k) = (v 2d |k| − μ 2d )σ 0 with v 2d the Fermi velocity, μ 2d the chemical potential, and σ 0 the 2 × 2 identity Pauli matrix in spin space. For the TI's surface we haveĤ l =Ĥ TI = kss ψ † TI,ks H TIss (k)ψ TI,ks , where ψ † TI,ks (ψ † TI,ks ) creates (annihilates) a surface Dirac fermion with spin s and momentum k, H TI (k) = −v TI (k × σ ) z − μ TI , v TI being the Fermi velocity on the TI's surface, μ TI the TI's surface chemical potential, and σ i , i = x, y the Pauli matrices in spin space.
For the disorder potential in layer l, V (D) , where the angle brackets denote average over disorder realizations, and W D l (r 1 − r 2 ) is the disorder-averaged spatial correlation. In momentum space we have W D l (q) = n l imp |U (q)| 2 , where n l imp is the impurity density in layer l, and U l (q) the Fourier transform of the potential profile U l (r) of a single impurity. Without loss of generality we can set V (D) l (r) = 0. Assuming the tunneling to be spin-conserving we havê wherel = l, and t 0 and T (q) are the spatially uniform and random components of the tunneling amplitude, respectively. We consider the limit where the uniform component is negligible and the remaining random component can be characterized by the spatial average of the tunneling matrix element T (r 1 )T (r 2 ) = W t (r 1 − r 2 ). This regime could arise, for example, due to random variations in the interlayer distance or in the alignment of orbitals on the surface. In the remainder we assume both the intralayer disorder and interlayer tunneling to be short range so that U l (q) = const = U l , and be the bare retarded (advanced) real-time Green's function for layer l. The total self-energy for layer l, l , has contributions from scattering with impurities 0 l , and random tunneling events t l . We have where q ≡ d 2 q/(2π ) 2 . In the self-consistent Born approximation, G l is the disorder-dressed Green's function for layer l. For the 2DEG, apart from an overall unimportant real constant, we have 0 imp U 2 2d , and ρ 2d is the density of states (DOS) at the Fermi energy. For the TI's surface, because the electrons behave as massless Dirac fermions, for U TI (q) = const, we have that the integral in the expression for 0 l has an ultraviolet divergence [35]. After properly regularizing such divergence [36] one finds that the intralayer disorder, in addition to generating an imaginary part of the self-energy, −i 0 TI σ 0 , with 0 TI = 1/τ 0 TI = πρ TI n TI imp U 2 TI and ρ TI , the TI's DOS at the Fermi energy, causes a renormalization of the Fermi velocity that we incorporate in the definition of v TI . The same ultraviolet divergence appears for the self-energy correction for the 2DEG due to tunneling events into the TI, t 2d . The proper renormalization of such divergence, consistent with the Ward identities, causes t 2d to have a nontrivial real part so that This result shows that even when the interlayer tunneling processes are random, a spinorbit coupling term is induced in the 2DEG due to TI's surface proximity. This term of the self-energy qualitatively affects the diffusive transport in the 2DEG, but it is not necessary to induce spin-charge transport in the 2DEG as we will show below. The self-energy correction for the TI due to tunneling events into the 2DEG, t TI , does not require any special care and simply results in an additional broadening of the quasipar- With the self-energy contributions, the dressed 2D system Green's functions take the form where 2d ≡ 0 2d + t 2d and TI ≡ 0 TI + t TI . In the diffusive regime, to leading order in 1/( F τ ), the retarded dynamical part of the spin-density response function for layer l, χ dyn l is obtained by summing all ladder vertex corrections to the bare spin-density response. In our case we have two types of ladder diagrams: the ones due to random interlayer tunneling and the ones due to intralayer disorder. In most experimentally relevant situations we expect the scattering time due to intralayer disorder to be much smaller than the relaxation time due to the interlayer random tunneling processes. For this reason in the remainder we assume t 0 . The main building block for the calculation of χ dyn 2d is the diffuson D 2d , which includes both interlayer tunneling and intralayer ladder diagrams. It satisfies the self-consistent equation [37,38] In this equation, the auxiliary intralayer diffuson for layer l,D l [l = (2d, TI)] is obtained by taking into account only intralayer disorder and the junctions J describe the transition between the layers. The constant κ collects disorderdependent normalizations with κ −1 = n 2d imp n TI imp U 2 2d U 2 TI . The self-consistency equation (7) is shown diagrammatically in Fig. 1(a).
Mathematically, the auxiliary diffusonD l satisfies the Bethe-Salpeter equation [see Fig. 1 Here, the quantum probability P l is defined as The junctions J l l = P l t 2 P l account for the tunneling processes. The expressions of P 2d and P TI are given in Appendix A.
For the purpose of finding χ dyn l , it is convenient to solve Eq. (7) in the spin-charge representation. To this end the diffusons, as well as the junctions, are contracted with the Pauli matrices as D αβ x, y, z) correspond to the charge and x, y, z components of the spin, respectively. With the knowledge of D l , the dynamical part of the spin-density response function can be found by introducing charge and spin vertices as illustrated in Fig. 1(c). The full response function is then obtained by adding the static part, χ l = χ st l + χ dyn l , where χ st,αβ l ∝ ρ l δ αβ . For systems with conserved particle number, the density response function χ 00 must satisfy the condition lim ω→0 [lim q→0 χ 00 (q, ω)] = 0. In the problem under consideration, electrons can move from one layer to the other. Therefore, a complete description of the time evolution of the charge and spin densities must include the mixed response function χ ll with l = l , i.e., the response of densities in layer l to perturbations in layer l . χ ll can be found in analogy to χ l .

III. DIFFUSION EQUATIONS
In section, we derive the charge and spin coupled diffusion equations and discuss the main general implications. In the 2DEG, the charge and spin response to external perturbations in the form of electric potentials or Zeeman fields may be conveniently cast in the form of coupled transport equations. In the diffusive limit, we find ∂ t n 2d =D∇ 2ñ 2d + ns l TI (ẑ × ∇)s 2d − ν∂ t (V 2d − V TI ), (10) where the effective charge diffusion constant is a weighted average of the diffusion constants D 2d = v 2 F τ 0 2d /2, and D TI = v 2 TI τ 0 TI in the 2DEG and TI, respectively. Moreover, l TI = v TI τ 0 TI is the TI mean free path. The spincharge coupling in the 2DEG is characterized by ns = 2 t 2d t TI /( t 2d + t TI ). The term containing the dimensionless constant α = F τ 0 2d /(2π 2 ρ TI D TI ) originates from the induced spin-orbit coupling in the 2DEG. The charge and spin densitiesñ ands appearing on the right-hand side of the diffusion equations include external driving potentials for the charge, V 2d , and spin, h 2d , respectively, asñ 2d = n 2d + 2ρ 2d V 2d and s = s − 2ρ 2d h 2d . The last term in Eq. (10) accounts for a potential loss of electrons in the 2DEG for a dynamically driven system, with coefficient ν = 2ρ 2d t 2d /( t 2d + t TI ). Equations (10)- (12) are the main result of this work. They show that in a 2DEG-TI system charge transport and spin transport are coupled even when the tunneling between the two systems is random. Notice that Eqs. (10)-(12) were obtained in the limit in which t l / 0 l 1, and ωτ 1, τ being the longest relaxation time: τ = max(τ t 2d , τ t TI ). Equations (10)-(12) can only describe transport over time scales much 033085-3 larger than τ and therefore are not valid in the limit t = 0 for which τ → ∞. For t = 0 the two systems are decoupled and for the 2DEG the diffusive transport of charge and spin are independent with D 2d = v 2 2d τ 0 2d /2. It is instructive to note that there are two mechanisms responsible for the spin-charge and spin-spin coupling in Eqs. (10) and (11). The term with coefficient α in Eq. (11) results from the real part of the self-energy in Eq. (4), i.e., from the tunneling-induced spin-orbit coupling in the effective single-particle Hamiltonian of the 2DEG. This term couples in-plane and out-of-plane spin components. The spin-charge coupling in Eqs. (10) and (11) has a different origin. The surface of the TI hosts a single gapless diffusion mode in the absence of tunneling, as can be seen by diagonalizing the diffuson [27,37,38]. For finite q, this mode has a nontrivial spin structure. By means of the random tunneling, this mode and the gapless modes in the 2DEG hybridize. The hybridization gives rise to spin-charge coupling via the term with coefficient ns in Eq. (10) and the final term in Eq. (11), as well as to anisotropic spin-diffusion encoded in the first term of the second line in Eq. (11). To leading order in tunneling, the two described mechanisms for spin-charge coupling are independent of each other. As follows from Ref. [26], spinorbit coupling eventually also leads to spin-charge coupled transport at higher orders in the coupling strength. A separate consequence of the tunneling in Eq. (11) is that, since spin is not conserved in the coupled system, a gap of size t 2d opens for the spin diffusion modes.
Equations (10) and (11) show that the strength of the coupling between charge transport and spin transport, and the spin-diffusion anisotropy, are proportional to the ratio t 2d / 0 2d . Given that t 2d = t 2 ρ TI π and that ρ TI scales linearly with μ TI , we see that in the 2DEG both the spin-charge coupling and the spin-diffusion anisotropy can be tuned simply by changing the doping of the TI's surface.

IV. APPLICATIONS
We now study the solution of Eqs. (10) and (11) for a simple setup, as in Refs. [26,27], to highlight some of the transports effects due to the coupling between spin and charge transport described by Eqs. (10) and (11), and to highlight some of the main similarities and differences between a 2DEG-TI system, a TI's surface, and a 2DEG with Rashba SOC.
We consider a system of size L along x, −L/2 < x < L/2, and in which all the quantities are uniform along y. In the stationary limit, due to the uniformity along y, Eqs. (10) and (11) separate into two independent sets of equations: one set describing the coupled transport of n and s y , one set describing the coupled transport of s x and s z . Given that we are interested in the coupling between charge and spin transport, we focus on the first set. Due to the assumption that all the quantities are homogenous along y, the coupled equations for n and s y for a 2DEG-TI, a TI, and a 2DEG with Rashba SOC have the same structure: where D n , D s , β n , and β s are constants whose expressions in terms of the parameters characterizing the system are given in Table I for a 2DEG-TI, a TI, and a 2DEG with Rashba SOC. From charge conservation, using Eq. (10), we find that the charge current takes the form and for the simple case described by Eq. (13), J = Jx, J = −D n dn/dx + 2β s s y , with D n and β s given in Table I. Similarly from Eq. (14) we can obtain an expression for the current of s y . This expression has the term β n ∂ x n; however, as pointed out before [39][40][41][42][43], such a term describes an equilibrium spin current and therefore should not be included in the definition of an externally driven spin current. Knowing the expression of J and of the spin current allows us to write the boundary conditions for Eqs. (13) and (14), corresponding to the situation when a charge current I is injected at x = −L/2 via a ferromagnetic electrode so that the incoming electrons have a net spin polarization φ along s y : Recalling that the voltage drop V (x) [44] at position x is given by V (x) = −(1/2eρ) x −L/2 dx (dn/dx ), and solving Eqs. (13) and (14) with the boundary conditions (16) we find (17) and the voltage drop between the leads In Eqs. (17) and (18), l −2 * ≡ 1/(τ s D s ) + 2β n β s /(D n D s ). Using the expressions given in Table I for D n , D s , τ s , β n , and β s , Eqs. (17) and (18) for a 2D-TI system become, to leading order in the tunneling amplitude (with l * ≈ D 2d τ t 2d ), The second term on the right-hand side of Eq. (19) shows that, as in the case of 2DEG with Rashba SOC [26] and a TI [27], an Edelstein [45] effect is present, i.e., a constant 033085-4 nonequilibrium spin polarization generated by a charge current I. This effect is present due to the "mirroring" into the 2DEG of the TI's gapless diffusion mode characterized by the coupling of charge and spin. It is interesting to notice that for a 2DEG-TI system such a term, as long as τ t 2d τ 0 2d to remain in the regime of validity of the diffusion equations (10) and (11), is independent of the interlayer tunneling strength. This is because in the 2DEG-TI van der Waals structure, in the 2DEG layer, both the spin relaxation rate 1/τ s and the spin-charge coupling β n in Eq. (14) scale as t 2 . As a consequence we expect that even in the limit of very small t a significant Edelstein effect should be present in a metallic 2D layer placed in proximity of a system with significant SOC such as a TI's surface. In addition, we see that for a 2DEG-TI system, contrary to a TI, the strength of the Edelstein effect can be tuned by varying the doping, and therefore ρ TI , of the TI's surface. The other important result is that the decay length of s y is l * which can also be tuned by varying the doping in the TI and can be very long in the weak tunneling regime, for which τ t 2d τ 0 2d . The last term on the right-hand side of Eq. (20) is a magnetoresistance contribution to the voltage drop due to the coupling of the charge and spin transport. For a 2DEG-TI system this term is therefore dependent on the relative strength of the disorder in the TI and 2DEG.

V. CONCLUSIONS
In conclusion, we have studied the electron and spin transport in a van der Waals system formed by one layer with strong spin-orbit coupling and a second layer without spinorbit coupling, in the regime when the interlayer tunneling is random. We have shown that in the layer without intrinsic spin-orbit coupling, spin-charge coupled transport can be induced by the hybridization of the diffusion modes of the two isolated layers. To exemplify the mechanism we have studied a van der Waals system formed by a 2DEG and TI's surface and shown how the coupling of the spin and charge transport in the TI is "mirrored" into the 2DEG. In addition, for the specific case of a 2DEG-TI van der Waals system, we show that a spin-orbit coupling term is induced into the 2DEG, and that the induced coupling of spin and charge transport in the 2DEG can be tuned by varying the TI's doping. Finally we showed how the coupled spin-charge transport described by the diffusive equations that we obtain for the 2DEG leads to a current-induced nonequilibrium spin accumulation and a magnetoresistance effect that are also tunable by changing the TI's doping.

APPENDIX B: SPIN-CHARGE DIFFUSION EQUATION FOR TI's SURFACE
To facilitate the comparison between the results that we obtain in the main text for a 2DEG-TI system and an isolated TI's surface we report here the diffusion equations for a TI's surface, first derived in Ref. [27]: where n TI is the carrier density on the TI's surface, and s TI = (s x TI , s y TI ). Notice that the spin densities are damped by scattering with nonmagnetic impurities due to spin-orbit coupling. Due to a typographic error in Ref. [27] the terms with mix derivatives have opposite sign compared to Eqs. (B2) and (B2). We can see that the negative sign in front of the terms ∂ 2 xy s x TI and ∂ 2 xy s y TI in Eqs. (B2) and (B2) is correct by considering that when n TI is uniform in time and space so that Eq. (B1) implies ∂ y s x TI = ∂ x s y TI , Eqs. (B2) and (B3) lead to ∂ t s α TI = [(1/2)D TI ∇ 2 − 1/τ 0 TI ]s α TI , the expected spin-diffusion equation in this simple limit.

APPENDIX C: DIFFUSION EQUATIONS FOR TWO COUPLED 2DEGs
In this Appendix, we review the density diffusion equation of a 2DEG-2DEG heterostructure. The effect of the coupling in the quantum interference has been studied before [46]. Each layer l poses its own diffusion constant D l and density of states ρ l , where l = T, B label the top and bottom 2DEG layers, respectively. We obtain ∂ t n (T ) 2d =D∇ 2 n (T ) 2d , where we have defined The renormalized diffusion constant contains corrections proportional to the diffusion constant in the bottom layer. The leading correction to the diffusion constant is given by a term proportional the ratio of the DOS in each layer. Given that there is no spin-orbit coupling, the spin follows analogous diffusion equations in each direction.