Intervalley coupling by quantum dot confinement potentials in monolayer transition metal dichalcogenides

Semiconductor quantum dots (QDs) have been widely explored in the past several decades for technological applications such as lasers and medical markers [1]. Spin and other quantum degrees of freedom of a single electron or hole confined in QDs have also been extensively researched as potential carriers of information for quantum computing and quantum spintronics [2–6]. QDs in conventional semiconductors are realized either by forming nanocrystals or by lateral confinements in two-dimensional (2D) heterostructures. Such lateral confinements can be provided by the thickness variations of the heterostructures or patterned electrodes [2, 3]. The emergence of graphene [7, 8], a 2D crystal with single atom thickness, has offered a new system to host QDs. Because of the zero gap nature of graphene, QD confinement has to be facilitated by terminations of the monolayer, and the various geometries explored include the graphene nano-island connected to source and drain contacts via narrow graphene constrictions [9–11], and the graphene nanoribbon with armchair edges with patterned electrodes [12]. The properties of these graphene QDs are affected significantly by the edges, and precise control of edge terminations is thus desired.

Monolayer group-VIB transition metal dichalcogenides (TMDs) are newly emerged members of the 2D crystal family [13]. These TMDs have the chemical composition of MX₂ (M = Mo,W; X = S, Se), and the monolayer has a structure of a X–M–X covalently bonded hexagonal quasi-2D network. The 2D bulk has a direct bandgap in the visible frequency range [14, 15], ideal for optoelectronic applications. This bandgap also makes possible the realization of QDs without the assistance of termination edges [16]. QDs can be defined by lateral confinement potentials on an extended monolayer, e.g. by patterned electrodes, similar to the QDs in III–V heterostructures [2]. As the conduction and valence band edges are shifted in the same way by the electrodes, the QDs will have a unique type-II band alignment with its surroundings (see figure 1). Moreover, QD confinement potentials can also be realized by the lateral heterojunctions between different TMDs within a single uniform crystalline monolayer. Various band alignments can form between different TMDs with different bandgaps and workfunctions. Lateral heterostructures with MoSe₂ islands surrounded by WSe₂ on a crystalline monolayer have been realized very recently [17]. This can realize QD confinement, also with type-II band alignment, for electrons at the MoSe₂ region. The band edge discontinuity at the heterojunction, with an order of magnitude of 0.2–0.4 eV [18–22], forms the potential well (with a vertical wall).

Zoom In Zoom Out Reset image size

Besides the unique 2D geometries, monolayer TMD QDs are highly appealing because of the interesting properties of the 2D bulk. The conduction and valence band edges are both at the degenerate K and $-K$ valleys at the corners of the hexagonal Brillouin zone. Remarkably, the direct-gap optical transitions are associated with a valley dependent selection rule: left- (right-) handed circular polarized light excites interband transitions in the K ($-K$) valley only [23, 24]. Based on this selection rule, optical pumping of valley polarization [25–27], and optical generation of valley coherence [28], have been demonstrated. Another remarkable property of the 2D TMDs is the strong interplay between spin and the valley pseudospin [29–31]. These quantum degrees of freedom with versatile controllability suggest well that a single electron in monolayer TMD QDs can be a promising carrier of information for quantum spintronics, provided that the bulk properties of interest can be inherited by the QDs. Since these interesting bulk properties are associated with the valley pseudospin, a central problem is whether the valley is still a good quantum number and whether the bulk valley physics is preserved in the QDs. A natural concern is that the lateral QD confinement may result in valley hybridization, like in silicon QDs [32–39], which may completely change the valley physics in QDs.

In this paper, we present a systematic study on the intervalley coupling strength by QD confinement potentials on an extended TMD monolayer. The numerical calculations used two methods: (i) the envelope function method (EFM) [40–43], in conjunction with first-principles wavefunctions of the band-edge Bloch states of the 2D bulk; and (ii) the real-space tight-binding (RSTB) approach based on a three-band model for monolayer TMDs [44]. The EFM is limited to circular shaped QDs and is used primarily as a benchmark for the RSTB approach which can handle arbitrarily shaped QDs. The intervalley coupling is defined here as the off-diagonal matrix element, due to the QD potentials, between states with the same spin but opposite valley index. For circular QD potentials with various lateral size, potential depth and smoothness, the intervalley couplings from the two different approaches agree well, which justify both methods. We then use the RSTB method to address QD potentials of lower symmetry, including the triangular, hexagonal, and square shaped QDs, where the intervalley coupling strength is investigated as functions of lateral size, depth and smoothness of confinement potentials.

Our main findings are summarized below. For confinement potentials with the C₃ rotational symmetry (like that of the lattice), both the numerical results and the symmetry analysis show that the intervalley coupling strength depends sensitively on the central position of the potential: the coupling is maximized (zero) if the potential is centered at an M (X) site (see figure 1). The results stated below refer to the M-centered potentials. The intervalley coupling as a function of the lateral size of the QDs exhibits fast oscillations with an nearly exponentially decaying envelop. When the wall of the potential well changes from a vertical one to a sloping one, the intervalley coupling strength has a fast decrease by 2–3 orders of magnitude, and saturates when the length scale of the slope is beyond five lattice constants. For potentials with sloping walls, the intervalley coupling increases with increasing the potential depth. For potentials with vertical walls, the dependence on the potential depth is non-monotonic. Interestingly, when all parameters are comparable, the intervalley coupling can be smaller by several orders of magnitude in certain QD potentials. These include the circular shaped QDs, and the triangular and hexagonal shaped QDs with all sides along the zigzag crystalline axes. The latter two types of confinement potentials can be relevant in QDs formed by lateral heterostructures of different TMDs [17]. In contrast to graphene QDs [12] and silicon QDs [32–39], the intervalley coupling here is much smaller as the electron wavefunction vanishes at the boundary of the QDs. For all QDs studied, the largest intervalley coupling is upper bounded by 0.1 meV, found for small QDs with diameters of 20 nm and with confinement potentials of vertical walls. Such intervalley coupling is much smaller compared to the diagonal energy difference between states with the same spin but opposite valley index, which arises from the spin–valley coupling and ranges from several to several tens meV in different monolayer TMDs [44]. Therefore, our results mean that valley hybridization is in general negligible, and valley pseudospin in monolayer TMD QDs is a good quantum degree of freedom as in the 2D bulk, with the optical addressability for quantum controls.

The rest of the paper is organized as follows. In section 2, we formulate the EFM for the two-band ${\boldsymbol{k}} \cdot {\boldsymbol{p}}$ model of the monolayer TMDs [24], which can be exactly solved for circular shaped QDs. We then present the numerical results for intervalley coupling, calculated in conjunction with the first-principles wavefunctions. In section 3, we present the symmetry analysis on the M-centered and X-centered potentials with C₃ rotational symmetry. Section 4 is on the RSTB approach. Numerical results by RSTB method for circular QDs are first compared with the EFM results, and then the RSTB results for various shaped QDs is presented. The results are qualitatively the same and quantitatively comparable for different MX₂, and MoS₂ is presented as the example throughout the paper.

2.1. The model and method

We consider first a QD defined by a circular shaped potential U (r) in a TMD monolayer. Such a confinement potential can be characterized by three parameters: the lateral size, the depth, and the smoothness of the potential well. Without loss of generality, we assume the potential has the form

