Momentum dependence of spin–orbit interaction effects in single-layer and multi-layer transition metal dichalcogenides

The TMDs MX₂ are composed, in their bulk configuration, of two-dimensional X–M–X layers stacked on top of each other, coupled by weak van der Waals forces. The M atoms are ordered in a triangular lattice, each of them bonded to six X atoms located in the top and bottom layers, forming a sandwiched material. Our starting point will be a 11-band TB spinless model which, for the single-layer, considers the five d orbitals of the metal atom M and the three p orbitals for each of the two chalcogen atoms X in the top and bottom layer [3]. We can introduce a Hilbert base defined by the 11-fold vector:

$$\begin{eqnarray}&&\begin{array}{l} \phi _{i}^{\dagger }=(p_{i,x,t}^{\dagger },p_{i,y,t}^{\dagger },p_{i,z,t}^{\dagger },d_{i,3{{z}^{2}}-{{r}^{2}}}^{\dagger },d_{i,{{x}^{2}}-{{y}^{2}}}^{\dagger }, \\ \quad d_{i,xy}^{\dagger },d_{i,xz}^{\dagger },d_{i,yz}^{\dagger },p_{i,x,b}^{\dagger },p_{i,y,b}^{\dagger },p_{i,z,b}^{\dagger }), \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where $d_{i,\alpha }^{\dagger }$ creates an electron in the orbital α of the M atom in the i-unit cell, $p_{i,\alpha ,t}^{\dagger }$ creates an electron in the orbital α of the top (t) layer atom X in the i-unit cell, and $p_{i,\alpha ,b}^{\dagger }$ creates an electron in the orbital α of the bottom (b) layer atom X in the i-unit cell. After an appropriate unitary transformation, the spinless (sl) representation of the single-layer (1L) Hamiltonian can be expressed in the block form

$$\begin{eqnarray}\hat{H}_{1{\rm L}}^{{\rm sl}}(k)=\left( \begin{array}{ccccccccccccccc} {{{\hat{H}}}_{{\rm E}}} & {{{\hat{0}}}_{6\times 5}} \\ {{{\hat{0}}}_{5\times 6}} & {{{\hat{H}}}_{{\rm O}}} \\ \end{array} \right),\end{eqnarray} \tag{ 2 }$$

where ${{\hat{H}}_{{\rm E}}}$ and ${{\hat{H}}_{{\rm O}}}$ are a 6 × 6 and 5 × 5 blocks with even (E) and odd (O) parity respectively upon the mirror inversion $z\to -z$, and ${{\hat{0}}_{m\times n}}$ denotes $m\times n$ zero matrices [3]. In particular, ${{\hat{H}}_{{\rm E}}}$ is built from hybridizations of the d_xy, ${{d}_{{{x}^{2}}-{{y}^{2}}}}$, ${{d}_{3{{z}^{2}}-{{r}^{2}}}}$ orbitals of the metal M with the symmetric (antisymmetric) combinations of the p_x, p_y (p_z) orbitals of the top and bottom chalcogen atoms X. On the other hand, the odd block, ${{\hat{H}}_{{\rm O}}}$, is made by hybridizations of the d_xz and d_yz orbitals of M with the antisymmetric (symmetric) combinations of the p_x, p_y (p_z) orbitals of the X atom in the top and bottom layers. Explicit expressions for all the matrix elements in terms of the Slater–Koster parameters were obtained in [3], and we notice that the 6 × 6 even block ${{\hat{H}}_{{\rm E}}}$ contains the relevant orbital contribution for the states of the upper valence band and the lower conduction band.

In the context of the present TB model, we include the SOC term in the Hamiltonian by means of a pure atomic SO interaction acting on both the metal and chalcogen atoms. Explicitly we consider here the SOC given by:

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm SO}}}=\mathop{\sum }\limits_{a}\frac{{{\lambda }_{a}}}{\hbar }{{\hat{L}}_{a}}\cdot {{\hat{S}}_{a}},\end{eqnarray} \tag{ 3 }$$

