Twist-angle dependence of the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides

TL;DR

This paper theoretically investigates how the twist angle between graphene and transition-metal dichalcogenides affects proximity-induced spin-orbit coupling, revealing significant enhancement and tunability of spin splitting with potential for twist-angle engineering.

Contribution

It introduces a theoretical method applicable to incommensurate bilayers, demonstrating twist-angle control over spin-orbit coupling in graphene on TMDCs.

Findings

Spin splitting is enhanced by a factor of 3-10 with twist angles.

Maximum spin splitting occurs around 20° twist angle.

The type of SOC can be switched from Zeeman to Rashba by rotation.

Abstract

We theoretically study the proximity spin-orbit coupling in graphene on transition-metal dichalcogenides monolayer stacked with arbitrary twist angles. We find that the relative rotation greatly enhances the spin splitting of graphene, typically by a few to ten times compared to the non-rotated geometry,and the maximum splitting is achieved around $20^\circ$. The induced SOC can be changed from the Zeeman-type to the Rashba-type by rotation. The spin-splitting is also quite sensitive to the gate-induced potential, and it sharply rises when the graphene's Dirac point is shifted toward the TMDC band. The theoretical method does not need the exact lattice matching and it is applicable to any incommensurate bilayer systems. It is useful for the twist-angle engineering of a variety of van der Waals proximity effects.

Figures

Tables

Table 1. Table 1: List of parameters for TMDCs and graphene-TMDC bilayers used in this work (see the text).

	MoS2	MoSe2	WS2	WSe2
$a_T$(Å)	3.18Fang et al. (2015)	3.32Fang et al. (2015)	3.18Fang et al. (2015)	3.32Fang et al. (2015)
$w$(Å)	3.13Fang et al. (2015)	3.34Fang et al. (2015)	3.14Fang et al. (2015)	3.35Fang et al. (2015)
$d$(Å)	3.37Fang et al. (2015)	3.41Ma et al. (2011)	3.41Kaloni et al. (2014)	3.42Kaloni et al. (2014)
$E_T-E_G$(eV)	0.02Gmitra and Fabian (2015)	0.6Ma et al. (2011)	0.12Kaloni et al. (2014)	1.06Agnoli et al. (2018); Kaloni et al. (2014)

Equations

$$\displaystyle{\bf R}_{X}=n_{1}{\bf a}_{1}+n_{2}{\bf a}_{2}+\mbox{\boldmath$\tau$}_{X}\quad(\mbox{graphene}),$$

$$\displaystyle{\bf R}_{\tilde{X}}=\tilde{n}_{1}\tilde{{\bf a}}_{1}+\tilde{n}_{2}\tilde{{\bf a}}_{2}+\mbox{\boldmath$\tau$}_{\tilde{X}}\quad(\mbox{TMDC}),$$

$$|\mathbf{k},X,s⟩=\frac{1}{\sqrt{N}}\sum _{{\mathbf{R}}_X}{e}^{i\mathbf{k}\cdot {\mathbf{R}}_X}|{\mathbf{R}}_X,s⟩(\text{mbox}graphene),$$

$$|\tilde{\mathbf{k}},\tilde{X},\tilde{s}⟩=\frac{1}{\sqrt{\tilde{N}}}\sum _{{\mathbf{R}}_{\tilde{X}}}{e}^{i\tilde{\mathbf{k}}\cdot {\mathbf{R}}_{\tilde{X}}}|{\mathbf{R}}_{\tilde{X}},\tilde{s}⟩(\text{mbox}TMDC),$$

$$⟨\tilde{\mathbf{k}},\tilde{X},\tilde{s}⟩H_{\mathrm{int}}\mathbf{k},X,s=$$