$$\begin{eqnarray}U(r)=\left\{ \begin{array}{ccccccccccccccc} 0, & r\lt R-\frac{L}{2} \\ \frac{V}{2}\left[ \frac{{\rm tanh} \frac{C(r-R)}{L/2}}{{\rm tanh} C}+1 \right], & R-\frac{L}{2}\leqslant r\leqslant R+\frac{L}{2} \\ V, & r\gt R+\frac{L}{2}, \\ \end{array} \right.\end{eqnarray} \tag{ 1 }$$

where R and V are respectively the lateral size and the depth of the potential well. L characterizes the smoothness of the potential well boundary (see figure 1(c)). Large L corresponds to smooth potentials (e.g. in confinement generated by patterned electrodes), while L = 0 corresponds to a potential well with a vertical wall (e.g. in lateral heterostructures of different TMDs). The parameter C = 2.5 is used here.

For typical confinement with R much larger than the lattice constants, the bound states in the QDs are formed predominantly from the band-edge Bloch states in the $\pm K$ valleys of the 2D bulk. In the EFM, we start off with the two-band ${\boldsymbol{k}} \cdot {\boldsymbol{p}}$ Hamiltonian derived for the band edges at the $\pm K$ valleys of monolayer TMDs [24, 44], where valley is a discrete index. Thus the EFM leads to confined wavefunctions formed with the K and $-K$ valley Bloch states, respectively, denoted as ${{\Psi }^{\tau ,s}}$ where $\tau =\pm$ is the valley index denoting the $\pm K$ valley, and $s\;=\;+(\uparrow )$ or $-(\downarrow )$ is the spin index. Taking these states as the basis, intervalley coupling here refers to the off-diagonal matrix elements between these states due to the confinement potential U (r) [32].

As confinement potential is spin-independent, intervalley coupling vanishes between states with opposite spin index. For holes, the band edges of the 2D bulk are spin–valley locked because of the giant spin–orbit coupling (SOC) [24, 45], i.e. valley K ($-K$) has spin up (down) only. Thus, intervalley coupling is absent for holes. We concentrate on the intervalley coupling of confined electron states ${{\Psi }^{+,s}}$ and ${{\Psi }^{-,s}}$ in the two valleys with the same spin index, which is defined as

$$\begin{eqnarray}&&V_{{\rm inter}}^{s}=\langle {{\Psi }^{+,s}}|U(r)|{{\Psi }^{-,s}}\rangle =\int {{\Psi }^{+,s}}{{\left( {\boldsymbol{r}} \right)}^{*}}U(r){{\Psi }^{-,s}}({\boldsymbol{r}} ){\rm d}{\boldsymbol{r}} .\end{eqnarray} \tag{ 2 }$$

$V_{{\rm inter}}^{s}$ here is the off-diagonal matrix element, due to the QD potentials, between states with the same spin but opposite valley index. This quantity determines to what degree the two bulk valleys can hybridize in the QD confinement, by competing with the diagonal energy difference between ${{\Psi }^{-,s}}$ and ${{\Psi }^{+,s}}$. If $V_{{\rm inter}}^{s}$ is much smaller than the diagonal energy difference, then the bulk valley index is still a good quantum number in QD. If the coupling is comparable or larger than the diagonal energy difference, then the QD eigenstates will be superpositions of the two bulk valleys. In the limit that $V_{{\rm inter}}^{s}$ is much larger than the diagonal energy difference, the QD eigenstates are the symmetric and antisymmetric superpositions of the two valleys, and their energy splitting is given by $2|V_{{\rm inter}}^{s}|$. Due to the time-reversal symmetry, $|V_{{\rm inter}}^{\uparrow }|=|V_{{\rm inter}}^{\downarrow }|$. Hence we consider here only spin-up case and drop the spin index s in superscripts for brevity.

To obtain ${{\Psi }^{\tau }}$ in EFM, we substitute $-{\rm i}\nabla$ for ${\boldsymbol{k}}$ in the ${\boldsymbol{k}} \cdot {\boldsymbol{p}}$ Hamiltonian of the 2D bulk[24, 44], and the confinement potential U (r) is added as an onsite energy term. The QD Hamiltonian in each valley is then [46–52]

$$\begin{eqnarray}{{H}^{\tau }}=\left[ \begin{array}{ccccccccccccccc} \frac{\Delta }{2}+U(r) & at(-{\rm i}\tau \frac{\partial }{\partial x}-\frac{\partial }{\partial y}) \\ at(-{\rm i}\tau \frac{\partial }{\partial x}+\frac{\partial }{\partial y}) & -\frac{\Delta }{2}+\tau s\lambda +U(r) \\ \end{array} \right],\end{eqnarray} \tag{ 3 }$$

where $\Delta \sim 2$ eV is the bulk bandgap of the TMD monolayer, a is the lattice constant, and t is the effective hopping. The bases of the Hamiltonian (3) are the conduction and valence Bloch states at $\pm K$ point, $\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )={{{\rm e}}^{{\rm i}\tau {\boldsymbol{K}} \cdot {\boldsymbol{r}} }}u_{\alpha }^{\tau }({\boldsymbol{r}} )$ ($\alpha ={\rm c}$, v and $u_{\alpha }^{\tau }$ is the cell periodic part), and the eigenfunction of (3) is a spinor of envelope function

$$\begin{eqnarray}&&{{\psi }^{\tau }}=\left[ \begin{array}{ccccccccccccccc} \psi _{{\rm c}}^{\tau }({\boldsymbol{r}} ) \\ \psi _{{\rm v}}^{\tau }({\boldsymbol{r}} ) \\ \end{array} \right],\end{eqnarray} \tag{ 4 }$$

which satisfy the Schrödinger equation ${{H}^{\tau }}{{\psi }^{\tau }}=E{{\psi }^{\tau }}$. The overall wavefunctions of the QD states are then

$$\begin{eqnarray}&&{{\Psi }^{\tau }}({\boldsymbol{r}} )=\psi _{{\rm c}}^{\tau }({\boldsymbol{r}} )\varphi _{{\rm c}}^{\tau }({\boldsymbol{r}} )+\psi _{{\rm v}}^{\tau }({\boldsymbol{r}} )\varphi _{{\rm v}}^{\tau }({\boldsymbol{r}} ).\end{eqnarray} \tag{ 5 }$$

Since U (r) has the rotational symmetry, it is easy to verify that the eigen spinor ${{\psi }^{\tau }}$ is of the form

$$\begin{eqnarray}&&{{\psi }^{\tau }}=\left[ \begin{array}{ccccccccccccccc} b_{{\rm c}}^{\tau }(r){{{\rm e}}^{{\rm i}m\theta }} \\ {\rm i}b_{{\rm v}}^{\tau }(r){{{\rm e}}^{{\rm i}(m+\tau )\theta }} \\ \end{array} \right],\end{eqnarray} \tag{ 6 }$$

where m is an integer, and $(r,\theta )$ are the polar coordinates of ${\boldsymbol{r}}$. $b_{{\rm c}}^{\tau }$ and $b_{{\rm v}}^{\tau }$ are real radial functions that satisfy

