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

Monolayer transition metal dichalcogenides (TMDs) offer new opportunities for realizing quantum dots (QDs) in the ultimate two-dimensional (2D) limit. Given the rich control possibilities of electron valley pseudospin discovered in the monolayers, this quantum degree of freedom can be a promising carrier of information for potential quantum spintronics exploiting single electrons in TMD QDs. An outstanding issue is to identify the degree of valley hybridization, due to the QD confinement, which may significantly change the valley physics in QDs from its form in the 2D bulk. Here we perform a systematic study of the intervalley coupling by QD confinement potentials on extended TMD monolayers. We find that the intervalley coupling in such geometry is generically weak due to the vanishing amplitude of the electron wavefunction at the QD boundary, and hence valley hybridization shall be well quenched by the much stronger spin-valley coupling in monolayer TMDs and the QDs can well inherit the valley physics of the 2D bulk. We also discover sensitive dependence of intervalley coupling strength on the central position and the lateral length scales of the confinement potentials, which may possibly allow tuning of intervalley coupling by external controls

Monolayer transition metal dichalcogenides (TMDs) offer new opportunities for realizing quantum dots (QDs) in the ultimate two-dimensional (2D) limit. Given the rich control possibilities of electron valley pseudospin discovered in the monolayers, this quantum degree of freedom can be a promising carrier of information for potential quantum spintronics exploiting single electrons in TMD QDs. An outstanding issue is to identify the degree of valley hybridization, due to the QD confinement, which may significantly change the valley physics in QDs from its form in the 2D bulk. Here we perform a systematic study of the intervalley coupling by QD confinement potentials on extended TMD monolayers. We find that the intervalley coupling in such geometry is generically weak due to the vanishing amplitude of the electron wavefunction at the QD boundary, and hence valley hybridization shall be well quenched by the much stronger spin-valley coupling in monolayer TMDs and the QDs can well inherit the valley physics of the 2D bulk. We also discover sensitive dependence of intervalley coupling strength on the central position and the lateral length scales of the confinement potentials, which may possibly allow tuning of intervalley coupling by external controls.