$$-\sum _{\mathbf{G},\tilde{\mathbf{G}}}{t}_{\tilde{X}X}(\mathbf{k}+\mathbf{G})e^{-i\mathbf{G}\cdot {\mathbf{\tau }}_X+i\tilde{\mathbf{G}}\cdot {\mathbf{\tau }}_{\tilde{X}}}{\delta }_{\mathbf{k}+\mathbf{G},\tilde{\mathbf{k}}+\tilde{\mathbf{G}}}{\delta }_{\tilde{s}s}.$$

$$t_{\tilde{X}X}(\mathbf{q})=\frac{1}{\sqrt{S\tilde{S}}}\int T_{\tilde{X}X}(\mathbf{r}+z_{\tilde{X}X}{\mathbf{e}}_z)e^{-i\mathbf{q}\cdot \mathbf{r}}{d}^{2}r,$$

$$[V_{\mathrm{eff}}(\mathbf{k}){]}_{X^{\prime }{s}^{\prime },Xs}=$$

$$\sum _{\tilde{n},\tilde{\mathbf{k}}}\frac{⟨\mathbf{k},X^{\prime },s^{\prime }⟩H_{\mathrm{int}}\tilde{n},\tilde{\mathbf{k}}⟨\tilde{n},\tilde{\mathbf{k}}⟩H_{\mathrm{int}}\mathbf{k},X,s}{E_G-E_{\tilde{n},\tilde{\mathbf{k}}}}.$$

$$V_{\mathrm{eff}}^{(\xi )}=\frac{\lambda }{2}\xi s_z+\frac{{\lambda }_R}{2}{e}^{-i\varphi s_z/2}(\xi {\sigma }_x{s}_y-{\sigma }_y{s}_x)e^{i\varphi s_z/2},$$

$$\displaystyle V^{(+)}_{\rm eff}=\left(\begin{array}[]{cccc}\lambda/2&&&\\ &\lambda/2&-i\lambda_{R}e^{-i\phi}&\\ &i\lambda_{R}e^{i\phi}&-\lambda/2&\\ &&&-\lambda/2\end{array}\right),$$

$$\displaystyle V^{(-)}_{\rm eff}=\left(\begin{array}[]{cccc}-\lambda/2&&&i\lambda_{R}e^{-i\phi}\\ &-\lambda/2&&\\ &&\lambda/2&\\ -i\lambda_{R}e^{i\phi}&&&\lambda/2\end{array}\right),$$

$${\sigma }_{xy}^{(\xi )}=\frac{e^2}{h}\sum _{n\in occ.}\int \frac{d^2\mathbf{k}}{2\pi }{\nabla }_{\mathbf{k}}\times {\mathbf{a}}_n(\mathbf{k})=-\frac{e^2}{h}\xi \mathrm{sgn}(\lambda ),$$

Full text

††preprint: APS/123-QED

The physical properties of 2D material are generally sensitive to the interference with other materials placed in contact. In recent years, a great deal of experimental and theoretical efforts have been made to explore the proximity-induced phenomena in van der Waals heterostructures consisting of different 2D crystals.Geim and Grigorieva (2013) In particular, it was shown that the negligibly small spin-orbit coupling (SOC) of intrinsic graphene can be significantly enhanced by superimposing on the surface of transition-metal dichalcogenides (TMDC), Avsar et al. (2014); Kaloni et al. (2014); Gmitra and Fabian (2015); Wang et al. (2015a) and it is expected to be useful to realize spintronic manipulation in graphene.