$$\begin{eqnarray}&&\left\{ \begin{array}{ccccccccccccccc} \tau \frac{\partial b_{{\rm c}}^{\tau }}{\partial r}-\frac{m}{r}b_{{\rm c}}^{\tau }=[-\frac{\Delta }{2}+s\tau \lambda +U(r)-E]b_{{\rm v}}^{\tau } \\ \tau \frac{\partial b_{{\rm v}}^{\tau }}{\partial r}+\frac{m+\tau }{r}b_{{\rm v}}^{\tau }=-[\frac{\Delta }{2}+U(r)-E]b_{{\rm c}}^{\tau } \\ \end{array} \right..\end{eqnarray} \tag{ 7 }$$

$b_{{\rm c}}^{\tau }$ and $b_{{\rm v}}^{\tau }$ are Bessel functions in the region where U (r) is a constant ($r\lt R-\frac{L}{2}$ or $r\gt R+\frac{L}{2}$), and need to be solved numerically in the region $R-\frac{L}{2}\leqslant r\leqslant R+\frac{L}{2}$. Boundary conditions are then matched at $r=R-\frac{L}{2}$ and $r=R+\frac{L}{2}$ to determine the solutions and eigenvalues. Bound states are found to exist in the energy range $(\frac{\Delta }{2},\frac{\Delta }{2}+V)$, as expected. Here t is used as the unit of energy and a as the unit of length.

The intervalley coupling defined in equation (2) can then be expressed as

$$\begin{eqnarray}&&{{V}_{{\rm inter}}}=\langle \mathop{\sum }\limits_{\alpha }\psi _{\alpha }^{+}{{{\rm e}}^{{\rm i}{\boldsymbol{K}} \cdot {\boldsymbol{r}} }}u_{\alpha }^{+}|U|\mathop{\sum }\limits_{\beta }\psi _{\beta }^{-}{{{\rm e}}^{-{\rm i}{\boldsymbol{K}} \cdot {\boldsymbol{r}} }}u_{\beta }^{-}\rangle =\mathop{\sum }\limits_{\alpha ,\beta }{{V}_{\alpha \beta }},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{{V}_{\alpha \beta }}=\mathop{\sum }\limits_{\Delta {\boldsymbol{G}} }{{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} ){{I}_{\alpha \beta }}(\Delta {\boldsymbol{G}} ),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{I}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )=\int \psi _{\alpha }^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\psi _{\beta }^{-}({\boldsymbol{r}} ){{{\rm e}}^{{\rm i}(\Delta {\boldsymbol{G}} -2{\boldsymbol{K}} )\cdot {\boldsymbol{r}} }}U(r){\rm d}{\boldsymbol{r}} ,\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&{{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )=\mathop{\sum }\limits_{{\boldsymbol{G}} }c_{\alpha }^{+}{{\left( {\boldsymbol{G}} \right)}^{*}}c_{\beta }^{-}({\boldsymbol{G}} +\Delta {\boldsymbol{G}} ).\end{eqnarray} \tag{ 11 }$$

In the above equations, $\alpha ,\beta \;\in \;${c, v}, and $c_{\alpha }^{\tau }({\boldsymbol{G}} )$ are the Fourier coefficients of the periodic part of the Bloch function $u_{\alpha }^{\tau }({\boldsymbol{r}} )={{\sum }_{{\boldsymbol{G}} }}c_{\alpha }^{\tau }({\boldsymbol{G}} ){{{\rm e}}^{{\rm i}{\boldsymbol{G}} \cdot {\boldsymbol{r}} }}$. ${\boldsymbol{G}}$ is a reciprocal lattice vector (RLV), and $\Delta {\boldsymbol{G}} ={\boldsymbol{G}} ^{\prime} -{\boldsymbol{G}}$ is also a RLV. Here we are interested in the QD ground states only, which correspond to the m = 0 case in equation (10) (see figure 2(a)). The coefficients $c_{\alpha }^{\tau }({\boldsymbol{G}} )$ are obtained here from the first-principles all-electron wavefunctions calculated by the ABINIT package [53]³ .

Zoom In Zoom Out Reset image size

We note that the Hamiltonian in equation (3) does not contain the conduction band spin–valley interaction [44, 54], which is of the Ising form that changes only the diagonal energies but not the wavefunctions of the states ${{\Psi }^{+,s}}$ and ${{\Psi }^{-,s}}$. Thus, this spin–valley interaction does not affect the intervalley coupling strength defined here, but it leads to an energy difference between states ${{\Psi }^{+,s}}$ and ${{\Psi }^{-,s}}$ that ranges from several to several tens of meV [44, 54]. Valley hybridization in the QDs is determined by the competition between $|{{V}_{{\rm inter}}}|$ and this diagonal energy difference.

2.2. Numerical results

In this subsection we present the numerical calculation results. The parameters from [24] are used: $t=1.1\;$ eV, $a=3.193\;$Å, $\Delta =1.66\;$ eV, and $2\lambda =0.15\;$ eV. We consider potentials centered at a Mo site, as shown in figure 1(a). As we will prove in section 3, the intervalley coupling is strongest for such Mo-centered confinement potential, and will vanish if the circular confinement potential is centered at a S site.

2.2.1. Energy levels and bound states

We first look at the energy eigenvalues and eigenstates of equation (3), which give the energy spectrum in the circular QD confinement, and the radial envelope functions to calculate the intervalley coupling (see equation (2)). The energy spectra are shown in figure 2(a) for R = 80, L = 10 and three well depths $V=$ 0.2, 0.05, and 0.005. For V = 0.2 and 0.05, only the low-energy parts of the spectra are shown. m is the azimuthal quantum number of the envelope function defined in equation (6). Obviously, the ground state has m = 0. The spectra shown are for the K valley ($\tau =1$) and it is obvious from equation (7) that states of τ, s, m are degenerate with states of $-\tau ,$ $-s$, $-m$. For the ground states, the radial envelope functions $b_{{\rm c}}^{+}(r)$ and $b_{{\rm v}}^{+}(r)$ are shown in figure 2(b), from which we can see that they are smooth and localized within the well and the $b_{{\rm v}}^{+}(r)$ components are orders of magnitude smaller compared to $b_{{\rm c}}^{+}(r)$. This means that for electron confined in the quantum dot, the wavefunction is predominantly contributed by the conduction band edge Bloch functions, and the contribution from the valence band edge is negligible. In general, the Hamiltonian (3) will be solved for addressing confinement of Dirac fermions, as we did here. Nevertheless, as shown here, the much simplified approach of effective mass approximation can be well justified for monolayer TMDs with the large bandgap. Importantly, the wavefunction has vanishing amplitude at the QD boundary (green shaded region in figure 2(b)). This is in contrast to the wavefunctions in graphene QDs [12] and silicon QDs [32–39] where significant valley hybridization is found.

2.2.2. The intervalley coupling

With the envelope functions given above, and the band-edge Bloch functions from the ABINIT package [53, see footnote 3], we can calculate the numerical values for the intervalley coupling (see equation (8)),

Figure 3 shows the intervalley coupling strength $|{{V}_{{\rm inter}}}|$ versus the radius of the potential R for circular QDs with fixed smoothness (L = 0 and 10) and depth (V = 0.005) of potential (see equation (1)). The curve with L = 10 is more than two orders of magnitude smaller than the one with L = 0, which can also be seen in figure 4 and is discussed in detail later. From figure 3(b) and (c), we can see that the dominant contribution to the intervalley coupling is from $|{{V}_{{\rm cc}}}|$, which is orders of magnitude larger than $|{{V}_{{\rm vv}}}|$, $|{{V}_{{\rm cv}}}|$, $|{{V}_{{\rm vc}}}|$ (see equation (9)), consistent with the fact that the wavefunction of the confined electron is predominantly from the conduction band edge Bloch functions. Figures 3(b) and (c) also clearly show that $|{{V}_{{\rm inter}}}|$ oscillates versus R in the scale of the lattice constant a. The oscillation with the size of QD is a generic feature of intervalley coupling, due to the large momentum space separation of the two valleys, which have been noted in previous works [37–39]. Apart from these fast oscillations, the intervalley coupling strength decreases with the increase of R in a nearly exponential way.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the $|{{V}_{{\rm inter}}}|$ as a function of L, the smoothness of the potential, for circular QDs with R = 80 and V = 0.005. The largest intervalley coupling is found at L = 0 which corresponds to a circular square well, where the potential has a sudden jump at the QD boundary. It can be seen that $|{{V}_{{\rm inter}}}|$ decreases monotonically as L increases up to a value of 5, and then it starts to oscillate rapidly around a constant value. The inset is the zoom in for $5\leqslant L\leqslant 10$, which shows that the oscillation is also in the length scale of lattice constant, similar to the $|{{V}_{{\rm inter}}}|$ versus R relation.

Figure 5 shows the dependence of $|{{V}_{{\rm inter}}}|$ on the potential depth V, when R = 80 and L = 0 and 10. We can see that $|{{V}_{{\rm inter}}}|$ increases with increasing V when L = 10, and it first increases and then decreases with increasing V when L = 0. In this figure, we also compare the intervalley couplings calculated with and without SOC. For the calculation without SOC, we set s = 0 in equation (3), and the plane-wave coefficients in equation (11) are also calculated without including the SOC in the ABINIT package. In figure 5, it is clear that intervalley couplings with and without SOC coincide when L = 10, and have small difference when L = 0. We also compared the $|{{V}_{{\rm inter}}}|$ versus R curve (fixing L = 10, V = 0.005) and $|{{V}_{{\rm inter}}}|$ versus L curve (fixing R = 80, V = 0.005), for the calculations with and without SOC, and the differences are also negligible. This means that SOC has negligible effects on $|{{V}_{{\rm inter}}}|$. This is expected as $|{{V}_{{\rm inter}}}|$ here is just the off-diagonal matrix element of the confinement potential, between the QD wavefunctions in the K and $-K$ valleys. The SOC in the K valleys is of the longitudinal form [24], which does not affect the wavefunction. The SOC, however, will cause spin splittings and hence affect the spectral of QD states in the K and $-K$ valleys (see figure 2(a)). When we consider the effect of intervalley coupling on the QD states, the off-diagonal matrix element $|{{V}_{{\rm inter}}}|$ will then compete with the differences in these diagonal energies, and the SOC induced spin splitting will play an important role.

Zoom In Zoom Out Reset image size

From the dependence of $|{{V}_{{\rm inter}}}|$ on R, L, and V discussed above, we find that the intervalley coupling strength in a circular confinement potential varies over several orders of magnitude, depending on the length scale, smoothness, and depth of the confinement. The magnitude of the intervalley coupling is small on the whole, in the order of μeV or below. This is much smaller compared to intervalley coupling in QDs formed in silicon inversion layers [32–39].

In the above section, we have seen that the intervalley coupling is very sensitive to the geometry of the confinement potential, as evident from the fast oscillations of $|{{V}_{{\rm inter}}}|$ versus R or L. Here, based on symmetry analysis, we will show that the intervalley coupling strength $|{{V}_{{\rm inter}}}|$ also depends on the position of the potential center, such that ${{V}_{{\rm inter}}}$ is largest when the potential is centered at a Mo atom (figure 1(a)) but vanishes if centered at a S atom (figure 1(b)). To explain this, we analyze ${{V}_{\alpha \beta }}$ (see equation (9)) by examining the behavior of ${{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )$ and ${{I}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )$ under the C₃ rotation.

Using the Fourier relation $c_{\alpha }^{\tau }({\boldsymbol{G}} )=\frac{1}{\Omega }{{\int }_{\Omega }}u_{\alpha }^{\tau }({\boldsymbol{r}} ){{{\rm e}}^{-{\rm i}{\boldsymbol{G}} \cdot {\boldsymbol{r}} }}\;{\rm d}{\boldsymbol{r}}$ in which Ω is the area of the 2D cell, ${{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )$ can be rewritten as

$$\begin{eqnarray}&&{{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )=\frac{1}{\Omega }{{\int }_{\Omega }}u_{\alpha }^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}u_{\beta }^{-}({\boldsymbol{r}} ){{{\rm e}}^{-{\rm i}\Delta {\boldsymbol{G}} \cdot {\boldsymbol{r}} }}\;{\rm d}{\boldsymbol{r}} \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\qquad \qquad =\;\frac{1}{\Omega }{{\int }_{\Omega }}\varphi _{\alpha }^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\varphi _{\beta }^{-}({\boldsymbol{r}} ){{{\rm e}}^{-{\rm i}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }}\;{\rm d}{\boldsymbol{r}} ,\end{eqnarray} \tag{ 13 }$$

where ${\boldsymbol{g}} \equiv \Delta {\boldsymbol{G}} -2{\boldsymbol{K}}$. Equation (12) shows that its integrand is a cell periodic function. This periodicity leads to the invariance of ${{S}_{\alpha \beta }}$ under the C₃ rotation of its integrand. For convenience, we rewrite ${{\mathbb{S}}_{\alpha \beta }}({\boldsymbol{g}} )\equiv {{S}_{\alpha \beta }}(\Delta {\boldsymbol{G}} )$. Using the form of integrand in equation (13), we get

$$\begin{eqnarray}\begin{array}{rcl} {{\mathbb{S}}_{\alpha \beta }}({\boldsymbol{g}} ) & = & \frac{1}{\Omega }{{\int }_{\Omega }}{{C}_{3}}\left[ \varphi _{\alpha }^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\varphi _{\beta }^{-}({\boldsymbol{r}} ){{{\rm e}}^{-{\rm i}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }} \right]\;{\rm d}{\boldsymbol{r}} \\ {} & = & \frac{1}{\Omega }{{\int }_{\Omega }}\left[ {{C}_{3}}{{\left( \varphi _{\alpha }^{+} \right)}^{*}} \right]\left[ {{C}_{3}}\varphi _{\beta }^{-} \right]{{{\rm e}}^{-{\rm i}{\boldsymbol{g}} \cdot C_{3}^{-1}{\boldsymbol{r}} }}\;{\rm d}{\boldsymbol{r}} \\ {} & = & {{\left( \gamma _{\alpha }^{+} \right)}^{*}}\gamma _{\beta }^{-}\frac{1}{\Omega }{{\int }_{\Omega }}{{\left( \varphi _{\alpha }^{+} \right)}^{*}}\varphi _{\beta }^{-}{{{\rm e}}^{-{\rm i}{{C}_{3}}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }}\;{\rm d}{\boldsymbol{r}} \\ {} & = & {{\left( \gamma _{\alpha }^{+} \right)}^{*}}\gamma _{\beta }^{-}{{\mathbb{S}}_{\alpha \beta }}({{C}_{3}}{\boldsymbol{g}} ), \\ \end{array}\end{eqnarray} \tag{ 14 }$$

where $\gamma _{\alpha }^{\tau }\equiv [{{C}_{3}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )]/\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )$ is the eigenvalue of the C₃ operation corresponding to the Bloch function $\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )$.

For the integral ${{I}_{\alpha \beta }}$, we also rewrite it as ${{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} )\equiv {{I}_{\alpha \beta }}(\Delta {\boldsymbol{G}} -2{\boldsymbol{K}} )$. Using the integral form of the Bessel function of order n (integer) ${{J}_{n}}(x)=\frac{1}{2\pi }\int _{-\pi }^{\pi }{{{\rm e}}^{{\rm i}(x{\rm sin} \phi -n\phi )}}\;{\rm d}\phi$, we can simplify ${{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} )$ to single integrals as

$$\begin{eqnarray}&&{{\mathbb{I}}_{{\rm cc}}}({\boldsymbol{g}} )=2\pi \int _{0}^{\infty }{{\left( b_{{\rm c}}^{+} \right)}^{*}}U(r)b_{{\rm c}}^{-}{{J}_{0}}(|{\boldsymbol{g}} |r)r\;{\rm d}r,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\mathbb{I}}_{{\rm cv}}}({\boldsymbol{g}} )=-2\pi {{{\rm e}}^{-{\rm i}\phi }}\int _{0}^{\infty }{{\left( b_{{\rm c}}^{+} \right)}^{*}}U(r)b_{{\rm v}}^{-}{{J}_{1}}(|{\boldsymbol{g}} |r)r\;{\rm d}r,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\mathbb{I}}_{{\rm vc}}}({\boldsymbol{g}} )=2\pi {{{\rm e}}^{-{\rm i}\phi }}\int _{0}^{\infty }{{\left( b_{{\rm v}}^{+} \right)}^{*}}U(r)b_{{\rm c}}^{-}{{J}_{1}}(|{\boldsymbol{g}} |r)r\;{\rm d}r,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{{\mathbb{I}}_{{\rm vv}}}({\boldsymbol{g}} )=-2\pi {{e}^{-i2\phi }}\int _{0}^{\infty }(b_{{\rm v}}^{+})U(r)b_{{\rm v}}^{-}{{J}_{2}}(|{\boldsymbol{g}} |r)r\;{\rm d}r,\end{eqnarray} \tag{ 18 }$$