where ${{\lambda }_{a}}$, the intra-atomic SOC constant, depends on the specific atom ($a=M,X$). ${{\hat{L}}_{a}}$ is the atomic orbital angular momentum operator and ${{\hat{S}}_{a}}$ is the electronic spin operator [33–35]. It is convenient to use the representation

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm SO}}}=\mathop{\sum }\limits_{a}\frac{{{\lambda }_{a}}}{\hbar }\left( \frac{\hat{L}_{a}^{+}\hat{S}_{a}^{-}+\hat{L}_{a}^{-}\hat{S}_{a}^{+}}{2}+\hat{L}_{a}^{z}\hat{S}_{a}^{z} \right),\end{eqnarray} \tag{ 4 }$$

where (omitting now for simplicity the atomic index a):

$$\begin{eqnarray}{{\hat{S}}^{+}}=\left( \begin{array}{ccccccccccccccc} 0 & 1 \\ 0 & 0 \\ \end{array} \right),\quad {{\hat{S}}^{-}}=\left( \begin{array}{ccccccccccccccc} 0 & 0 \\ 1 & 0 \\ \end{array} \right),\quad {{\hat{S}}^{z}}=\frac{1}{2}\left( \begin{array}{ccccccccccccccc} 1 & 0 \\ 0 & -1 \\ \end{array} \right).\end{eqnarray} \tag{ 5 }$$

In a similar way, the orbital angular momentum operator $\hat{L}$ acts on the states $|l,m\rangle$ as

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{L}}}^{\pm }}|l,m\rangle & = & \hbar \sqrt{l(l+1)-m(m\pm 1)}\ |l,m\pm 1\rangle , \\ {{{\hat{L}}}^{z}}|l,m\rangle & = & \hbar m\ |l,m\rangle , \\ \end{array}\end{eqnarray} \tag{ 6 }$$

where l refers to the orbital momentum quantum number and m to its z component.

We choose the orbital basis set in the following manner:

$$\begin{eqnarray}\begin{array}{rcl} |{{p}_{z}}\rangle & = & |1,0\rangle \\ |{{p}_{x}}\rangle & = & -\frac{1}{\sqrt{2}}[|1,1\rangle -|1,-1\rangle ] \\ |{{p}_{y}}\rangle & = & \frac{i}{\sqrt{2}}[|1,1\rangle +|1,-1\rangle ] \\ |{{d}_{3{{z}^{2}}-{{r}^{2}}}}\rangle & = & |2,0\rangle \\ |{{d}_{xz}}\rangle & = & -\frac{1}{\sqrt{2}}[|2,1\rangle -|2,-1\rangle ] \\ |{{d}_{yz}}\rangle & = & \frac{i}{\sqrt{2}}[|2,1\rangle +|2,-1\rangle ] \\ |{{d}_{{{x}^{2}}-{{y}^{2}}}}\rangle & = & \frac{1}{\sqrt{2}}[|2,2\rangle +|2,-2\rangle ] \\ |{{d}_{xy}}\rangle & = & -\frac{i}{\sqrt{2}}[|2,2\rangle -|2,-2\rangle ]. \\ \end{array}\end{eqnarray} \tag{ 7 }$$

We further simplify the problem by introducing the aforementioned symmetric (S) and antisymmetric (A) combination of the p orbitals of the top (t) and bottom (b) X layers:

$$\begin{eqnarray}\begin{array}{rcl} |{{p}_{\alpha ,{\rm S}}}\rangle & = & \frac{1}{\sqrt{2}}[|{{p}_{\alpha ,t}}\rangle +|{{p}_{\alpha ,b}}\rangle ], \\ |{{p}_{\alpha ,{\rm A}}}\rangle & = & \frac{1}{\sqrt{2}}[|{{p}_{\alpha ,t}}\rangle -|{{p}_{\alpha ,b}}\rangle ]. \\ \end{array}\end{eqnarray} \tag{ 8 }$$

The total Hamiltonian, including the SO interaction for the single-layer, can be now written as