In the studies on the proximity effect on 2D materials, however, the importance of the relative lattice orientation has often been overlooked. The previous theoretical calculations of proximity spin-orbit effects of graphene/TMDC system are limited to the non-rotated geometry. Kaloni et al. (2014); Gmitra and Fabian (2015); Wang et al. (2015a) On the other hand, the sensitive dependence on the relative twist angle $\theta$ was noticed in various 2D hetrostructures, and controlling $\theta$ is expected to be powerful means of manipulating their electronic properties. Carr et al. (2017); Ribeiro-Palau et al. (2018) In graphene on hexagonal BN system, for instance, the moiré interference pattern gives rise to the formation of the secondary Dirac points and the miniband structure. Kindermann et al. (2012); Wallbank et al. (2013); Mucha-Kruczyński et al. (2013); Jung et al. (2014); Moon and Koshino (2014); Dean et al. (2013); Ponomarenko et al. (2013); Hunt et al. (2013); Yu et al. (2014) The twisted bilayer graphene also exhibits the dramatic angle-dependent phenomena, such as the flat band formation Lopes dos Santos et al. (2007); Mele (2010); Trambly de Laissardiere et al. (2010); Shallcross et al. (2010); Morell et al. (2010); Bistritzer and MacDonald (2011); Moon and Koshino (2012); de Laissardiere et al. (2012) and the emergent superconductivity. Cao et al. (2018a, b) For graphene/TMDC hetrostructure, the twist-angle dependent band structure was theoretically simulated for several commensurate angles by the density functional theory (DFT) Wang et al. (2015b); Di Felice et al. (2017); Hou et al. (2017), and it is also experimentally probed. Jin et al. (2015); Pierucci et al. (2016); Du and Yu (2017); Lu et al. (2017) However, the $\theta$-dependence of spin-orbit coupling induced on graphene remains still unclear. It is generally hard to consider arbitrary twist angles in the DFT calculation, because it requires exact lattice matching to have a finite unit cell.

In this letter, we theoretically study the proximity SOC effect in graphene-TMDC heterostructures with arbitrary twist angles $\theta$, and reveal the angle dependence of SOC for various different TMDCs. Using the tight-binding model and the perturbational approach, which do not need the commensurate lattice matching, we obtain the effective Hamiltonian of graphene as a continous function of $\theta$. We find that the relative rotation greatly enhances the spin splitting, typically by a few to ten times compared to the non-rotated geometry $(\theta=0)$, and the maximum splitting is achieved around $\theta\sim 20^{\circ}$. We also show that the induced SOC is composed of the Zeeman-like term and the rotated Rashba-like term, and the relative magnitude can be controlled by rotation. Finally, we demonstrate that the spin-splitting is quite sensitive to the relative band energy between graphene and TMDC, and it sharply rises when the graphene’s Dirac point is shifted toward TMDC band by applying the gate voltage. The theoretical method proposed here is applicable to any incommensurate bilayer systems where the DFT calculation cannot be used, and therefore it considerably extends the applicability of the theoretical framework to a wide variety of van der Waals heterostructures.

We consider monolayer graphene placed on the top of a TMDC monolayer. Graphene and TMDC are two dimensional honeycomb lattices with different lattice periods, $a_{G}=2.46$Å  for graphene and $a_{T}$ for TMDC given in table 1. We define the stacking geometry starting from non-rotated arrangement with parallel bond directions, and then rotating TMDC by the twist angle $\theta$ around the common center of hexagon as in Fig. 1(a). The lattice structure has the $C_{3}$ (120∘) rotational symmetry with respect to the rotation center. We neglect the degree of freedom of the in-plane parallel translation between TMDC and graphene, since in an incommensurate system it can always be incorporated with the shift of the origin. Koshino (2015)

The lattice vectors of graphene are then given by $\mathbf{a}_{1}=a_{G}(1,0)$ and $\mathbf{a}_{2}=a_{G}(1/2,\sqrt{3}/2)$, and those of TMDC are by $\tilde{\mathbf{a}}_{1}=Ra_{T}(1,0)$ and $\tilde{\mathbf{a}}_{2}=Ra_{T}(1/2,\sqrt{3}/2)$, where $R=R(\theta)$ is the rotation matrix. The unit cell area is $S=|\mathbf{a}_{1}\times\mathbf{a}_{2}|$ and $\tilde{S}=|\tilde{\mathbf{a}}_{1}\times\tilde{\mathbf{a}}_{2}|$ for graphene and TMDC, respectively. The reciprocal lattice vectors, $\mathbf{a}^{*}_{1},\mathbf{a}^{*}_{2},\tilde{\mathbf{a}}^{*}_{1},\tilde{\mathbf{a}}^{*}_{2}$, are defined by $\mathbf{a}_{i}\cdot\mathbf{a}^{*}_{j}=\tilde{\mathbf{a}}_{i}\cdot\tilde{\mathbf{a}}^{*}_{j}=2\pi\delta_{ij}$. We define $d$ as the distance between the graphene layer and the top chalcogen layer, and $w$ as the distance between top and bottom chalcogen layers. The values of $d$ and $w$ depend on TMDCs as shown in table 1.