where ϕ is the polar angle of the 2D vector ${\boldsymbol{g}}$. It is obvious that under ${\boldsymbol{g}} \to {{C}_{3}}{\boldsymbol{g}}$, the integral parts of equations (15)–(18) do not change, and hence the ratio ${{\eta }_{\alpha \beta }}\equiv {{\mathbb{I}}_{\alpha \beta }}({{C}_{3}}{\boldsymbol{g}} )/{{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} )$ is determined only by the prefactor in front of the integrals:

$$\begin{eqnarray}&&{{\eta }_{{\rm cc}}}=1,\ \ \ \ {{\eta }_{{\rm cv}}}={{\eta }_{{\rm vc}}}={{{\rm e}}^{-{\rm i}\frac{2\pi }{3}}},\ \ \ \ {{\eta }_{{\rm vv}}}={{{\rm e}}^{{\rm i}\frac{2\pi }{3}}}.\end{eqnarray} \tag{ 19 }$$

Because the set of ${\boldsymbol{g}}$, $\mathbb{G}$, also has the C₃ symmetry, we can choose one third elements of $\mathbb{G}$, denoted as ${{\mathbb{G}}_{\frac{1}{3}}}$, such that $\mathbb{G}=\sum _{n=0}^{2}C_{3}^{n}{{\mathbb{G}}_{\frac{1}{3}}}$. Using equations (14) and (19), we have

$$\begin{eqnarray}\begin{array}{rcl} {{V}_{\alpha \beta }} & = & \mathop{\sum }\limits_{{\boldsymbol{g}} \in {{\mathbb{G}}_{\frac{1}{3}}}}\mathop{\sum }\limits_{n=0}^{2}\left[ {{\mathbb{S}}_{\alpha \beta }}(C_{3}^{n}{\boldsymbol{g}} ){{\mathbb{I}}_{\alpha \beta }}(C_{3}^{n}{\boldsymbol{g}} ) \right] \\ {} & = & \mathop{\sum }\limits_{{\boldsymbol{g}} \in {{\mathbb{G}}_{\frac{1}{3}}}}\mathop{\sum }\limits_{n=0}^{2}\left\{ {{\left[ {{\left( \gamma _{\alpha }^{+} \right)}^{*}}\gamma _{\beta }^{-} \right]}^{-n}}{{\mathbb{S}}_{\alpha \beta }}({\boldsymbol{g}} )\eta _{\alpha \beta }^{n}{{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} ) \right\} \\ {} & = & (1+{{\lambda }_{\alpha \beta }}+\lambda _{\alpha \beta }^{2})\mathop{\sum }\limits_{{\boldsymbol{g}} \in {{\mathbb{G}}_{\frac{1}{3}}}}{{\mathbb{S}}_{\alpha \beta }}({\boldsymbol{g}} ){{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} ), \\ \end{array}\end{eqnarray} \tag{ 20 }$$

where

$$\begin{eqnarray}&&{{\lambda }_{\alpha \beta }}={{\left[ {{\left( \gamma _{\alpha }^{+} \right)}^{*}}\gamma _{\beta }^{-} \right]}^{-1}}{{\eta }_{\alpha \beta }}.\end{eqnarray} \tag{ 21 }$$

${{\eta }_{\alpha \beta }}$ is already given in equation (19) and the rest is to determine $\gamma _{\alpha }^{\tau }$.