$$\begin{eqnarray}&&{{\hat{H}}_{1{\rm L}}}(k)=\hat{H}_{1{\rm L}}^{{\rm sl}}(k)\otimes {{\mathbb{1}}_{2}}+\hat{H}_{1{\rm L}}^{{\rm SO}},\end{eqnarray} \tag{ 9 }$$

where the SOC term $\hat{H}_{1{\rm L}}^{{\rm SO}}$ is

$$\begin{eqnarray}\hat{H}_{1{\rm L}}^{{\rm SO}}=\left( \begin{array}{ccccccccccccccc} {{{\hat{M}}}^{\uparrow \uparrow }} & {{{\hat{M}}}^{\uparrow \downarrow }} \\ {{{\hat{M}}}^{\downarrow \uparrow }} & {{{\hat{M}}}^{\downarrow \downarrow }} \\ \end{array} \right),\end{eqnarray} \tag{ 10 }$$

and where

$$\begin{eqnarray}{{\hat{M}}^{\sigma \sigma }}=\left( \begin{array}{ccccccccccccccc} \hat{M}_{{\rm EE}}^{\sigma \sigma } & {{{\hat{0}}}_{6\times 5}} \\ {{{\hat{0}}}_{5\times 6}} & \hat{M}_{{\rm OO}}^{\sigma \sigma } \\ \end{array} \right)\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}{{\hat{M}}^{\sigma \bar{\sigma }}}=\left( \begin{array}{ccccccccccccccc} {{{\hat{0}}}_{6\times 6}} & \hat{M}_{{\rm EO}}^{\sigma \bar{\sigma }} \\ \hat{M}_{{\rm OE}}^{\sigma \bar{\sigma }} & {{{\hat{0}}}_{5\times 5}} \\ \end{array} \right).\end{eqnarray} \tag{ 12 }$$

Here we have chosen the spin notation $\bar{\sigma }=\downarrow$ ($\bar{\sigma }=\uparrow$) when $\sigma =\uparrow$ ($\sigma =\downarrow$).

The different blocks $\hat{M}_{{\rm EE}}^{\sigma \sigma }$, $\hat{M}_{{\rm OO}}^{\sigma \sigma }$, $\hat{M}_{{\rm EO}}^{\sigma \bar{\sigma }}$, $\hat{M}_{{\rm OE}}^{\sigma \bar{\sigma }}$, that constitute the above 22 × 22 matrix, are explicitly reported in the appendix. We notice here that, in the most general case, the SO interaction couples the E and O sectors of the 22 × 22 TB matrix. Such mixing arises in particular from the spin-flip/spin–orbital processes associated with the transverse quantum fluctuation described by the first two terms of equation (4). The effective relevance of these terms can now be directly investigated in a simple way, pointing out the advantages of a TB model with respect to first-principles calculations. The explicit analysis of this issue is discussed in section 3. We anticipate here that the effects of the off-diagonal spin-flip terms result to be negligible for all the cases of interest here. This is essentially due to the fact that such processes involve virtual transitions towards high-order energy states [17]. At a very high degree of accuracy, we are thus justified in neglecting the spin-flip terms and retaining in (4) only the spin-conserving terms $\propto \;{{\lambda }_{a}}\hat{L}_{a}^{z}\hat{S}_{a}^{z}$. An immediate consequence of that is that the even and odd sectors of the Hamiltonian remain uncoupled, allowing us to restrict our analysis, for the low-energy states of the valence and conduction bands, only to the E sector.

Once introduced the TB model for a single-layer in the presence of SOC, it is quite straightforward to construct a corresponding theory for the bulk and bilayer systems by including the relevant inter-layer hopping terms in the Hamiltonian. Considering that the unit cell is now doubled, we can thus write the Hamiltonian for bulk MX₂ in the presence of SOC in the matrix form:

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm Bulk}}}(k)=\hat{H}_{{\rm Bulk}}^{{\rm sl}}(k)\otimes {{\mathbb{1}}_{2}}+\hat{H}_{{\rm Bulk}}^{{\rm SO}},\end{eqnarray} \tag{ 13 }$$