We model graphene by the tight-binding model of carbon $p_{z}$ orbitals, where the sublattice is labeled as $X=p_{z}^{A},p_{z}^{B}$ for A and B sites, respectively. For TMDC, we adopt the tight-binding model including three $p$ orbitals for a chalcogen atom and five $d$ orbitals for a transition metal atom Fang et al. (2015). The orbitals in a TMDC unit cell is labeled by $\tilde{X}=d_{z^{2}},d_{xy},d_{x^{2}-y^{2}},d_{xz},d_{yz},p_{x}^{t},p_{y}^{t},p_{z}^{t},p_{x}^{b},p_{y}^{b},p_{z}^{b}$, where $t$ and $b$ represent top and bottom chalcogen layers. The positions of the orbitals are given by

[TABLE]

where $n_{i}$ and $\tilde{n}_{i}$ are integers, and \mbox{\boldmath\tau}_{X} and \mbox{\boldmath\tau}_{\tilde{X}} are the sublattice position inside the unit cell. Specifically, they are expressed as $\bm{\tau}_{p_{z}^{A}}=-\bm{\tau}_{1}$, $\bm{\tau}_{p_{z}^{B}}=\bm{\tau}_{1}$ for graphene, and $\bm{\tau}_{\tilde{X}}=-\tilde{\bm{\tau}}_{1}-(d+w/2){\bf e}_{z}$ for the transition metal $d$-orbitals and $\bm{\tau}_{\tilde{X}}=\tilde{\bm{\tau}}_{1}-d{\bf e}_{z},\tilde{\bm{\tau}}_{1}-(d+w){\bf e}_{z}$ for the top and bottom charcogen $p$-orbitals, respectively, where $\bm{\tau}_{1}=(-\mathbf{a}_{1}+2\mathbf{a}_{2})/3$ and $\tilde{\bm{\tau}}_{1}=(-\tilde{\mathbf{a}}_{1}+2\tilde{\mathbf{a}}_{2})/3$.

The Hamiltonian is spanned by the Bloch bases,

[TABLE]

where $s,\tilde{s}$ are the spin indexes, ${\bf k}$ and $\tilde{{\bf k}}$ are the two-dimensional Bloch wave vectors parallel to the layer, and $N=S_{\rm tot}/S$ and $\tilde{N}=S_{\rm tot}/\tilde{S}$ are the number of unit cells of TMDC and graphene, respectively, in the total system area $S_{\rm tot}$.

The total tight-binding Hamiltonian is expressed as $H=H_{\rm G}+H_{\rm T}+H_{\rm int},$ where $H_{\rm G}$ and $H_{\rm T}$ are the Hamiltonian for the intrinsic graphene monolayer and TMDC monolayer, respectively, and $H_{\rm int}$ is for the coupling between graphene and TMDC. For $H_{\rm T}$, we adopt the hopping parameters based on the first principles calculation Fang et al. (2015) where the spin-orbit coupling is included by on-site ${\bf L}\cdot{\bf S}$ term for each atom. The on-site energy of the TMDC atoms relative to the carbon atoms is extracted from the relative energy $E_{T}-E_{G}$ from the graphene Dirac point to TMDC conduction band edges in the first principles calculationsGmitra and Fabian (2015); Agnoli et al. (2018); Ma et al. (2011, 2011); Kaloni et al. (2014), which are listed in Table 1.