I. INTRODUCTION
Semiconductor quantum dots (QDs) have been widely explored in the past several decades for technological applications such as lasers and medical markers [? ]. Spin and other quantum degrees of freedom of single electron or hole confined in QDs are also extensively researched as potential carriers of information for quantum computing and quantum spintronics [? ? ? ? ? ]. 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 [? ? ]. The emergence of graphene [? ? ], the 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 [? ? ? ], and graphene nanoribbon of armchair edges with patterned electrodes [? ]. The properties of these graphene QDs are affected significantly by the edges, and precise control on edge terminations are thus desired.
Monolayer group-VIB transition metal dichalcogenides (TMDs) are newly emerged members of the 2D crystal family [? ]. These TMDs have the chemical composition of MX 2 (M = Mo,W; X = S,Se), and the monolayer has a structure of X-M-X covalently bonded hexagonal quasi-2D network. The 2D bulk has a direct bandgap in the visible frequency range [? ? ], ideal for optoelectronic applications. This bandgap also makes possible the realization of QDs without the assistance of termination edges [? ]. 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 [? ]. 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 (cf. 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 alignment can form between different TMDs with different bandgaps and workfunctions. Lateral heterostructures with MoSe 2 islands surrounded by WSe 2 on a crystalline monolayer has been realized very recently [? ]. This can realize QD confinement, also with type-II band alignment, for electrons at the MoSe 2 region. The band edge discontinuity at the heterojunction, with an order of magnitude of 0. 2-0.4 eV [? ? ? ? ? ], forms the potential well (with a vertical wall).
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 Another remarkable property of the 2D TMDs is the strong interplay between spin and the valley pseudospin [? ? ? ]. These quantum degrees of freedom with versatile controllability well suggest that single electron in monolayer TMD QDs can be 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 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 [? ? ? ? ? ? ? ? ], 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 extended TMD monolayer. The numerical calculations used two methods: (i) the envelope function method (EFM) [? ? ? ? ], in conjunction with first-principles wavefunctions of the band-edge Bloch states of the 2D bulk; and (ii) the real-space tightbinding (RSTB) approach based on a three-band model for monolayer TMDs [? ]. 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 3 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 a M (X) site (cf. 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, 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 [? ]. In contrast to graphene QDs [? ] and silicon QDs [? ? ? ? ? ? ? ? ], the intervalley coupling here is much smaller as the electron wavefunction vanishes at the boundary of 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 [? ]. 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 II, we formulate the EFM for the two-band k·p model of the monolayer TMDs [? ], 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 III, we present the symmetry analysis on the Mcentered and X-centered potentials with C 3 rotational symmetry. Section IV 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 2 , and MoS 2 is presented as the example throughout the paper.

A. The model and method
We consider first a QD defined by a circular shaped potential U (r) in a TMD monolayer. Such 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 . ω ≡ e i 2π 3 , ω * and 1 are respectively the projection of the plane wave component of the Bloch function at K on each Mo site (i.e. e iK·r Mo ), referred as the lattice phase factor in the text. (c) The spatial profile of the conduction and valence band edges for the circular potential well with radius R and depth V . The length scale L characterizes the smoothness of the potential.
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 potential well with 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 ±K valleys of the 2D bulk. In the EFM, we start off with the two-band k · p Hamiltonian derived for the band edges at the ±K valleys of monolayer TMDs [? ? ], 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 Ψ τ,s where τ = ± is the valley index denoting the ±K valley, and s = +(↑) or −(↓) 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) [? ].
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 [? ? ], 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 Ψ +,s and Ψ −,s in the two valleys with the same spin index, which is defined as V s inter 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 Ψ −,s and Ψ +,s . If V s inter 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 s inter 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 s inter |. Due to the time-reversal symmetry, |V ↑ inter | = |V ↓ inter |. Hence we consider here only spin-up case and drop the spin index s in superscripts for brevity.
To obtain Ψ τ in EFM, we substitute −i∇ for k in the k · p Hamiltonian of the 2D bulk[? ? ], and the confinement potential U (r) is added as an onsite energy term. The QD Hamiltonian in each valley is then [? ? ? ? ? ? ? ] where ∆ ∼ 2 eV is the bulk bandgap of 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 ±K point, ϕ τ α (r) = e iτ K·r u τ α (r) (α =c, v and u τ α is the cell periodic part), and the eigenfunction of (3) is a spinor of envelope function which satisfy the Schrodinger equation H τ ψ τ = Eψ τ . The overall wavefunctions of the QD states are then Since U (r) has the rotational symmetry, it is easy to verify that the eigen spinor ψ τ is of the form where m is integer, and (r, θ) are the polar coordinates of r. b τ c and b τ v are real radial functions that satisfy , and need to be solved numerically in the region R − L 2 ≤ r ≤ R + L 2 . Boundary conditions are then matched at r = R − L 2 and r = R + L 2 to determine the solutions and eigenvalues. Bound states are found to exist in the energy range ( ∆ 2 , ∆ 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 Eq.
(2) can then be expressed as where and In the above equations, α, β ∈{c, v}, and c τ α (G) are the Fourier coefficients of the periodic part of the Bloch function u τ α (r) = G c τ α (G)e iG·r . G is a reciprocal lattice vector (RLV), and ∆G = G − 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)  We note that the Hamiltonian in equation (3) does not contain the conduction band spin-valley interaction [? ? ], which is of the Ising form that changes only the diagonal energies but not the wavefunctions of the states Ψ +,s and Ψ −,s . Thus, this spin-valley interaction does not affect the intervalley coupling strength defined here, but it leads to an energy difference between states Ψ +,s and Ψ −,s that ranges from several to several tens meV [? ? ]. Valley hybridization in the QDs is determined by the competition between |V inter | and this diagonal energy difference.

B. Numerical results
In this subsection we present the numerical calculation results. The parameters from [? ] are used: t = 1.1 eV, a = 3.193 Å, ∆ = 1.66 eV, and 2λ = 0.15 eV. We consider potentials centered at a Mo site, as shown in figure 1(a). As we will prove in Section III, 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.

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 [cf. equation (2)]. The energy spectra are shown in figure 2(a) for R = 80, L = 10 and three well depthsV =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 (τ = 1) and it is obvious from equation (7) (6)] for the ground states in (a), i.e. the lowest-energy ones with m = 0. Note that b + v (r) are magnified by a factor of 80 for clarity. The green shaded region corresponds to the slopping wall of the potential well. Units for energy and length are t and a respectively. Other symbols denote the contributions |V αβ | [αβ =cc (red cross), cv (blue triangle), vc (magenta square), and vv (green diamond)] to the intervalley coupling, for L = 0 and 10 respectively. See equations (8), (9), and (10) in the text for definitions. All lengths are in the unit of a and energies in the unit of t. and b + 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 + v (r) components are orders of magnitude smaller compared to b + 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) shall 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 [? ] and silicon QDs [? ? ? ? ? ? ? ? ] where significant valley hybridization is found.