which is a 44 × 44 matrix due to the doubling of the unit cell with respect to the single-layer case discussed in section 2.1.

Here $\hat{H}_{{\rm Bulk}}^{{\rm sl}}(k)$ represents the spinless Hamiltonian for the bulk system,

$$\begin{eqnarray}\hat{H}_{{\rm Bulk}}^{{\rm sl}}(k)=\left( \begin{array}{ccccccccccccccc} \hat{H}_{1}^{{\rm sl}} & {{{\hat{H}}}_{\bot ,{\rm Bulk}}} \\ \hat{H}_{\bot ,{\rm Bulk}}^{\dagger } & \hat{H}_{2}^{{\rm sl}} \\ \end{array} \right),\end{eqnarray} \tag{ 14 }$$

where $\hat{H}_{i}^{{\rm sl}}$ describes the spinless Hamiltonian (i.e. in the absence of SOC) for the layer i = 1,2, while ${{\hat{H}}_{\bot ,{\rm Bulk}}}$ accounts for the 11 × 11 Hamiltonian describing interlayer hopping between X atoms beloging to different layers. We remind that $\hat{H}_{2}^{{\rm sl}}$ is related to $\hat{H}_{1}^{{\rm sl}}$ through the following relation dictated by the lattice structure [3]:

$$\begin{eqnarray}&&H_{2,\alpha ,\beta }^{{\rm sl}}({{k}_{x}},{{k}_{y}})={{P}_{\alpha }}{{P}_{\beta }}H_{1,\alpha ,\beta }^{{\rm sl}}({{k}_{x}},-{{k}_{y}}),\end{eqnarray} \tag{ 15 }$$

where ${{P}_{\alpha }}\;=\;+(-)1$ if the orbital α has even (odd) symmetry with respect to $y\to -y$. Furthermore, the (spin-diagonal) interlayer term ${{\hat{H}}_{\bot ,{\rm Bulk}}}$ can be written as:

$$\begin{eqnarray}{{\hat{H}}_{\bot ,{\rm Bulk}}}(k)=\left( \begin{array}{ccccccccccccccc} {{{\hat{I}}}_{{\rm E}}}{\rm cos} \zeta & {{{\hat{I}}}_{{\rm EO}}}{\rm sin} \zeta \\ -\hat{I}_{{\rm EO}}^{{\rm T}}{\rm sin} \zeta & {{{\hat{I}}}_{{\rm O}}}{\rm cos} \zeta \\ \end{array} \right),\end{eqnarray} \tag{ 16 }$$

where $\zeta ={{k}_{z}}c/2$ (c being the vertical size of the unit cell in the bulk system), and where the matrices ${{\hat{I}}_{{\rm E}}}$, ${{\hat{I}}_{{\rm O}}}$ and ${{\hat{I}}_{\;{\rm EO}}}$ describe the inter-layer hopping between the p orbitals of the adjacent chalcogen atoms. One can notice that interlayer hopping leads, for an arbitrary wave-vector $k$, to a mixture of the E and O sectors of the Hamiltonian, which is accounted for by the term ${{\hat{I}}_{{\rm EO}}}$ in (16) [3]. The analysis is however simplified at specific high-symmetry points of the BZ, as we discuss below⁵ .

Finally, $\hat{H}_{{\rm Bulk}}^{{\rm SO}}$ in equation (13) accounts for the spin–orbit coupling (SOC) in the bulk system, and it can be written as:

$$\begin{eqnarray}\hat{H}_{{\rm Bulk}}^{{\rm SO}}=\left( \begin{array}{ccccccccccccccc} {{{\hat{M}}}^{\uparrow \uparrow }} & 0 & {{{\hat{M}}}^{\uparrow \downarrow }} & 0 \\ 0 & {{{\hat{M}}}^{\uparrow \uparrow }} & 0 & {{{\hat{M}}}^{\uparrow \downarrow }} \\ {{{\hat{M}}}^{\downarrow \uparrow }} & 0 & {{{\hat{M}}}^{\downarrow \downarrow }} & 0 \\ 0 & {{{\hat{M}}}^{\downarrow \uparrow }} & 0 & {{{\hat{M}}}^{\downarrow \downarrow }} \\ \end{array} \right),\end{eqnarray} \tag{ 17 }$$