For the interlayer interaction, we assume that the transfer integral from ${\bf R}_{X}$ to ${\bf R}_{\tilde{X}}$ is expressed as $-T_{\tilde{X}X}({\bf R}_{\tilde{X}}-{\bf R}_{X})$, with the Slater-Koster parameterization Slater and Koster (1954) and the exponential decay in the distance. Here the hopping amplitude and the decay length are determined to fit the first principles calculations. The detailed method is described in the supplementary materials. The coupling between the Bloch state of graphene and that of TMDC is then given by Bistritzer and MacDonald (2011); Koshino (2015)

[TABLE]

Here ${\bf G}=m_{1}{\bf a}^{*}_{1}+m_{2}{\bf a}^{*}_{2}$ and $\tilde{{\bf G}}=\tilde{m}_{1}\tilde{{\bf a}}^{*}_{1}+\tilde{m}_{2}\tilde{{\bf a}}^{*}_{2}$ are reciprocal lattice vectors of graphene and TMDC, respectively, ${t}_{\tilde{X}X}({\bf q})$ is the in-plane Fourier transform of the transfer integral defined by

[TABLE]

where z_{\tilde{X}X}=(\mbox{\boldmath\tau}_{\tilde{X}}-\mbox{\boldmath\tau}_{X})\cdot{\bf e}_{z}.

The Hamiltonian of graphene including the TMDC proximity effect can be obtained by the second order perturbation as $H_{\rm eff}({\bf k})=H_{G}({\bf k})+V_{\rm eff}({\bf k})$, where

[TABLE]

Here $E_{G}$ is the energy of the graphene’s Dirac point, and $E_{\tilde{n},\tilde{{\bf k}}}$ and $|\tilde{n},\tilde{{\bf k}}\rangle$ are the eigen energy and eigen state of $H_{\rm T}$, respectively, with the band index $\tilde{n}$ (including the spin degree of freedom) and the Bloch vector $\tilde{{\bf k}}$. Note that $|\tilde{n},\tilde{{\bf k}}\rangle$ is written as a linear combination of $|\tilde{\mathbf{k}},\tilde{X},\tilde{s}\rangle$ of the same $\tilde{{\bf k}}$. The summation over $\tilde{{\bf k}}$ in Eq. (5) is taken according to the condition Eq. (3).

The low-energy Hamiltonian is obtained by expanding $H_{\rm eff}({\bf k})$ around the valley center $\mathbf{K}_{\xi}\equiv-\xi(2\mathbf{a}^{*}_{1}+\mathbf{a}^{*}_{2})/3$, where $\xi=\pm 1$ is the valley index. Within the linear term, $H_{\rm G}$ is approximated by $H^{(\xi)}_{G}({\bf k})=-\hbar v(\mathbf{k}-\mathbf{K}_{\xi})\cdot(\xi\sigma_{x},\sigma_{y})$, where $v$ is the band velocity of graphene, and $\sigma_{x}$ and $\sigma_{y}$ are Pauli matrices for the sublattice space $X=p_{z}^{A},p_{z}^{B}$. Ando (2009). For the proximity SOC term, we only take the zero-th order $V_{\rm eff}({\bf K}_{\xi})\equiv V^{(\xi)}_{\rm eff}$. Now that the transfer integral $T_{\tilde{X}X}({\bf R})$ attenuates exponentially and so does its Fourier transform $t_{\tilde{X}X}(\mathbf{q})$, it suffices to keep only a few $\tilde{{\bf k}}$’s in the summation of Eq. (5). For ${\bf k}={\bf K}_{\xi}$, the dominant contribution comes from three points, $\tilde{\mathbf{k}}=\mathbf{K}_{\xi}+\xi\tilde{{\bf a}}^{*}_{1},\mathbf{K}_{\xi}+\xi({\bf a}^{*}_{1}+\tilde{{\bf a}}^{*}_{2}),\mathbf{K}_{\xi}+\xi({\bf a}^{*}_{1}+{\bf a}^{*}_{2}-\tilde{{\bf a}}^{*}_{1}-\tilde{{\bf a}}^{*}_{2})$, while the effect of other $\tilde{{\bf k}}$’s are negligibly small. In this way, the effective proximity potential can be obtained by considering TMDC Bloch states at only three wave points, and the corresponding computing cost is considerably low.