The intervalley coupling
With the envelope functions given above, and the band-edge Bloch functions from the ABINIT package [? ? ], we can calculate the numerical values for the intervalley coupling [cf. equation (8)], Figure 3 shows the intervalley coupling strength |V 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 [cf. 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 discussed in detail later.   From figure 3(b) and (c), we can see that the dominant contribution to the intervalley coupling is from |V cc |, which is orders of magnitude larger than |V vv |, |V cv |, |V vc | [cf. equation (9)], consistent with the fact that the wavefunction of confined electron is predominantly from the conduction band edge Bloch functions. Figures 3(b) and (c) also clearly show that |V inter | oscillates vs 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 [? ? ? ]. Apart from these fast oscillations, the intervalley coupling strength decreases with the increase of R in a nearly exponential way. Figure 4 shows the |V inter | as a function of L, the smoothness of the potential, for circular QDs with R = 80 and V = 0.005. Largest intervalley coupling is found at L = 0 that corresponds to a circular square well, where the potential has a sudden jump at the QD boundary. It can be seen that |V 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 ≤ L ≤ 10, which shows that the oscillation is also in the length scale of lattice constant, similar to the |V inter | vs. R relation. Figure 5 shows the dependence of |V inter | on the potential depth V , when R = 80 and L = 0 and 10. We can see that |V 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 spin-orbit coupling (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 inter | vs. R curve (fixing L = 10, V = 0.005) and |V inter | vs. 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 inter |. This is well expected as |V 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 [? ], which does not affect the wavefunction. The SOC, however, will cause spin splittings and hence affect the spectral of quantum dot states in the K and −K valleys [cf. figure 2(a)]. When we consider the effect of intervalley coupling on the quantum dot states, the off-diagonal matrix element |V inter | will then compete with the differences in these diagonal energies, and the SOC induced spin splitting will play an important role.
From the dependence of |V 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 [? ? ? ? ? ? ? ? ].

III. SYMMETRY ANALYSIS FOR QDS WITH C3 ROTATIONAL SYMMETRY
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 inter | versus R or L. Here, based on symmetry analysis, we will show that the intervalley coupling strength |V inter | also depends on the position of the potential center, in a way that V inter is largest when the potential is centered at Mo atom [figure 1(a)] but vanishes if centered at S atom [ figure 1(b)]. To explain this, we analyze V αβ [see equation (9)] by examining the behavior of S αβ (∆G) and I αβ (∆G) under the C 3 rotation.
Using the Fourier relation c τ α (G) = 1 Ω´Ω u τ α (r)e −iG·r dr in which Ω is the area of the 2D cell, S αβ (∆G) can be rewritten as where g ≡ ∆G − 2K. Equation (12) shows that its integrand is a cell periodic function. This periodicity leads to the invariance of S αβ under the C 3 rotation of its integrand. For convenience, we rewrite S αβ (g) ≡ S αβ (∆G). Using the form of integrand in equation (13), we get where γ τ α ≡ [C 3 ϕ τ α (r)]/ϕ τ α (r) is the eigenvalue of C 3 operation corresponding to the Bloch function ϕ τ α (r). For the integral I αβ , we also rewrite it as I αβ (g) ≡ I αβ (∆G − 2K). Using the integral form of the Bessel function of order n (integer) J n (x) = 1 2π´π −π e i(x sin φ−nφ) dφ, we can simplify I αβ (g) to single integrals as where φ is the polar angle of the 2D vector g. It is obvious that under g → C 3 g, the integral parts of equations (15)-(18) do not change, and hence the ratio η αβ ≡ I αβ (C 3 g)/I αβ (g) is determined only by the prefactor in front of the integrals: Because the set of g's, G, also has the C 3 symmetry, we can choose one third elements of G, denoted as G 1 3 , such that G = 2 n=0 C n 3 G 1 3 . Using equations (14) and (19), we have where η αβ is already given in equation (19) and the rest is to determine γ τ α . For MoS 2 monolayer, the conduction and valence band edge Bloch functions at ±K points are formed predominantly by the d 0 (i.e. d z 2 ) and d ±2 [i.e. 1 √ 2 (d x 2 −y 2 ± id xy )] orbitals of Mo atoms respectively [? ? ? ]. There are two contributions to the eigenvalue γ τ α of the C 3 rotation on the Bloch function ϕ τ α (r) = e iτ K·r u τ α (r) (details in Appendix A): (i) the rotation of atomic orbital around its own center, i.e. C 3 d 0 (r) = d 0 (r) and C 3 d ±2 (r) = e ±i 2π 3 d ±2 (r); (ii) the change of lattice phase, defined as the value of the planewave phase factor e iτ K·r at each lattice site (cf. figure 1 for the case of τ = +1). The change of the lattice phase factor under C 3 rotation depends on the rotation center. If the rotation center is at Mo atom as shown in figure 1(a), the C 3 rotation does not change the lattice phase factors leaving only contribution (i) taking effect, and we obtain γ τ c = 1 and γ τ v = e iτ 2π 3 (for Mo center).
However, if the rotation center is at S atom, as shown in figure 1(b), both contributions (i) and (ii) take effect and we have γ τ c = e iτ 2π 3 and γ τ v = e −iτ 2π 3 (for S center).
Putting equations (19) and (22) [or (23)] into equation (21), we obtain that λ αβ = 1 for Mo centered potential and λ αβ = e −i 2π 3 for S centered potential, and the final intervalley coupling is 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 Mo atom, while it vanishes if the potential is centered at S atom. This is due to the dependence of γ τ α , the eigenvalue of the Bloch function ϕ τ α under the C 3 rotation, on the location of rotation center. Below, we show that the same conclusions can be drawn as long as the confinement potential has the C 3 symmetry, i.e. intervalley coupling is strongest (zero) if the center of the potential is at Mo (S) site.
A. Noncircular QD with C3 symmetry Consider a noncircular QD confinement potential with C 3 symmetry only, C 3 U (r) = U (r). In such potential that lacks the circular symmetry, we do not have exact solution in general for the massive Dirac Fermion Hamiltonian in equation (3). Nevertheless, the numerical results presented in section II B have well justified 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 ψ τ c can be obtained by solving the following effective-mass Hamiltonian where m eff is the effective mass of conduction band at ±K. With the same m eff at ±K, the envelope functions in the two valleys are the same: ψ + c = ψ − c . To analyze the symmetry, we need to know γ τ c and η cc [see equations (20) and (21)]. γ τ c is independent of the shape of the QD, and equations (22) and (23) are still valid here. η cc is determined by the symmetry of the envelope function ψ τ c . Since U (r) is symmetric under C 3 rotation, the non-degenerate ground-state envelope function ψ τ c in each valley has to be an eigenstate of C 3 , i.e. C 3 ψ τ c = ρψ τ c where |ρ| = 1. Using equation (10) and ρ * ρ = 1, we have Thus η cc = 1, consistent with equation (19). The conclusion of equation (24) is therefore valid for QD confinement potential with C 3 symmetry. In summary, the band edge Bloch functions at K and −K have the C 3 rotational symmetry of the lattice. However, they transform differently when the rotation center is Mo site or S site, as shown in Figure 1a and 1b respectively. When the QD confinement potential also has the C 3 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.

IV. REAL-SPACE TIGHT-BINDING METHOD
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 real-space tight-bind (RSTB) method. Using RSTB, we first calculate intervalley couplings of circular shaped QDs and compare with the results from EFM. The agreement between the two methods well justifies 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.

A. The RSTB model
We use supercell and periodic boundary condition to model a QD. The supercell is a square with side length L sc (cf. figure  9). Each supercell has N = L sc * round(2L sc / √ 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 [? ], 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 r n (n = 1, 2, · · · , N ). The tight-binding Hamiltonian of the QD is where U (r i ) is the QD confinement potential at r i , α, β = d z 2 , d xy , d x 2 −y 2 are orbital indices, c † iα is the creation operator for electron on site i with orbital α. t iα,jβ is the hopping between α orbital at position i and β orbital at position j, and , denotes summation over nearest-neighbor pairs only. In practice, we write the Hamiltonian H in form of a 3N × 3N matrix. The part of H involving hoppings between site i and j can be written as a 3 × 3 block H ij = h(r j − r i ), where r j − r i is either 0 or the six nearest-neighbour vectors R 1 to R 6 defined in [? ]. Using the real hopping parameters { 1 , 2 , t 0 , t 1 , t 2 , t 11 , t 12 , t 22 } defined in [? ] (the GGA-version parameters therein), we have, The 3N × 3N RSTB Hamiltonian matrix for a QD is then diagonalized to find the eigenstates and eigenenergies. Unlike the EFM approach where the bound state wavefunctions are first given in each valley, and the intervalley coupling then calculated as the off-diagonal matrix element between them [cf. 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 inter |. Thus, the magnitude of the intervalley coupling V inter is read out from the difference between the lowest two conduction-band energy eigenvalues of H in the RSTB solution [? ].

B. Circular shaped QDs
We study here 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 inter | vs. L sc relation for different circular QDs. We can see that |V inter | converges fast with increasing L sc , which means that the finite-size effects introduced by the method can be well eliminated by having large enough supercell. Figure 6 also shows that |V inter | converges faster with increasing L sc for larger potential depth V , just as expected.
|V inter | calculated using RSTB method with a large supercell (L sc = 240) is compared with the results from EFM in figure  7. For the |V inter | vs. 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.
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 inter | vs. R curves of the same circular shaped potential centered at a Mo or S site. Consistent with the symmetry analysis in Section III, |V inter | is three orders of magnitude larger for the Mo-centered case compared to the S-centered ate, the latter can be regarded as the numerical errors from zero. Figure 8(b) plots |V inter | as a function of the central position of the confinement potential. Clearly, |V inter | decreases when the potential center moves away from the Mo site and eventually vanishes when the potential center is at the S site.
C. 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 [? ]. 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 [? ? ? ? ? ? ? ? ]. Smaller QDs formed by the lateral heterostructures can be well 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 RSTB method to investigate intervalley coupling of these noncircular QDs. We studied QDs with triangular, hexagonal as well as the square shapes for comparison.
As shown in figure 9, the triangular, hexagonal, and square QDs are defined by the confinement potentials U T (r), U H (r), and U S (r) respectively, where U T/H/S equals to 0 within the QDs (yellow regions), and equals to 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 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 (cf. figure 5). We note that in heterostructures formed between different monolayer dichalcogenides (e.g. MoSe 2 and WSe 2 ), the conduction and valence band edge discontinuities are discovered to be in the range of 0.2-0.4 eV [? ]. For this range of values for V , the intervalley coupling is well negligible for the triangular and hexagonal QDs.
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 examine the change of the intervalley coupling when we rotate the QD confinement potential by an angle θ 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 θ rot significantly, and reaches a maximum at θ rot = 30 • where the intervalley coupling becomes comparable with the square QDs. Note that at θ rot = 30 • , the sides of the confinement potential of the triangular and hexagonal QDs are along the armchair crystalline axes [cf. figure 12(a) and (b)]. Figure 13 shows the dependence of |V 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 inter | strongly depends on the potential center. For triangular and hexagonal QDs with the C 3 rotational symmetry, the dependence is similar to the circular QDs [cf. 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 lacks the C 3 rotational symmetry.  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 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 (cf. figure 4), i.e. |V inter | drops fast by several orders of magnitude when L increases from 0 to 5, and becomes largely independent of L for L > 0. From figure 14, we can see that the |V 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.

V. CONCLUSIONS AND DISCUSSIONS
To conclude, we have investigated intervalley coupling in QDs defined by confinement potentials of various shapes on extended TMD monolayer. The numerical results obtained using two completely different approaches agree well with each other. For confinement potentials with the C 3 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 [? ] and in silicon QDs [? ? ? ? ? ? ? ? ]. 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 [cf. 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 [? ]. 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 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 [? ? ? ? ? ? ? ? ]. 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 well inherit the valley physics of the 2D bulk such as the valley optical selection rules, making the valley pseudospin of single confined electron a perfect candidate of 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.