where both the spin-diagonal (${{\hat{M}}^{\sigma \sigma }}$) and spin-flip (${{\hat{M}}^{\sigma \bar{\sigma }}}$) processes induced by the atomic SO interaction are present.

Equations (13)–(17) provide the general basic framework for a deeper analysis in more specific cases. In particular, as already mentioned above, the spin-flip terms triggered by SOC can be substantially neglected for all the cases of interest without loosing accuracy. The total Hamiltonian (13) can thus be divided in two 22 × 22 blocks $\hat{H}_{{\rm Bulk}}^{\sigma \sigma }(k)$ related by the symmetry $\hat{H}_{{\rm Bulk}}^{\uparrow \uparrow }(k)=\hat{H}_{{\rm Bulk}}^{\downarrow \downarrow }(-k)$. Further simplifications are available at specific symmetry points of the BZ. More specifically, we can notice that for k_z = 0 the E and O sectors remain uncoupled. Focusing, at low-energies for the conduction and valence bands, only on the E sector, we can write

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm Bulk},{\rm E}}}(k,{{k}_{z}}=0)=\hat{H}_{{\rm Bulk},{\rm E}}^{{\rm sl}}(k)+\hat{H}_{{\rm Bulk},{\rm E}}^{{\rm SO}},\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}\hat{H}_{{\rm Bulk},{\rm E}}^{{\rm sl}}(k)=\left( \begin{array}{ccccccccccccccc} {{{\hat{H}}}_{{\rm E},1}} & {{{\hat{I}}}_{{\rm E}}} & 0 & 0 \\ \hat{I}_{{\rm E}}^{\dagger } & {{{\hat{H}}}_{{\rm E},2}} & 0 & 0 \\ 0 & 0 & {{{\hat{H}}}_{{\rm E},1}} & {{{\hat{I}}}_{{\rm E}}} \\ 0 & 0 & \hat{I}_{{\rm E}}^{\dagger } & {{{\hat{H}}}_{{\rm E},2}} \\ \end{array} \right)\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}\hat{H}_{{\rm Bulk},{\rm E}}^{{\rm SO}}=\left( \begin{array}{ccccccccccccccc} \hat{M}_{{\rm EE}}^{\uparrow \uparrow } & 0 & 0 & 0 \\ 0 & \hat{M}_{{\rm EE}}^{\uparrow \uparrow } & 0 & 0 \\ 0 & 0 & \hat{M}_{{\rm EE}}^{\downarrow \downarrow } & 0 \\ 0 & 0 & 0 & \hat{M}_{{\rm EE}}^{\downarrow \downarrow } \\ \end{array} \right),\end{eqnarray} \tag{ 20 }$$

where the explicit expression of each block Hamiltonian is also reported in the appendix.

The Hamiltonian for the bilayer can also be derived in a very similar form as in the bulk case. In particular, we can write:

$$\begin{eqnarray}&&{{\hat{H}}_{2{\rm L}}}(k)=\hat{H}_{2{\rm L}}^{{\rm sl}}(k)+\hat{H}_{2{\rm L}}^{{\rm SO}}.\end{eqnarray} \tag{ 21 }$$

Since we are considering intrinsic SOC, thus it is not affected by the interlayer coupling. Therefore, we have $\hat{H}_{2{\rm L}}^{{\rm SO}}=\hat{H}_{{\rm Bulk}}^{{\rm SO}}$, where $\hat{H}_{{\rm Bulk}}^{{\rm SO}}$ is defined in equation (20).

On the other hand, similar to the bulk case in equation (14), the spinless TB term $\hat{H}_{2{\rm L}}^{{\rm sl}}(k)$ for the bilayer case can be written as:

$$\begin{eqnarray}\hat{H}_{2{\rm L}}^{{\rm sl}}(k)=\left( \begin{array}{ccccccccccccccc} \hat{H}_{1}^{{\rm sl}} & {{{\hat{H}}}_{\bot ,2{\rm L}}} \\ \hat{H}_{\bot ,2{\rm L}}^{\dagger } & \hat{H}_{2}^{{\rm sl}} \\ \end{array} \right),\end{eqnarray} \tag{ 22 }$$

where now

$$\begin{eqnarray}{{\hat{H}}_{\bot ,2L}}(k)=\frac{1}{2}\left( \begin{array}{ccccccccccccccc} {{{\hat{I}}}_{{\rm E}}} & {{{\hat{I}}}_{{\rm EO}}} \\ -\hat{I}_{{\rm EO}}^{{\rm T}} & {{{\hat{I}}}_{{\rm O}}} \\ \end{array} \right).\end{eqnarray} \tag{ 23 }$$

Note that equation (23) can be obtained as limiting case of equation (16) by setting $\zeta =\pi /4$, corresponding to the effective uncoupling of bilayer blocks.

The role of the SOC on the spin–orbital-valley entanglement at the band edge at K of the single-layer and bilayer compounds has been previously discussed in the literature, using mainly low-energy effective Hamiltonians focused on the role of the transition metal M d-orbitals and of their corresponding SOC [5, 8, 28–32]. Such scenario can be now well reproduced by the present TB model and generalized to the whole BZ.

The SOC, in particular, is expected to be most relevant for the band edges of the valence band at the K point, whose orbital content is mainly associated with the d_xy and ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ orbitals of the transition metal. A large band splitting induced by the SOC is thus predicted in this case. Such feature is indeed well captured by the TB model. In figures 3(a) and (b) we show the Fermi surfaces obtained with the present TB model, including atomic SOC, for a finite hole-doping probing the valence band of both single-layer and bulk compounds. In order to point out the different physics occurring close to the different band edges at K and Γ points, we show here Fermi surfaces corresponding to a sizable negative Fermi energy cutting both edges at K and Γ. In particular, the central Fermi pocket located around Γ appears to be spin degenerate, for both single-layer and bulk systems since its orbital character is mainly due to the ${{d}_{3{{z}^{2}}-{{r}^{2}}}}$ orbitals of M and to the p_z orbitals of X [3], both of them with L_z = 0. On the other hand, the pockets around K and K' are mainly due to the ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ and d_xy orbitals of the metal M (with $|m|=2$), plus a minor component of p_x and p_y orbitals of the chalcogen X (with $|m|=1$). This results in a finite SOC splitting of the valence band at the K and K' points, due mainly to first order SOC on the d orbitals of M. Furthermore, because of the lack of inversion symmetry in single layer samples (or in multi-layer samples with an odd number of layers), the spin degeneracy is lifted, presenting an opposite spin polarization on different valleys [5]. This feature is well reproduced by our model and shown in figure 3(a), where Fermi surfaces with main ${{S}_{z}}=\uparrow$ character are denoted by solid blue lines, while Fermi surfaces with main ${{S}_{z}}=\downarrow$ character are denoted by dashed blue lines. On the other hand, the Fermi surfaces of hole-doped bulk MoS₂, for the same E_F, are shown in figure 3(b). Since the maximum of the valence band for the bulk compound, because of the interlayer coupling, is located at the Γ point (see the band structure of figure 1(c)), the central pocket in figure 3(b) is considerably larger than in figure 3(a) for single layer samples. The the double Fermi surfaces around the K and K' points in 3(b) are spin degenerate, as impose by inversion symmetry. A recent set of ARPES measurements for MoS₂ and MoSe₂ [44] have shown the importance of the SOC in the band structure, obtaining experimental constant energy contours in very much agreement with those presented in figures 3(a) and (b).

Zoom In Zoom Out Reset image size