We can show that $V^{(\xi)}_{\rm eff}$ can be written as,

[TABLE]

where $s_{i}\,(i=x,y,z)$ is the Pauli matrix for spin. It is explicitly written in a matrix form,

[TABLE]

where the bases are arranged by order of $(X,s)=(A,\uparrow)$, $(B,\uparrow)$, $(A,\downarrow)$ and $(B,\downarrow)$. The difference in the diagonal elements $\lambda$ leads to the spin splitting between spin up and spin down, and the off-diagonal term $\lambda_{R}$ mixes the different spins. The term with $\lambda_{R}$ is similar to the Rashba spin-orbit coupling Kane and Mele (2005); Wang et al. (2015a) but here the spin axis can be rotated by an angle $\phi$ on $xy$-plane. The energy gap at the charge neutral point is given by $E_{\rm gap}=|\lambda\lambda_{R}|/(\lambda^{2}+\lambda_{R}^{2})^{1/2}$. The spin splitting width in large $k$ is given by $E_{\rm split}=(\lambda^{2}+\lambda_{R}^{2})^{1/2}$. The effective Hamiltonian $H^{(\xi)}_{\rm eff}=H^{(\xi)}_{G}({\bf k})+V^{(\xi)}_{\rm eff}$ is formally equivalent with that of the asymmetric bilayer graphene, McCann (2006); McCann and Koshino (2013) where the spin up and down correspond to layer 1 and 2, respectively, and $\lambda_{R}$ and $\lambda$ to the interlayer coupling and the interlayer asymmetric potential, respectively.

The form of Eq. (6) is forced by the symmetry of the system. The terms in Eq. (6) are generally allowed in the time reversal symmetry $\cal{T}$ and the $C_{3}$ symmetry. Actually, the term proprotional to $\sigma_{z}$ (different on-site energies at $A$ and $B$ sites) is also possible under $\cal{T}$ and $C_{3}$ Wang et al. (2015a), while it is prohibited by the incommesurability between graphene and TMDC, as explained in Supplementary Material. An additional space symmetry imposes a constraint on $V^{(\xi)}_{\rm eff}$. At $\theta=0$, the reflection symmetry $R_{x}:(x,y,z)\to(-x,y,z)$ requires $e^{i\phi}$ is real. At $\theta=30^{\circ}$, the reflection symmetry $R_{y}:(x,y,z)\to(x,-y,z)$ requires real $e^{i\phi}$ and also $\lambda=0$, i.e., the SOC is dominated by the Rashba term. The detailed argument of the symmetry consideration is presented in Supplementary Material.

We numerically calculate $V^{(\xi)}_{\rm eff}$ for MoS2, WS2, MoSe2 and WSe2. Figure 2 summarizes the results, where $\lambda$, $\lambda_{R}$, the central energy gap $E_{\rm gap}$ and the spin splitting $E_{\rm split}$ are plotted against the twist angle $\theta$. In Fig. 3, we present the band structures for each system at the rotation angles $\theta=0^{\circ},15^{\circ}$, and $30^{\circ}$. In the band plots of MoS2 and WS2 at $\theta=0^{\circ}$, the dotted green line indicates the first DFT calculations, from which we extract the interlayer hopping parameters. For the DFT calculation, we assume the approximate commensurate lattice structure of which unit cell is comprised of 3$\times$3 supercell of MoS2 and 4$\times$4 of graphene, and use Quantum EspressoGiannozzi et al. (2009, 2017) with the generalized gradient approximationPerdew et al. (1996). We can see that the effective model well reproduces the DFT band structure.