For the MoS₂ monolayer, the conduction and valence band edge Bloch functions at $\pm K$ points are formed predominantly by the d₀ (i.e. ${{d}_{{{z}^{2}}}}$) and ${{d}_{\pm 2}}$ (i.e. $\frac{1}{\sqrt{2}}({{d}_{{{x}^{2}}-{{y}^{2}}}}\pm {\rm i}{{d}_{xy}})$) orbitals of Mo atoms respectively [24, 44, see footnote 3]. There are two contributions to the eigenvalue $\gamma _{\alpha }^{\tau }$ of the C₃ rotation on the Bloch function $\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )={{{\rm e}}^{{\rm i}\tau {\boldsymbol{K}} \cdot {\boldsymbol{r}} }}u_{\alpha }^{\tau }({\boldsymbol{r}} )$ (details in appendix A): (i) the rotation of the atomic orbital around its own center, i.e. ${{C}_{3}}{{d}_{0}}({\boldsymbol{r}} )={{d}_{0}}({\boldsymbol{r}} )$ and ${{C}_{3}}{{d}_{\pm 2}}({\boldsymbol{r}} )={{{\rm e}}^{\pm }}{\rm i}\frac{2\pi }{3}{{d}_{\pm 2}}({\boldsymbol{r}} )$; (ii) the change of lattice phase, defined as the value of the planewave phase factor ${{{\rm e}}^{{\rm i}\tau {\boldsymbol{K}} \cdot {\boldsymbol{r}} }}$ at each lattice site (see figure 1 for the case of $\tau =+1$). The change of the lattice phase factor under C₃ rotation depends on the rotation center. If the rotation center is at a Mo atom as shown in figure 1(a), the C₃ rotation does not change the lattice phase factors leaving only contribution (i) taking effect, and we obtain

$$\begin{eqnarray}&&\gamma _{{\rm c}}^{\tau }=1\ \ \ {\rm and}\ \ \ \gamma _{{\rm v}}^{\tau }={{{\rm e}}^{{\rm i}\tau \frac{2\pi }{3}}}\ \ \ \ ({\rm for}\ {\rm a}\ {\rm Mo}\ {\rm center}).\end{eqnarray} \tag{ 22 }$$

However, if the rotation center is at a S atom, as shown in figure 1(b), both contributions (i) and (ii) take effect and we have

$$\begin{eqnarray}&&\gamma _{{\rm c}}^{\tau }={{{\rm e}}^{{\rm i}\tau \frac{2\pi }{3}}}\ \ \ {\rm and}\ \ \ \gamma _{{\rm v}}^{\tau }={{{\rm e}}^{-{\rm i}\tau \frac{2\pi }{3}}}\ \ \ \ ({\rm for}\ {\rm a}\ {\rm S}\ {\rm center}).\end{eqnarray} \tag{ 23 }$$

Putting equations (19) and (22) (or (23)) into equation (21), we obtain that ${{\lambda }_{\alpha \beta }}=1$ for Mo centered potential and ${{\lambda }_{\alpha \beta }}={{{\rm e}}^{-{\rm i}\frac{2\pi }{3}}}$ for S centered potential, and the final intervalley coupling is