Although smaller and less noticed [17–19, 22–24], a spin-valley coupling is present also for the conduction band edge of the single-layer systems at the K and K' points. It is important to remind here that the orbital character in these points of the BZ is mainly associated with the ${{d}_{3{{z}^{2}}-{{r}^{2}}}}$ orbital (with m = 0) of the transition metal M, but with a finite contribution from the p_x and p_y orbitals of the chalcogen, with $m=\pm 1$) [3]. The SOC of the chalcogen atom X, mainly through the diagonal term $L_{X}^{z}S_{X}^{z}$, results thus in a smaller but finite splitting of the conduction band edge, as it can be also inferred by the Fermi surfaces for electron-doped single-layer compounds, as shown in figure 3(c). It is worth to stress that, although the resulting spin-induced splitting can be quite small, the entanglement between band splitting, spin and valley degrees appears to be quite strong, so that the lower band is $\uparrow$ polarized and the upper band $\downarrow$ polarized (or viceversa, depending on the valley). Note also that, although the atomic SOC due to the sulfur in MoS₂ or WS₂ is not very large, it can be of importance for Se compounds (with a larger atomic mass than sulfur), as MoSe₂ or WSe₂. The role of the SOC on the chalcogen atom will be analyzed in more detail in section 4.2.

Finally, we can note that, as previously discussed in [4] using first principles calculations, the SOC induces a finite band splitting in single-layer systems also at the Q point, with a corresponding spin-polarization. Also this feature is nicely captured by our TB model in the presence of atomic SOC on both chalcogen and transition metal atoms, as shown in figure 3(c) where we plot the Fermi surfaces of an electron-doped system with a Fermi level cutting only the lower conduction band at Q. As we can see, the TB model is able not only to reproduce the band splitting, but also to point out a strong degree of entanglement in this point of the BZ, with Fermi pockets with a strong spin polarization, and with an alternating polarization of the entangled spin/valley/orbital degrees of freedom along the six inequivalent valleys [45]. On the microscopic ground, we can notice that the main orbital character of the conduction bands at the Q point is due to a roughly equal distribution of the ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ and d_xy orbitals of the transition metal M, and of the p_x and p_y orbitals of the chalcogen atom X. Given the presence of a large contribution from both p- and d-orbitals, we expect these states to stem from a strong hybridization between X and M atoms, and hence to be highly sensitive to uniform and local strains and lattice distortions [46]. In addition, it should be kept in mind that the minimum of the conduction band at Q becomes the effective band edge in bilayer and multilayer compounds (as well as in strained single-layer systems). These considerations thus suggest the minima of the conduction band at the Q point as the most promising states for tuning the spin/orbital/valley entanglement in these materials by means of strain engineering [46] or (in multilayer systems) by means of electric fields [14].

Most of the existing theoretical works have focused on the effects of the SO interaction associated with the transition metal atom. Less attention has been paid, in general, to the SOC induced by the chalcogen atom. As we have seen in the previous section, however, the role of the SOC can be remarkably relevant also at the Q point of the BZ, resulting in a strong spin/orbital/valley entanglement also in this point, with the advantage to be extremely sensitive to the M-X hybridization and hence to the lattice effects. In addition, since the orbital content in this point is a mixture of d and p orbitals of the metal and the chalcogen atoms, the SOC is expected to be significantly driven not only by the d-orbital of the transition metal M, but also by the p orbitals of the chalcogen X atom. TB models can be quite useful to investigate this issue since we can easily tune the atomic SOC, keeping all the remaining Hamiltonian (Slater–Koster) parameters fixed, which permits to isolate the effects of the modified SOC without involving other structural and electronic changes. Figure 4 shows the effect of removing the SOC on either the Mo or the S atoms for the case of the single-layer of MoS₂. While the splitting of the valence band at the K point is fully caused by the SOC on the transition metal, the contributions to the splitting of the conduction band at Q from Mo and S are comparable.

Zoom In Zoom Out Reset image size

We can validate these findings by performing DFT calculations on the four compounds MX₂ with $M\;=$ Mo, W and $X\;=$ S, Se (all of them done using the experimental structure). In the DFT calculation, we can also turn on and off the SOC on a particular species, by removing the SO component of the pseudopotential [47]. Figure 5 shows the DFT results for the four compounds, including the SOC on all the atoms, and removing this coupling on either the chalcogen or the transition metal. In particular, the DFT results for MoS₂ shown in figure 5(a) agrees reasonably well with those of figure 4, signaling that the SOC splitting of the bands in the sulfur compounds is dominated by the contribution due to the transition metal atom.