For the angle dependence, we find that $\lambda$ and $\lambda_{R}$ are greatly enhanced by rotation, and they take the maximum around $\theta\sim 20^{\circ}$. For WS2, in particular, the maximum splitting is about 5 times as large as that of 0∘. At 30*∘*, the parameter $\lambda$ vanishes and the $V_{\rm eff}$ is dominated by $\lambda_{R}$ as expected. There the band structure is formally equivalent to the symmetric AB-stacked bilayer graphene, and the expectation value of spin lies on the $xy$-plane.

The enhancement of the spin-splitting near $20^{\circ}$ can be understood by considering the second-order process, Eq. (5). The amplitude of $V_{\rm eff}$ is related to the spin splitting of the TMDC bands at $\tilde{\mathbf{k}}$ points which are hybridized with the graphene’s Dirac point. Figure 4(a) illustrates the positions of the three dominant $\tilde{\mathbf{k}}$’s for $\xi=+$, in WS2 with $\theta=0^{\circ},17.9^{\circ}$, and $30^{\circ}$. Figure 4(b) presents the band structure of WS2 with the vertical dashed lines indicating the $\tilde{\mathbf{k}}$’s for the three rotation angles. Now the lowest valence band of TMDC makes the greatest contribution to $V_{\rm eff}$, as it is the closest to the graphene’s Dirac point energy (black horizontal line), leading to a small denominator in Eq. (5). We can see that the $\tilde{\mathbf{k}}$ point for $17.9^{\circ}$ happens to be very close to the $Q$-valley, where the magnitude of the spin splitting much greater than in other angles. This qualitatively explains the sharp rise of $\lambda$ and $\lambda_{R}$ around $20^{\circ}$. Actually, the spin splitting can be even enhanced by shifting the relative energy between graphene and TMDC. Figure 4(c) plots the angle dependence of the spin-splitting of WS2 with $E_{T}-E_{G}=-0.08$ eV, where the graphene’s Fermi energy (blue horizontal line) hits the bottom of the Q-valley. Although the Fermi energy is just a little higher than in Fig. 2(b), the maximum spin splitting sharply increases to 20 meV, about 10 times as big as in $\theta=0$, because the denominator in Eq. (5) becomes very small. This suggests that tuning of the spin-orbit coupling would be possible using the external gate voltage.

Finally, the graphene under the proximity potential has the non-zero valley Hall conductivity when the Fermi energy lies in the central gap. The Hall conductivity of each valley sector can be calculated using the Berry curvature as

[TABLE]

where ${\bf a}_{n}({\bf k})=-i\langle u_{n{\bf k}}|\nabla_{\bf k}|u_{n{\bf k}}\rangle$ is the Berry connection, $u_{n{\bf k}}$ is the Bloch function (eigenvector of $H^{(\xi)}_{\rm eff}$) of the band $n$, and occ. stands for the occupied valence bands ($n=1,2$). As a result, the valley Hall conductivity becomes $\sigma^{(+)}_{xy}-\sigma^{(-)}_{xy}=-(2e^{2}/h){\rm sgn}(\lambda)$.

To conclude, we have studied the proximity spin-orbit interaction in graphene-TMDC bilayers stacked with arbitrary twist angles. By using the perturbational approach based on the tight-binding model, we derived the effective Hamiltonian of graphene as a continuous function of the twist angle $\theta$, and found that the magnitude of SOC is greatly enhanced by the rotation. We also show that the SOC sharply rises when the graphene’s Dirac point is shifted toward TMDC band, by applying the gate voltage. The theoretical method proposed here does not need the exact lattice matching, so that it is applicable to any incommensurate bilayer systems which cannot be treated by the DFT calculation. It would be useful for the twist-angle engineering of a wide variety of van der Waals proximity effects, including ferromagnetism and superconductivty. M. K. acknowledges the financial support of JSPS KAKENHI Grant No. JP17K05496.