$$\begin{eqnarray}{{V}_{{\rm inter}}}=\left\{ \begin{array}{ccccccccccccccc} 3\mathop{\sum }\limits_{\alpha \beta }\mathop{\sum }\limits_{{\boldsymbol{g}} \in {{\mathbb{G}}_{\frac{1}{3}}}}{{\mathbb{S}}_{\alpha \beta }}({\boldsymbol{g}} ){{\mathbb{I}}_{\alpha \beta }}({\boldsymbol{g}} ), & {\rm for}\ {\rm Mo}\ {\rm centered}\ {\rm potential} \\ 0, & {\rm for}\ {\rm S}\ {\rm centered}\ {\rm potential} \\ \end{array} \right..\end{eqnarray} \tag{ 24 }$$

From the above symmetry analysis, we conclude that the intervalley coupling induced by a circular confinement potential is sensitive to the center of the potential. The intervalley coupling is enhanced when the potential is centered at a Mo atom, while it vanishes if the potential is centered at a S atom. This is due to the dependence of $\gamma _{\alpha }^{\tau }$, the eigenvalue of the Bloch function $\varphi _{\alpha }^{\tau }$ under the C₃ rotation, on the location of the rotation center. Below, we show that the same conclusions can be drawn as long as the confinement potential has the C₃ symmetry, i.e. intervalley coupling is strongest (zero) if the center of the potential is at a Mo (S) site.

3.1. Noncircular QD with C₃ symmetry

Consider a noncircular QD confinement potential with C₃ symmetry only, ${{C}_{3}}U({\boldsymbol{r}} )=U({\boldsymbol{r}} )$. In such a potential which lacks the circular symmetry, we do not have an exact solution in general for the massive Dirac fermion Hamiltonian in equation (3). Nevertheless, the numerical results presented in section 2.2 have justified well the effective mass approximation, where we can construct the QD electron wavefunction based on the conduction band edge Bloch function only. The ground-state envelope function $\psi _{{\rm c}}^{\tau }$ can be obtained by solving the following effective-mass Hamiltonian

$$\begin{eqnarray}&&{{H}_{{\rm em}}}=-\frac{{{\hbar }^{2}}}{2{{m}_{{\rm eff}}}}{{\nabla }^{2}}+U({\boldsymbol{r}} ),\end{eqnarray} \tag{ 25 }$$

where ${{m}_{{\rm eff}}}$ is the effective mass of conduction band at $\pm K$. With the same ${{m}_{{\rm eff}}}$ at $\pm K$, the envelope functions in the two valleys are the same: $\psi _{{\rm c}}^{+}=\psi _{{\rm c}}^{-}$. To analyze the symmetry, we need to know $\gamma _{{\rm c}}^{\tau }$ and ${{\eta }_{{\rm cc}}}$ (see equations (20) and (21)). $\gamma _{{\rm c}}^{\tau }$ is independent of the shape of the QD, and equations (22) and (23) are still valid here. ${{\eta }_{{\rm cc}}}$ is determined by the symmetry of the envelope function $\psi _{{\rm c}}^{\tau }$. Since $U({\boldsymbol{r}} )$ is symmetric under C₃ rotation, the non-degenerate ground-state envelope function $\psi _{{\rm c}}^{\tau }$ in each valley has to be an eigenstate of C₃, i.e. ${{C}_{3}}\psi _{{\rm c}}^{\tau }=\rho \psi _{{\rm c}}^{\tau }$ where $|\rho |=1$. Using equation (10) and ${{\rho }^{*}}\rho =1$, we have

$$\begin{eqnarray}\begin{array}{rcl} {{\mathbb{I}}_{{\rm cc}}}({{C}_{3}}{\boldsymbol{g}} ) & = & \int \psi _{{\rm c}}^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\psi _{{\rm c}}^{-}({\boldsymbol{r}} ){{{\rm e}}^{{\rm i}{{C}_{3}}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }}U({\boldsymbol{r}} ){\rm d}{\boldsymbol{r}} =\int {{\left[ \rho \psi _{{\rm c}}^{+}({\boldsymbol{r}} ) \right]}^{*}}[\rho \psi _{{\rm c}}^{-}({\boldsymbol{r}} )]{{{\rm e}}^{{\rm i}{\boldsymbol{g}} \cdot C_{3}^{-1}{\boldsymbol{r}} }}U({\boldsymbol{r}} ){\rm d}{\boldsymbol{r}} \\ {} & = & \int {{\left[ {{C}_{3}}\psi _{{\rm c}}^{+}({\boldsymbol{r}} ) \right]}^{*}}[{{C}_{3}}\psi _{{\rm c}}^{-}({\boldsymbol{r}} )]{{{\rm e}}^{{\rm i}{\boldsymbol{g}} \cdot C_{3}^{-1}{\boldsymbol{r}} }}[{{C}_{3}}U({\boldsymbol{r}} )]\;{\rm d}{\boldsymbol{r}} =\int {{C}_{3}}[\psi _{{\rm c}}^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\psi _{{\rm c}}^{-}({\boldsymbol{r}} ){{{\rm e}}^{{\rm i}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }}U({\boldsymbol{r}} )]\;{\rm d}{\boldsymbol{r}} \\ {} & = & \int \psi _{{\rm c}}^{+}{{\left( {\boldsymbol{r}} \right)}^{*}}\psi _{{\rm c}}^{-}({\boldsymbol{r}} ){{{\rm e}}^{{\rm i}{\boldsymbol{g}} \cdot {\boldsymbol{r}} }}U({\boldsymbol{r}} ){\rm d}{\boldsymbol{r}} ={{\mathbb{I}}_{{\rm cc}}}({\boldsymbol{g}} ) \\ \end{array}\end{eqnarray} \tag{ 26 }$$

Thus ${{\eta }_{{\rm cc}}}=1$, consistent with equation (19). The conclusion of equation (24) is therefore valid for a QD confinement potential with C₃ symmetry.

In summary, the band edge Bloch functions at K and $-K$ have the C₃ rotational symmetry of the lattice. However, they transform differently when the rotation center is a Mo site or a S site, as shown in figures 1(a) and (b), respectively. When the QD confinement potential also has the C₃ rotational symmetry about a S site, the rotational symmetry of the band edge Bloch functions about the S site leads to destructive interference and hence vanishing intervalley coupling. This is the physical origin of the dependence of intervalley coupling on the central position of QD potential.

In the EFM, except for the circular shaped confinement potential, it is difficult to find the exact wavefunction of the bound states in the QD. To investigate the intervalley coupling in the more general case of the QDs with noncircular shape, we consider here the alternative approach of the RSTB method. Using RSTB, we first calculate intervalley couplings of circular shaped QDs and compare these with the results from EFM. The agreement between the two methods justifies well the validity of the RSTB method. Then we calculate intervalley couplings of noncircular shaped QDs and analyze the results. Since SOC has negligible effects on intervalley coupling, we did not take into account SOC in the RSTB calculation.

4.1. The RSTB model

We use supercell and periodic boundary condition to model a QD. The supercell is a square with side length ${{L}_{{\rm sc}}}$ (see figure 9). Each supercell has $N={{L}_{{\rm sc}}}*{\rm round}\;(2{{L}_{{\rm sc}}}/\sqrt{3})$ lattice points for primitive cells. All lengths are in unit of lattice constant a. We adopt the three-band tight-binding model developed in [44], considering hoppings between nearest-neighbour Mo-${{d}_{{{z}^{2}}}}$, d_xy, and ${{d}_{{{x}^{2}}-{{y}^{2}}}}$ orbitals. The coordinates of the lattice sites of Mo atoms are denoted as ${{{\boldsymbol{r}} }_{n}}$ $(n=1,2,\cdots ,N)$. The tight-binding Hamiltonian of the QD is

$$\begin{eqnarray}&&H=\mathop{\sum }\limits_{i=1}^{N}\mathop{\sum }\limits_{\alpha }\left[ {{\epsilon }_{\alpha }}+U({{{\boldsymbol{r}} }_{i}}) \right]c_{i\alpha }^{\dagger }{{c}_{i\alpha }}+\mathop{\sum }\limits_{\langle i,j\rangle }\mathop{\sum }\limits_{\alpha ,\beta }{{t}_{i\alpha ,j\beta }}c_{i\alpha }^{\dagger }{{c}_{j\beta }},\end{eqnarray} \tag{ 27 }$$

where $U({{{\boldsymbol{r}} }_{i}})$ is the QD confinement potential at ${{{\boldsymbol{r}} }_{i}}$, $\alpha ,\beta ={{d}_{{{z}^{2}}}},\;{{d}_{xy}},\;{{d}_{{{x}^{2}}-{{y}^{2}}}}$ are orbital indices, $c_{i\alpha }^{\dagger }$ is the creation operator for electron on site i with orbital α. ${{t}_{i\alpha ,j\beta }}$ is the hopping between α orbital at position i and β orbital at position j, and $\langle ,\rangle$ denotes summation over nearest-neighbor pairs only. In practice, we write the Hamiltonian H in form of a $3N\times 3N$ matrix. The part of H involving hoppings between site i and j can be written as a 3 × 3 block ${{H}_{ij}}=h({{{\boldsymbol{r}} }_{j}}-{{{\boldsymbol{r}} }_{i}})$, where ${{{\boldsymbol{r}} }_{j}}-{{{\boldsymbol{r}} }_{i}}$ is either ${\boldsymbol{0}}$ or the six nearest-neighbour vectors ${{{\boldsymbol{R}} }_{1}}$ to ${{{\boldsymbol{R}} }_{6}}$ defined in [44]. Using the real hopping parameters $\{{{\epsilon }_{1}},{{\epsilon }_{2}},{{t}_{0}},{{t}_{1}},{{t}_{2}},{{t}_{11}},{{t}_{12}},{{t}_{22}}\}$ defined in [44] (the GGA-version parameters therein), we have,

$$\begin{eqnarray}h({\boldsymbol{0}} )=\left[ \begin{array}{ccccccccccccccc} {{\epsilon }_{1}} & 0 & 0 \\ 0 & {{\epsilon }_{2}} & 0 \\ 0 & 0 & {{\epsilon }_{2}} \\ \end{array} \right],\ \ \ \ \ \ \ h({{{\boldsymbol{R}} }_{1}})=\left[ \begin{array}{ccccccccccccccc} {{t}_{0}} & {{t}_{1}} & {{t}_{2}} \\ -{{t}_{1}} & {{t}_{11}} & {{t}_{12}} \\ {{t}_{2}} & -{{t}_{12}} & {{t}_{22}} \\ \end{array} \right],\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}h({{{\boldsymbol{R}} }_{2}})=\left[ \begin{array}{ccccccccccccccc} {{t}_{0}} & \frac{1}{2}\left( {{t}_{1}}-\sqrt{3}{{t}_{2}} \right) & -\frac{1}{2}\left( \sqrt{3}{{t}_{1}}+{{t}_{2}} \right) \\ -\frac{1}{2}\left( {{t}_{1}}+\sqrt{3}{{t}_{2}} \right) & \frac{1}{4}\left( {{t}_{11}}+3{{t}_{22}} \right) & -\frac{\sqrt{3}}{4}({{t}_{11}}-{{t}_{22}})-{{t}_{12}} \\ \frac{1}{2}\left( \sqrt{3}{{t}_{1}}-{{t}_{2}} \right) & -\frac{\sqrt{3}}{4}({{t}_{11}}-{{t}_{22}})+{{t}_{12}} & \frac{1}{4}\left( 3{{t}_{11}}+{{t}_{22}} \right) \\ \end{array} \right],\quad \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}h({{{\boldsymbol{R}} }_{3}})=\left[ \begin{array}{ccccccccccccccc} {{t}_{0}} & -\frac{1}{2}\left( {{t}_{1}}-\sqrt{3}{{t}_{2}} \right) & -\frac{1}{2}\left( \sqrt{3}{{t}_{1}}+{{t}_{2}} \right) \\ \frac{1}{2}\left( {{t}_{1}}+\sqrt{3}{{t}_{2}} \right) & \frac{1}{4}\left( {{t}_{11}}+3{{t}_{22}} \right) & \frac{\sqrt{3}}{4}({{t}_{11}}-{{t}_{22}})+{{t}_{12}} \\ \frac{1}{2}\left( \sqrt{3}{{t}_{1}}-{{t}_{2}} \right) & \frac{\sqrt{3}}{4}({{t}_{11}}-{{t}_{22}})-{{t}_{12}} & \frac{1}{4}\left( 3{{t}_{11}}+{{t}_{22}} \right) \\ \end{array} \right],\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&h({{{\boldsymbol{R}} }_{4}})=h{{\left( {{{\boldsymbol{R}} }_{1}} \right)}^{\dagger }},h({{{\boldsymbol{R}} }_{5}})=h{{({{{\boldsymbol{R}} }_{2}})}^{\dagger }},\ \ \ \ h({{{\boldsymbol{R}} }_{6}})=h{{\left( {{{\boldsymbol{R}} }_{3}} \right)}^{\dagger }}.\end{eqnarray} \tag{ 31 }$$

The $3N\times 3N$ RSTB Hamiltonian matrix for a QD is then diagonalized to find the eigenstates and eigenenergies.

Zoom In Zoom Out Reset image size

Unlike the EFM approach where the bound state wavefunctions are first given in each valley, and the intervalley coupling is then calculated as the off-diagonal matrix element between them (see equation (2)), the eigenstates and eigenenergies of the RSTB Hamiltonian H already include the effects of intervalley coupling. The intervalley coupling causes a fine splitting of the energy levels, which otherwise have the two-fold valley degeneracy in the absence of SOC. This energy splitting is just $2|{{V}_{{\rm inter}}}|$. Thus, the magnitude of the intervalley coupling ${{V}_{{\rm inter}}}$ is read out from the difference between the lowest two conduction-band energy eigenvalues of H in the RSTB solution⁴ .

4.2. Circular shaped QDs

Here we study circular shaped QDs using RSTB and compare with the results in EFM. First, we check the impact of supercell size on the intervalley coupling. Figure 6 shows the $|{{V}_{{\rm inter}}}|$ versus ${{L}_{{\rm sc}}}$ relation for different circular QDs. We can see that $|{{V}_{{\rm inter}}}|$ converges fast with increasing ${{L}_{{\rm sc}}}$, which means that the finite-size effects introduced by the method can be eliminated well by having a large enough supercell. Figure 6 also shows that $|{{V}_{{\rm inter}}}|$ converges faster with increasing ${{L}_{{\rm sc}}}$ for larger potential depth V, just as expected.

Zoom In Zoom Out Reset image size

$|{{V}_{{\rm inter}}}|$ calculated using the RSTB method with a large supercell (${{L}_{{\rm sc}}}=240$) is compared with the results from EFM in figure 7. For the $|{{V}_{{\rm inter}}}|$ versus R dependence (figure 7(a)), the $|{{V}_{{\rm inter}}}|$ versus L dependence (figure 7(b)), and the $|{{V}_{{\rm inter}}}|$ versus V dependence (figure 7(c)), the RSTB data, including the rapid local oscillations, agree quantitatively well with the EFM ones. The EFM results rely on the details of the first-principles wavefunctions of Bloch states from the ABINIT package, while the RSTB calculations are based on a few hopping matrix elements only. This agreement between the two entirely different approaches well justify both the RSTB method and EFM.