Zoom In Zoom Out Reset image size

The importance of the SOC of the chalcogen atom is expected to be more remarkable for heavier atoms, such as selenium, instead of sulfur. In figures 5(c) and (d) we show the DFT band structure for MoSe₂ and WSe₂, isolating the contribution of the SOC due to the metal and to the chalcogen atoms. As expected, we observe that a relevant contribution to the SOC splitting of the bands is due to the Se atom. This can be seen by a noticeable splitting of the blue lines in figures 5(c) and (d) (for which the SOC due to the metal M has been switched off) which is governed by the SO interaction of the Se atoms. Interestingly, this effect is not relevant only at the Q point of the conduction band, but also at the K point of the valence band, for which the orbital weight of the p_x and p_y orbitals of Se is of only $\sim 20\%$ [3]. We conclude that, although for the MoS₂ and WS₂ the effect of the SOC of the chalcogenides does not have much effect on the band structure, when S is changed by Se, the effect is much noticeable.

Of special interest is the case of bilayer TMD, corresponding to a stack of two single layers in-plane rotated by ${{180}^{{}^\circ }}$ with respect to each other, such that the transition metal atoms of one layer are above the chalcogen atoms of the other layer. The two layers are bound by means of weak van der Waals interactions. The inter-layer hopping of electrons between different layers leads to a strong modification of the band structure, driving a transition from a direct gap semiconductor in single-layer systems to an indirect gap semiconductor in bilayer and multi-layer compounds. The inter-layer hopping links mainly the p orbitals of the chalcogen atoms X of different layers [3]. The result of this hopping is a splitting of the maximum of the valence band at the Γ point, which becomes the effective valence band edge, as well as a splitting of the minimum of the conduction band at the Q point which becomes the absolute minimum of the conduction band. This situation is shown in figure 6, where we report the band structure of bilayer MoS₂ and WS₂ calculated by DFT methods.

Zoom In Zoom Out Reset image size

A qualitative similar feature is observed also in other bilayer compounds, as MoSe₂ or WSe₂.

Contrary to single-layer MX₂, bilayer MX₂ presents point-center inversion symmetry [14, 16, 48]. Therefore, as we have discussed for the bulk case, the corresponding band structure remains spin degenerate even in the presence of SOC. However, since the SOC Hamiltonian does not couple orbitals of different layers, each single band preserves a finite entanglement between spin, valley and the layer index. Such spin-valley-layer coupling has been discussed in [48], where the authors focused on the relevance of this effect at the K point of the valence band. Here we notice that the same effect occurs also for the conduction band, and it can be thus relevant for electron-doped samples. Indeed for slightly electron-doped bilayer MoS₂ and WS₂ the Fermi surface presents six pockets centered at the inequivalent Q valleys of the BZ, and no pockets at the K and K' valleys. Interestingly, the SOC for the TMD families with stronger SO interaction, like WS₂ and WSe₂, can be larger than the inter-layer hopping, enhancing the spin/layer/valley entanglement. Then, although inversion symmetry forces each Fermi pocket to be spin degenerate, the layer polarization makes that each layer contributes with opposite spin in alternating valleys. This property can be of interest for valleytronics devices: by partially filling only one of the two subbands at the Q point of the conduction band, one would have a situation in which the upper layer contributes to three of the six valleys with spin-$\uparrow$, and with spin-$\downarrow$ to the other three valleys, whereas the opposite contribution is inferred from the bottom layer. This spin-valley coupling scenario resembles that of single-layer and bilayer MX₂ discussed in the literature, but for electron-doped samples, which is the kind of doping most commonly reported for those materials. Although we have focused in this section in the most simple multi-layer compound, which is the bilayer MX₂, the physics discussed above applies also to any multi-layer TMD with an even number of layers, because they contain the same symmetry properties as that of bilayer MX₂ discussed here.