Zoom In Zoom Out Reset image size

In figure 8, we analyze the dependence of the intervalley coupling strength on the central position of the confinement potential. Figure 8(a) shows $|{{V}_{{\rm inter}}}|$ versus R curves of the same circular shaped potential centered at a Mo or S site. Consistent with the symmetry analysis in section 3, $|{{V}_{{\rm inter}}}|$ is three orders of magnitude larger for the Mo-centered case compared to the S-centered case, the latter can be regarded as the numerical errors from zero. Figure 8(b) plots $|{{V}_{{\rm inter}}}|$ as a function of the central position of the confinement potential. Clearly, $|{{V}_{{\rm inter}}}|$ decreases when the potential center moves away from the Mo site and eventually vanishes when the potential center is at the S site.

Zoom In Zoom Out Reset image size

4.3. Intervalley coupling in noncircular QD potentials with vertical wall (L = 0)

QDs can be formed by the lateral heterostructures between different monolayer TMDs. Currently such heterostructures are grown by chemical vapor deposition (CVD), where the shape of confinement is triangular [17]. The lateral size of the confinement in available heterostructures is several μm, still too large for QDs. Nevertheless, growth of monolayer TMDs flakes with sizes ranging from tens of nm to tens of μm and with both triangular and hexagonal shapes have been reported by various groups [55–62]. Smaller QDs formed by the lateral heterostructures can be expected. In such heterostructures, the confinement potential is due to the band edge difference between the TMDs. Thus the potential well has a vertical wall at the QD boundary, corresponding to the L = 0 limit discussed above. Here we use the RSTB method to investigate intervalley coupling of these noncircular QDs. We studied QDs with triangular, hexagonal as well as square shapes for comparison.

As shown in figure 9, the triangular, hexagonal, and square QDs are defined by the confinement potentials ${{U}_{{\rm T}}}({\boldsymbol{r}} )$, ${{U}_{{\rm H}}}({\boldsymbol{r}} )$, and ${{U}_{{\rm S}}}({\boldsymbol{r}} )$, respectively, where ${{U}_{{\rm T}/{\rm H}/{\rm S}}}$ equals 0 within the QDs (yellow regions), and equals the constant V outside (white regions). The 'radius' of the QDs is denoted as R (see figure 9). For the triangular and hexagonal QDs, we have chosen the orientation of the confinement potential such that the sides are all along the zigzag crystalline axes, corresponding to those realizable by the lateral heterostructures of different dichalcogenides. For the square QDs, two sides are along the zigzag and the others are along the armchair crystalline axes.

Figure 10 plots the intervalley coupling as function of the QDs size, for three different potential depths $V=0.01\;$ eV, $V=0.1\;$ eV and $V=0.2\;$ eV. With increasing R, the intervalley coupling strength has an overall trend to decrease, but with rapid local oscillations. These features are similar to the circular QDs (see figure 3). Interestingly, for the triangular and hexagonal QDs, intervalley coupling is found to be much smaller with the larger potential jump V at the QD boundary. In figure 11, the dependence of the intervalley coupling on V is investigated, in which all curves are non-monotonic. The non-monotonic dependence of intervalley coupling strength is also found for the circular QDs when L = 0 (see figure 5). We note that in heterostructures formed between different monolayer dichalcogenides (e.g. MoSe₂ and WSe₂), the conduction and valence band edge discontinuities are discovered to be in the range of 0.2–0.4 eV [18]. For this range of values for V, the intervalley coupling is negligible for the triangular and hexagonal QDs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Comparing QDs of different shapes, the intervalley coupling in square QDs is orders of magnitude larger than that in the triangular and hexagonal QDs. Such a sharp difference is due to the relative orientation of the sides of the confinement potential to the crystalline axes. Inhomogeneous junctions along the zigzag (armchair) crystalline axes preserves the momentum along (perpendicular to) the lines connecting the two valleys in the momentum space. Thus, the zigzag junctions introduce the minimum intervalley coupling, while the armchair junctions cause maximum intervalley coupling. For the triangular and hexagonal QDs, all sides are along the zigzag crystalline axes, while for square QDs two sides are along the armchair directions and thus intervalley coupling is much larger. This is further verified by examining the change of the intervalley coupling when we rotate the QD confinement potential by an angle ${{\theta }_{{\rm rot}}}$ relative to the orientation defined in figure 9, as shown in figure 12. For triangular and hexagonal QDs, we find the intervalley coupling increases with ${{\theta }_{{\rm rot}}}$ significantly, and reaches a maximum at ${{\theta }_{{\rm rot}}}={{30}^{{}^\circ }}$ where the intervalley coupling becomes comparable with the square QDs. Note that at ${{\theta }_{{\rm rot}}}={{30}^{{}^\circ }}$, the sides of the confinement potential of the triangular and hexagonal QDs are along the armchair crystalline axes (see figures 12(a)and (b)).

Zoom In Zoom Out Reset image size

Figure 13 shows the dependence of $|{{V}_{{\rm inter}}}|$ on the central position of the confinement potential in a unit cell for the three kinds of noncircular QDs. It is clear that $|{{V}_{{\rm inter}}}|$ strongly depends on the potential center. For triangular and hexagonal QDs with the C₃ rotational symmetry, the dependence is similar to the circular QDs (see figure 8(b)), i.e. intervalley coupling is strongest if the potential is centered at a Mo site, but vanishes if the potential is centered at a S site. Such behavior is absent for the square QDs that lack the C₃ rotational symmetry.

Zoom In Zoom Out Reset image size

4.4. Intervalley coupling in a smooth noncircular QD potential

Finally, we examine square shaped QDs but with a smooth slope for the potential change on the boundary (i.e. finite L), which represents a general situation for QDs defined by patterned electrodes. The calculated intervalley coupling strength $|{{V}_{{\rm inter}}}|$ as a function of QD size R is shown in figure 14. The peak value of intervalley coupling decreases nearly exponentially with the QD size. The dependence on the smoothness L is similar to the circular shaped QDs (see figure 4), i.e. $|{{V}_{{\rm inter}}}|$ drops fast by several orders of magnitude when L increases from 0 to 5, and becomes largely independent of L for $L\gt 5$. From figure 14, we can see that the $|{{V}_{{\rm inter}}}|$ curves as functions of R for L = 10, L = 20, and L = 30 are basically the same, except for the location of some fine dips.

Zoom In Zoom Out Reset image size

To conclude, we have investigated intervalley coupling in QDs defined by confinement potentials of various shapes on an extended TMD monolayer. The numerical results obtained using two completely different approaches agree well with each other. For confinement potentials with the C₃ rotational symmetry, the intervalley coupling is maximized (zero) if the potential is centered at a M (X) site. Comparing smooth confinement potentials with slopping walls and sharp confinement potentials with vertical walls, the intervalley coupling in the latter is much stronger. When the length scales of the QDs vary, the intervalley coupling exhibits fast oscillations where the maxima can be two orders of magnitude larger than the neighboring minima, and the envelop of the maxima has an overall nearly exponential decay with the increase of the QDs size. When all parameters are comparable, the intervalley coupling can be smaller by several orders of magnitude in circular shaped QDs, and in triangular and hexagonal shaped QDs with all sides along the zigzag crystalline axes. For all QDs studied, the largest intervalley coupling is upper bounded by 0.1 meV, found for small QDs with diameters of 20 nm and with sharp confinement potentials (i.e. vertical walls).

The intervalley coupling in these monolayer TMD QDs is much smaller compared to that in graphene nanoribbon QDs [12] and in silicon QDs [32–39]. In QDs formed in 2D crystals, intervalley coupling mainly occurs at the boundary where the translational invariance is lost, thus the coupling strength is proportional to the probability distribution of the electron at the QD boundary. In the TMD QDs studied here, the electron wavefunction has vanishing amplitude at the QD boundary (see green shaded region in figure 2(b)). This is a characteristic of the QDs defined by confinement potentials on an otherwise extended crystalline monolayer, as opposed to the nanoribbon QDs [12]. The geometry here is also qualitatively different from the QDs formed in silicon inversion layers. In those silicon QDs, the two valleys are separated by a wavevector along the z-direction (perpendicular to the inversion layer), and their coupling is caused by the sharp confinement with a length scale of nm in the z-direction [32–39]. Such confinement strongly squeezes the wavefunction, giving rise to its large amplitude at the interface where valley hybridization occurs, and hence there is a stronger intervalley coupling in the order of meV. In contrast, here it is the much larger lateral size R of the QDs that determines the wavefunction amplitude at the QD boundary. We can expect that intervalley coupling is generically small in QDs defined by confinement potentials on otherwise extended monolayer.

Because of the very small intervalley coupling in the monolayer TMD QDs, valley hybridization is then well quenched by the much stronger spin–valley coupling of the electron in monolayer TMDs. As a result, the QDs can easily inherit the valley physics of the 2D bulk such as the valley optical selection rules, making the valley pseudospin of a single confined electron a perfect candidate for a quantum bit. The sensitive dependence of intervalley coupling strength on the central position and the lateral length scales of the confinement potentials further makes possible the tuning of the intervalley coupling by external controls. The qualitative conclusions here are also applicable to QDs formed by confinement potentials on other 2D materials with finite gap.

The work was supported by the Croucher Foundation under the Croucher Innovation Award, the Research Grant Council of HKSAR under Grant No. HKU705513P and HKU9/CRF/13G (GBL, HP, and WY); the NSFC with Grant No. 11304014, the National Basic Research Program of China 973 Program with Grant No. 2013CB934500 and the Basic Research Funds of Beijing Institute of Technology with Grant No. 20131842001 and 20121842003 (GBL.); the MOST Project of China with Grants Nos. 2014CB920903 and 2011CBA00100, the NSFC with Grant Nos. 11174337 and 11225418, the Specialized Research Fund for the Doctoral Program of Higher Education of China with Grants No. 20121101110046 (YY).

Here we give a derivation of $\gamma _{\alpha }^{\tau }$, the eigenvalue of the Bloch function $\varphi _{\alpha }^{\tau }$ under the C₃ rotation. We first express $\varphi _{\alpha }^{\tau }$ in terms of linear combination of atomic orbitals ${{d}_{\alpha }}$ ($\alpha =$ c, v)

$$\begin{eqnarray}&&\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )=\frac{1}{\sqrt{N}}\mathop{\sum }\limits_{{\boldsymbol{R}} }{{{\rm e}}^{{\rm i}\tau K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}{{d}_{\alpha }}({\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} ),\end{eqnarray} \tag{ A.1 }$$

in which ${\boldsymbol{R}}$ is lattice vector, ${\boldsymbol{ \delta }}$ is the position of Mo atom in the unit cell, ${{d}_{{\rm c}}}={{d}_{0}}={{d}_{{{z}^{2}}}}$, and ${{d}_{{\rm v}}}={{d}_{\tau 2}}=\frac{1}{\sqrt{2}}({{d}_{{{x}^{2}}-{{y}^{2}}}}+{\rm i}\tau {{d}_{xy}})$. According the definition of $\gamma _{\alpha }^{\tau }\equiv [{{C}_{3}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )]/\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )$ and using the orthogonality of C₃ (${\boldsymbol{k}} \cdot {\boldsymbol{r}} ={{C}_{3}}{\boldsymbol{k}} \cdot {{C}_{3}}{\boldsymbol{r}}$), we have

$$\begin{eqnarray}\begin{array}{rcl} {{C}_{3}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} ) & = & \frac{1}{\sqrt{N}}\mathop{\sum }\limits_{{\boldsymbol{R}} }{{{\rm e}}^{{\rm i}\tau K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}{{d}_{\alpha }}(C_{3}^{-1}{\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} ) \\ {} & {} & =\frac{1}{\sqrt{N}}\mathop{\sum }\limits_{{\boldsymbol{R}} }{{{\rm e}}^{{\rm i}\tau {{C}_{3}}K\cdot {{C}_{3}}({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}{{d}_{\alpha }}(C_{3}^{-1}({\boldsymbol{r}} -{{C}_{3}}({\boldsymbol{R}} +{\boldsymbol{ \delta }} ))). \\ \end{array}\end{eqnarray} \tag{ A.2 }$$

Summing over all Mo positions $\{{\boldsymbol{R}} +{\boldsymbol{ \delta }} \}$ is equivalent to summing over all positions $\{{{C}_{3}}({\boldsymbol{R}} +{\boldsymbol{ \delta }} )\}$ when the rotation center locates at either a Mo or S site (see figure 1). Accordingly, with the substitution ${{C}_{3}}({\boldsymbol{R}} +{\boldsymbol{ \delta }} )\to ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )$, the above equation becomes

$$\begin{eqnarray}&&{{C}_{3}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )=\frac{1}{\sqrt{N}}\mathop{\sum }\limits_{{\boldsymbol{R}} }{{{\rm e}}^{{\rm i}\tau {{C}_{3}}K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}{{d}_{\alpha }}(C_{3}^{-1}({\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} )).\end{eqnarray} \tag{ A.3 }$$

Because C₃ is an element of the wave-vector group of $\tau K$, we have ${{{\rm e}}^{{\rm i}\tau {{C}_{3}}K\cdot {\boldsymbol{R}} }}={{{\rm e}}^{{\rm i}\tau K\cdot {\boldsymbol{R}} }}$ and consequently

$$\begin{eqnarray}&&{{{\rm e}}^{{\rm i}\tau {{C}_{3}}K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}={{{\rm e}}^{{\rm i}\tau K\cdot {\boldsymbol{R}} }}{{{\rm e}}^{{\rm i}\tau {{C}_{3}}K\cdot {\boldsymbol{ \delta }} }}={{{\rm e}}^{{\rm i}\tau K\cdot {\boldsymbol{R}} }}{{{\rm e}}^{{\rm i}\tau K\cdot C_{3}^{-1}{\boldsymbol{ \delta }} }}={{{\rm e}}^{{\rm i}\tau K\cdot (C_{3}^{-1}{\boldsymbol{ \delta }} -{\boldsymbol{ \delta }} )}}{{{\rm e}}^{{\rm i}\tau K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}.\end{eqnarray} \tag{ A.4 }$$

Define ${{\gamma }_{{{d}_{\alpha }}}}\equiv [{{C}_{3}}{{d}_{\alpha }}({\boldsymbol{r}} )]/{{d}_{\alpha }}({\boldsymbol{r}} )={{d}_{\alpha }}(C_{3}^{-1}{\boldsymbol{r}} )/{{d}_{\alpha }}({\boldsymbol{r}} )$, we have

$$\begin{eqnarray}&&{{d}_{\alpha }}(C_{3}^{-1}({\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} ))={{\gamma }_{{{d}_{\alpha }}}}{{d}_{\alpha }}({\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} ).\end{eqnarray} \tag{ A.5 }$$

Plugging equations (A.4) and (A.5) into (A.3), we obtain

$$\begin{eqnarray}&&{{C}_{3}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )=\frac{1}{\sqrt{N}}\mathop{\sum }\limits_{{\boldsymbol{R}} }{{{\rm e}}^{{\rm i}\tau K\cdot (C_{3}^{-1}{\boldsymbol{ \delta }} -{\boldsymbol{ \delta }} )}}{{{\rm e}}^{{\rm i}\tau K\cdot ({\boldsymbol{R}} +{\boldsymbol{ \delta }} )}}{{\gamma }_{{{d}_{\alpha }}}}{{d}_{\alpha }}({\boldsymbol{r}} -{\boldsymbol{R}} -{\boldsymbol{ \delta }} )={{{\rm e}}^{{\rm i}\tau K\cdot (C_{3}^{-1}{\boldsymbol{ \delta }} -{\boldsymbol{ \delta }} )}}{{\gamma }_{{{d}_{\alpha }}}}\varphi _{\alpha }^{\tau }({\boldsymbol{r}} )\end{eqnarray} \tag{ A.6 }$$

and therefore

$$\begin{eqnarray}&&\gamma _{\alpha }^{\tau }={{{\rm e}}^{{\rm i}\tau K\cdot (C_{3}^{-1}{\boldsymbol{ \delta }} -{\boldsymbol{ \delta }} )}}{{\gamma }_{{{d}_{\alpha }}}},\end{eqnarray} \tag{ A.7 }$$

in which the factor ${{\gamma }_{{{d}_{\alpha }}}}$ comes from the rotation of the atomic orbital around its own center and the factor ${{{\rm e}}^{{\rm i}\tau K\cdot (C_{3}^{-1}{\boldsymbol{ \delta }} -{\boldsymbol{ \delta }} )}}$ comes from the change of the lattice phase of Mo atom after rotation. According to equation (A.7), when the rotation center locates at Mo (figure 1(a)), we have ${\boldsymbol{ \delta }} =0$ and hence $\gamma _{\alpha }^{\tau }={{\gamma }_{{{d}_{\alpha }}}}$ (equation (22)); when the rotation center locates at S (figure 1(b)), we have ${\boldsymbol{ \delta }} =(1,\frac{1}{\sqrt{3}})a$ and hence $\gamma _{\alpha }^{\tau }={{{\rm e}}^{{\rm i}\tau \frac{2\pi }{3}}}{{\gamma }_{{{d}_{\alpha }}}}$ (equation (23)).