Spin-orbit coupling, quantum dots, and qubits in monolayer transition metal dichalcogenides

We derive an effective Hamiltonian which describes the dynamics of electrons in the conduction band of transition metal dichalcogenides (TMDC) in the presence of perpendicular electric and magnetic fields. We discuss in detail both the intrinsic and the Bychkov-Rashba spin-orbit coupling (SOC) induced by an external electric field. We point out interesting differences in the spin-split conduction band between different TMDC compounds. An important consequence of the strong intrinsic SOC is an effective out-of-plane $g$-factor for the electrons which differs from the free-electron g-factor $g\simeq 2$. We identify a new term in the Hamiltonian of the Bychkov-Rashba SOC which does not exist in III-V semiconductors. Using first-principles calculations, we give estimates of the various parameters appearing in the theory. Finally, we consider quantum dots (QDs) formed in TMDC materials and derive an effective Hamiltonian which allows us to calculate the magnetic field dependence of the bound states in the QDs. We find that all states are both valley and spin split, which suggests that these QDs could be used as valley-spin filters. We explore the possibility of using spin and valley states in TMDCs as quantum bits, and conclude that, due to the relatively strong intrinsic spin-orbit splitting in the conduction band, the most realistic option appears to be a combined spin-valley (Kramers) qubit at low magnetic fields.


I. INTRODUCTION
Monolayers of transition metal dichalcogenides [1] (TMDCs) posses a number of remarkable electrical and optical properties, which makes them an attractive research platform. Their material composition can be described by the formula MX 2 , where M = Mo or W and X = S or Se. They are atomically thin, two-dimensional materials, and in contrast to graphene [2], they have a finite direct optical band gap of ≈ 1.5 − 2 eV, which is in the visible frequency range [3,4]. This has facilitated the theoretical [5] and experimental [6][7][8][9][10][11] study of the rich physics related to the coupling of the spin and the valley degrees of freedom.
Very recently, there has also been a growing interest in the transport properties of these materials. Although contacting and gating monolayer TMDCs is not entirely straightforward experimentally, progress is being made in this respect [12][13][14][15][16][17][18]. Electric [17] and magnetic field [19,20] effects are also studied currently, both in monolayer and few-layer samples. In addition, a promising experimental work has recently appeared regarding spin-physics in these materials, showing, e.g., a viable method for spin-injection from ferromagnetic contacts [16].
The finite band gap in the TMDCs should also make it possible to confine the charge carriers with external gates and therefore to create, e.g., quantum dots. Together with the above mentioned progress in contacting and gating TMDCs, this raises the exciting question of * e-mail: andor.kormanyos@uni-konstanz.de whether these materials could be suitable platforms to host qubits [21]. Our work is motivated by this question. First, we are going to introduce an effective Hamiltonian which accurately describes the physics in the conduction band (CB) of TMDCs in the (degenerate) K and K ′ valleys of the Brillouin zone (BZ). We confine our attention to the CB while the effect of the valence band (VB) and other relevant bands are taken into account through an appropriate choice of the parameters appearing in the model. This approach is motivated by the facts that i) the band-gap energy E bg is large with respect to other energy scales appearing in the problem, and ii) according to experimental observations, the samples of TMDCs are often intrinsically n-doped [16,22] or show unipolar n-type behavior [23]. To obtain realistic values of the parameters appearing in the theory we have performed density functional theory (DFT) calculations. We discuss the important effects of the intrinsic SOC which manifest themselves both through the spin-splitting of the CB and the different effective masses associated with the spin-split bands. We also point out that a perpendicular magnetic field, in addition to the usual orbital effect, leads to the breaking of valley degeneracy. Moreover, due to the strong SOC, the coupling of the spin degree of freedom to the magnetic field is described by an out-of-plane effective g-factorg ⊥ sp . We then study the effect of an external electric field and derive the Bychkov-Rashba SOC Hamiltonian for TMDCs. This is motivated by recent experiments [11,22], where strong electric fields were created by back gates to study the charged excitons. In particular, we find that in contrast to III-V semiconductors and graphene, due to the lower symmetry of the system, the Bychkov-Rashba SOC Hamiltonian contains two terms, one of which has not yet been discussed in the literature. Using perturbation theory and first-principles (FP) calculations, we can estimate the magnitude of this effect for each TMDC material.
Finally, we consider QDs obtained by confining the charge carriers with gate electrodes. We study the dependence of the spectrum of such QDs on a perpendicularly applied external magnetic field. We show that while pure spin and pure valley qubits are possible, e.g., in small QDs in MoS 2 , but they require large magnetic fields because of the relatively strong intrinsic SOC in the CB. On the other hand, combined spin-valley qubits represented by a Kramers pair can be operated at small magnetic fields. QDs in nanowires consisting of a MoS 2 nanoribbon with armchair edges or crystallographically aligned confining gates have been recently discussed [24]. Our proposal does not require atomically sharp boundaries or a precise control of the placement of the confining gates; therefore it should be easier to fabricate experimentally. Moreover, we explicitly take into account the intrinsic spin-splitting of the CB.
The paper is organized as follows. In Sec. II we derive an effective Hamiltonian describing electrons in the CB. We take into account the effects of perpendicular external electric and magnetic fields. Using the results of FP calculations we obtain values for the important parameters appearing in our model. In Sec. III we use this model to study the magnetic field dependence of the bound states in a QD. We also discuss the possible types of qubits that QDs in TMDCs can host. We conclude in Sec. IV. In Appendices A and B we present the details of the derivation of the effective Hamiltonian. We collect some useful formulas in Appendix C and the details of our DFT calculations can be found in Appendix D.

II. EFFECTIVE HAMILTONIAN
We consider a monolayer TMDC and introduce a lowenergy effective Hamiltonian which captures the most important effects in the spin-split conduction band at the K (K ′ ) point. The detailed derivation of the model, which is based on a seven-band (without the spin degree of freedom) k · p Hamiltonian, is presented in Appendix A. It is important to note that, as pointed out in Refs. 25-27, there are several band extrema in the band structure of TMDCs which can be of importance: see Fig. 2, where we show the band structure of MoS 2 obtained from DFT calculations. Since we assume that the system is n-doped, the maximum at the Γ point of the VB is not relevant. More important are the secondary minima in the CB, which are usually called the Q (a.k.a. T ) points. The exact alignment of the Q point energy minimum with respect to the K point minimum is difficult to deduce from DFT and GW calculations, because it depends quite sensitively on the details of these computations [28]. We have found that using the local density approximation (LDA), all compounds, with the exception of MoS 2 , become indirect gap semiconductors if we take into account the SOC, because the Q point minimum is lower than the K point minimum. More advanced GW calculations also give somewhat conflicting results and are quite sensitive to the level of theory [29] (G 0 W 0 , GW 0 etc,) and the lattice constant used. Experimentally, monolayer TMDCs show a significant increase of photoluminescence [10,22,30,31] with respect to few-layer or bulk TMDCs, which is usually interpreted as evidence that they are direct gap semiconductors. Therefore we assume that for low densities it is enough to consider only the K and K ′ points of the CB. For the formation of QDs from states around the K point, the safest material appears to be MoS 2 , where the secondary minima are most likely above the K point minimum by a few hundred meV [26,32]. However, for operation at low temperatures, the other TMDCs may also be suitable, as long as the Q point lies a few meV higher than the K points. In cases where the Q point lies below the K point, one can envisage QDs formed within the Q valley, but this is beyond the scope of this paper.
A. Electronic part and intrinsic spin-orbit coupling Due to the absence of a center of inversion and strong SOC, the bands of monolayer TMDC materials are spinsplit everywhere in the Brillouin zone (BZ), except at the high-symmetry points Γ and M , where the bands remain degenerate. In addition, the projection of the spin onto the quantization axis perpendicular to the plane of the monolayer is also preserved. This is a consequence of another symmetry, namely, the presence of a horizontal mirror plane σ h . Therefore, a suitable basis to describe the CB is given by the eigenstates ↑, ↓ of the dimensionless spin Pauli matrix s z with eigenvalues s = ±1. In what follows, we will often use the shorthand notation ↑ for s = 1 and ↓ for s = −1.
In the absence of external magnetic and electric fields, the effective low-energy Hamiltonian which describes the spin-split CB at the K (K ′ ) point in the basis ↑, ↓ is Here, we introduce the inverse effective mass 1 for K (K ′ ) and the wavenumbers q ± = q x ± iq y are measured from the K (K ′ ) point. Leaving the discussion of the effects of magnetic field to Sec. II B, we set q + q − = q 2 x + q 2 y and therefore the dispersion described by the Hamiltonian (1) is parabolic and isotropic. The trigonal warping [26], which is much more pronounced in the VB than in the CB, is neglected here.
The strong spin-orbit coupling in TMDCs has two consequences: firstly, as already mentioned, the CB is spinsplit at the K (K ′ ) point and this is described by the parameter ∆ cb . Secondly, the effective mass is different for the ↑ and ↓ bands. Our sign convention for the effective mass assumes that the spin-up band is heavier than the spin-down band at the K point (for details on the effective mass calculations see Appendix B). The effective mass m K,s eff of different TMDCs, obtained from fitting the DFT band structure [33], is shown in Table I (note that m K ′ ,s eff = m K,−s eff ). As one can see, the difference between m K,↑ eff and m K,↓ eff is around 10 − 14% for MoS 2 and MoSe 2 , while it is 30% for the WX 2 compounds. In the sevenband k·p model this can be explained by the fact that the effective mass depends on the ratio of the spin splittings in other bands (most importantly, in the VB and the second band above the CB) and the band gap E bg . For the heavier compounds the spin-splittings are larger, but E bg remains roughly the same or even decreases, leading to a larger difference in the effective masses.
The results of DFT calculations also suggest that in the case of MoX 2 materials there are band crossings between the spin-split CB because the heavier band has higher energy. For WX 2 materials such a band crossing is absent. Taking MoS 2 and WS 2 as an example, the dispersion in the vicinity of the K point is shown in Fig. 3.  We note that a model Hamiltonian similar to Eq. (1), but without taking into account the difference in the effective masses, has been used in Refs. 34 and 35 to study spin-relaxation processes in MoS 2 . The effective mass difference and the sign of the effective SOC in the CB has also been discussed recently in Ref. 36.

B. Effects of a perpendicular magnetic field
We assume that a homogeneous, perpendicular magnetic field of strength B z is applied. The k·p Hamiltonian can be obtained by using the Kohn-Luttinger prescription, which amounts to replacing the numbers q x and q y in the above formulas with operators: q →q = 1 i ∇+ e A, where A is the vector potential in Landau gauge and e > 0 is the magnitude of the electron charge. Note that due to this replacementq + andq − become noncommuting operators: [q − ,q + ] = 2eBz , where |B z | is the strength of the magnetic field. Therefore their order has to be preserved when one folds down a multi-band Hamiltonian, which lies behind the low-energy effective Hamiltonian (1). As a consequence, for finite magnetic field further terms appear in the effective Hamiltonian. The derivation of these terms within a seven-band k · p model is given in Appendix B.
One finds that in an external magnetic field H τ,s el in Eq. (1) is replaced bỹ where ω τ,s c = e|B z |/m τ,s eff . The term ∼ ω τ,s c in the bulk case introduces a shift in the index of the Landau levels, so that there is an "unpaired" lowest Landau level in one of the valleys. The next term,H τ vl = −τg vl µ B B z , breaks the valley symmetry of Landau levels. Hereg vl is the "valley gfactor". Similar effects have also been found in gapped monolayer [37] and bilayer [38,39] graphene, and has recently been noted for MoS 2 as well [40][41][42]; therefore we do not discuss them here in detail.
A new term, to our knowledge not yet considered in the literature of monolayer TMDC, is due to the strong SOC in these materials. It can be written in terms of an outof-plane effective spin g-factor g ⊥ so :H s sp = 1 2 g ⊥ so µ B s z B z where µ B is the Bohr magneton. In addition, the wellknown Zeeman term H Z = 1 2 g e µ B s z B z also has to be taken into account [43]. Here g e ≈ 2 is the free-electron g-factor. The coupling of the spin to the magnetic field can therefore be described bỹ where the total g-factor in the CB isg ⊥ sp = g e + g ⊥ so . Values ofg vl and |g ⊥ so | obtained with the help of our DFT calculations are shown Table II. The sign of g ⊥ so cannot be obtained with our methods; it should be deduced either from experiments or from more advanced FP calculations. For the numerical calculations in Sec. III A we will assume that g ⊥ so > 0. In Sec. III A we will study the interplay of the magnetic field and the quantization due to confinement in QDs. While Eq. (4) is a convenient starting point to understand the Landau level physics, for relatively weak magnetic fields, when the effect of the confinement potential is important with respect to orbital effects due to the magnetic field, one may re-writeH τ,s el ,H τ vl , and H s sp,tot in a slightly different form: where g vl = (2m e /m 0 eff ) −g vl and g ⊥ sp =g ⊥ sp − (2m e /δm eff ). This form shows explicitly that in contrast to H τ,s el , which depends on the product of τ and s (through m τ,s eff ), H τ vl and H s sp,tot depend only on τ and s z , respectively. This can help to understand the level splittings patterns in QDs: see Sec. III A. In particular, for states which form a Kramers pair τ · s = 1 or −1, therefore H τ,s el , which only depends on the product of τ and s, would not lift their degeneracy in the presence of a magnetic field. Due toH τ vl , however, the degeneracy of the Kramers pair states will be lifted. Assuming g ⊥ so > 0 and B z > 0, as in the calculations that lead to Figs. 4 and 5, the values of g vl and g ⊥ sp are shown in Table II.

C. External electric field and the Bychkov-Rashba SOC
The effective Hamiltonian (1), describing the dispersion and the spin splitting of the CB is diagonal in spin space. An external electric field has two effects: i) it can induce Bychkov-Rashba type SOC which will couple the different spin states, and ii) it can change the energy of the band edge. We start with the discussion of the Bychkov-Rashba SOC.
For simplicity, we assume that the external electric field is homogeneous and that its strength is given by E z . Then the Bychkov-Rashba SOC in TMDCs is described by the Hamiltoniañ The first term, λ i BR (s y q x − s x q y ), is the well-known Bychkov-Rashba [44] Hamiltonian, which is also present in GaAs and other III-V semiconductor compounds. It is equivalent to the Bychkov-Rashba Hamiltonian recently discussed in Ref. 45 in the framework of an effective twoband model, which includes the VB. The second term, λ r BR (s x q x + s y q y ), is also allowed by symmetry (see Table I of Ref. 46) because the pertinent symmetry group at the K point in the presence of an external electric field is C 3 . A derivation of the Hamiltonian (5) is given in Appendices A and B. We note that the coupling constants λ r BR and λ i BR cannot be tuned independently, because both of them are proportional to the electric field but with different proportionality factors. Using our microscopic model and FP calculations similar to those in Ref. 47, we can estimate the magnitude of λ BR but not λ r BR and λ i BR separately. The |λ BR | values that we have obtained are shown in Table III. They give an upper limit for the real values because we have neglected, e.g., screening in these calculations (for details see Appendix B). More advanced DFT calculations, such as those recently done for bilayer graphene [48], would be certainly of interest here.  Table III to the values found in InAs [49] or InSb [50], one can see that for relatively small values of the electric field (E z 10 −2 V/Å), where the perturbation theory approach can be expected to work, |λ BR | is smaller by an order of magnitude than in these semiconductor quantum wells. Nevertheless, the Bychkov-Rashba SOC is important because it constitutes an intra-valley spin-relaxation channel, which does not require the simultaneous flip of spin and valley. Thus, it may play a role in the quantitative understanding of the relaxation processes in the recent experiment of Jones et al. [11], where a large back gate voltage was used.

Comparing the numbers shown in
The external electric field has a further effect, which, however, turns out to be less important for our purposes. Namely, it shifts up the band edge of the CB, and the shift is, in principle, spin dependent [see Eqs. (B2c), (B3c) in Appendix B]. The shift of the CB edge can be understood in terms of the electric field dependence of the band gap (we note that the band edge of the VB also depends on the electric field, and the shifts of the VB and CB edges together would describe the change of the band gap). In contrast to Ref. 40, however, in our model the shift of the band edge depends quadratically on the strength of the electric field and not linearly. We think this is due to the fact that in the model used in Ref. 40 the p orbitals of the sulfur atoms are admixed only to the CB. In fact, symmetry considerations [26,45] and our DFT calculations show that the p (or d) orbitals of the X atoms have a small weight at the K point both in the VB and in the CB . Taking this into account, as in the tight-binding model of Ref. 27, one would find that for weak electric field regime the dependence of the band gap is quadratic in the electric field. Moreover, both our perturbation theory and preliminary DFT results suggest that the shift of the band edge in the CB is actually very small, at least in the regime where the perturbation theory approach is applicable (see Appendix B for details). Therefore we neglect it in the rest of the paper. The spindependence of the band-edge shift, being a higher-order effect, is expected to be even smaller.

A. Quantum dots in TMDCs
QDs in novel low-dimensional structures, such as bilayer graphene [38,[51][52][53] and semiconductor nanowires with strong SOC [54,55], are actively studied and the applicability of these structures for hosting qubits has also been discussed. Motivated by the interesting physics revealed in these studies, we now consider QDs in twodimensional semiconducting TMDCs defined by external electrostatic gates. In particular, we will be interested in the magnetic field dependence of the spectrum and discuss which eigenstates can be used as two-level systems for qubits. We consider relatively small QDs which can be treated in the ballistic limit. The opposite limit, where disorder effects become important and the spectrum acquires certain universal characteristics, can be treated along the lines of Ref. 56, but this is beyond the scope of the present work.
Nevertheless, based on the findings of Sec. II A, the following general considerations can be made: assuming a chaotic QD with mean level spacing δ = 2π 2 /(m eff A), where A is the area of the dot, one can see that one needs relatively small QD in order to make δ larger than the thermal energy k B T . For instance, taking a dot of radius R ≈ 40 nm we find for, e.g., MoS 2 that δ ≈ 0.2 meV, corresponding to T = 2.3 K, whereas for WS 2 , due to its smaller effective mass, the mean level spacing is T ≈ 3.4 K. In this respect TMDCs with smaller m eff , such as WS 2 and WSe 2 , might be more advantageous. Although the required temperatures are smaller than in the case of GaAs (which has m eff ≈ 0.067m e ), they are still achievable with present-day techniques.
In the following, for simplicity, we will study circular QDs because their spectrum can be obtained relatively easily and can illustrate some important features of the spectrum of more general cases. In particular, we will consider QDs in MoS 2 and WS 2 . The total Hamiltonian in the K, K ′ valleys (τ = ±1) reads where V dot is the confinement potential for the QD. As we have shown,H τ BR is relatively small; therefore we treat it as a perturbation, whereas the stronger intrinsic SOI is treated exactly. The Hamiltonian of the nonperturbed system is given by i.e., it is diagonal both in valley and in spin space. We consider a circular QD with hard wall boundary condi- In cylindrical coordinates, the perpendicular magnetic field can be taken into account using the axial gauge, where A φ = B z r/2 and A r = 0. With this choice, since the rotational symmetry around the z axis is preserved, H dot commutes with the angular momentum operatorl z and they have common eigenfunctions. , where l B = eBz is the magnetic length, one finds for B z > 0 that The eigenfunctions of the operatorsα + andα − , which are (i) regular at ρ = 0 and (ii) also eigenfunctions ofl z , are g a,l (ρ, ϕ) = e ilϕ ρ |l| 2 e − ρ 2 M (a, |l| + 1, ρ), where l is an integer and M (a, |l|+1, ρ) is the confluent hypergeometric function of the first kind [58]. One can show that (For details see Appendix C.) Considering now the Schrödinger equation for the bulk problem, i.e., for V dot = 0 in valley τ for spin s, it reads [ ω τ,s cα+α− + where Θ(x) is the Heaviside step function. The wave |l| + 1, ρ) and   Tables I and II.
Here E τ,s = 1/2sgn(B z ) ω τ,s c + τ s ∆ cb + 1 2 (τ g vl µ vl + s g ⊥ sp µ B )B z − E. The bound state solutions of the QD problem are determined by the condition that the wave function has to vanish at r = R d , i.e., one has to find the energy E τ,s l for which M (a l , |l| + 1, ρ[r = R d ]) = 0. The task is therefore to find for a given magnetic field B z and quantum number l the roots of M (a l , |l|+ 1, ρ[r = R d ]) = 0 as a function of a l . The a l values can be calculated numerically. Once the nth root a n,l is known, the energy of the bound state E τ,s n,l can be expressed using Eq. (11). The numerically calculated spectrum for a QD with R d = 40 nm in MoS 2 is shown in Fig. 4(a). At zero magnetic field, because of the quadratic dispersion in our model, there is an effective time reversal symmetry acting within each valley and therefore states with angular momentum ±l within the same valley are degenerate. For finite magnetic field all levels are both valley and spin split. For even larger magnetic fields, when l B R d , the dot levels merge into Landau levels. Since ∆ cb is relatively small with respect to the cyclotron energy ω τ,s c , spin-split states ↓ and ↑ from the same valley can cross at some larger, but still finite magnetic field (see, e.g., the crossing between the black and green lines for E > 3 meV for states in valley K in Fig. 4a).  Table I, whereas g vl = 1.6, and g ⊥ sp = 1.99 (see Table II).
Taking into account the Bychkov-Rashba SOC turns the crossings between states |a, l, ↑ and |a, l + 1, ↓ , l ≥ 0 into avoided crossings. The selection rules for H τ BR can be derived by rewritingH τ BR in terms of the operators α − and α + and calculating their effect on the non-perturbed eigenstates (see Appendix C for details). For the lowlying energy states, in which we are primarily interested, the effect of the Bychkov-Rashba SOC is to introduce level repulsion between these states and higher energy ones allowed by the selection rules. Taking |λ BR |/l B as a characteristic energy scale of this coupling and using Table III one can see that for magnetic fields 10 T and electric fields E z 10 −2 V/Å the level repulsion is much smaller than the spin splitting ∆ cb and therefore we neglect it. Figure 4(b) shows the low-field and low-energy regime of Fig. 4(a). As one can see, for B z 1 T the lowest energy states reside in valley K. We emphasize that, in contrast to gapped monolayer [38,59,60] and bilayer [38,60] graphene, the energy states are also spin polarized. This suggest that QDs in MoS 2 can be used as simultaneous valley and spin filters. Figure 5 shows the low-energy spectrum of a WS 2 QD with radius R d = 40 nm. Qualitatively, it is similar to MoS 2 , but because the spin splitting ∆ cb between the ↑ and ↓ states belonging to the same valley is much larger than was the case for MoS 2 , they do not cross for the magnetic field range shown in Fig 5. One can also observe that the B z = 0 level spacing is somewhat larger than in the MoS 2 QD [see Fig. 4(b)]. Another important observation that can be made by comparing the results for MoS 2 and WS 2 is the following: for a given magnetic field, e.g., B z = 5 T, the splitting between states belonging to different valleys is significantly larger for the former material than for the latter (compare Figs. 4(b) and 5). This is due to the different sign of ∆ cb and hence different spin polarization of the lowest levels in the two materials: in the case of MoS 2 the valley split-  Table I and we used g vl = 2.31 and g ⊥ sp = 1.84 (c.f. Fig. 5). Black (red) lines show spin ↑ (↓) states from the K (K ′ ) valley. ting (described by H τ vl ) and the coupling of the spin to the magnetic field (given by H sp,tot ) reinforce each other, whereas for WS 2 they counteract, and since g vl and g ⊥ sp have similar magnitude, in the end the valley splitting of the levels at large magnetic fields is small. This suggests that for spin and valley filtering the MoX 2 compounds are better suited.
The qualitative difference between MoS 2 and WS 2 regarding the valley splitting does not depend crucially on the exact values of the bulk parametersg vl and g ⊥ so . However, on a more quantitative level, the valley splitting does depend on the exact values of the valley and spin gfactors, which were calculated using the DFT band gap and the k · p parameter γ 3 (see Sec. B for details). It is known that DFT underestimates the band gap, and the value of γ 3 depends to some extent on the way it is extracted from the FP computations. As a result, the values shown in Table II probably overestimateg vl and g ⊥ so . To illustrate this point, we show in Fig. 6 the lowenergy spectrum of the same WS 2 quantum dot as in Fig. 5 but using a g vl (g ⊥ sp ) which was obtained from ã g vl (g ⊥ so ) that is ∼ 20% smaller than the one shown in Table II. The valley splitting of the bound states can now barely be observed.

B. Qubits in TMDC quantum dots
Circular hard-wall QDs in two-dimensional semiconducting TMDCs have a spectrum similar to the characteristic Fock-Darwin spectrum for harmonically confined QDs (Fig. 4). Taking MoS 2 as an example, due to the intrinsic spin-orbit splitting of about 3 meV, each of the spin-and valley-degenerate states |l splits into two Kramers pairs at vanishing magnetic field B = 0, namely (|l, K, ↑ , |l, K ′ , ↓ ) and (|l, K ′ , ↑ , |l, K, ↓ ). Only at relatively high magnetic fields do we observe a crossing of two states with the same spin and opposite valley or within the same valley with opposite spin. These valley and spin pairs could serve as valley or spin qubits, respectively, but the required high magnetic field and the other overlapping levels with different l ′ quantum numbers complicate their realization. (The energy of higher angular momentum states can in principle be increased by making the QD smaller).
In view of the above, the most realistic approach seems to be to use the lowest Kramers pairs around B = 0, e.g., |l = 0, K ′ , ↑ and |l = 0, K, ↓ as a combined spin-valley qubit [54,61]. The energy splitting of these two-level systems could be tuned using the external magnetic field. The relaxation time of such spin-valley qubits in TMDC QDs will be limited only by the longer spin or valley relaxation time, while the pure dephasing time will be limited by the shorter of the two. The exchange interaction then provides the necessary coupling of adjacent spin-valley qubits for the realization of two-qubit gates.

IV. SUMMARY
In summary, we have studied TMDCs as possible host materials for QDs and qubits. We considered n-doped samples, which can be described by an effective model which involves only the CB. Using our FP calculations, we have obtained the parameters that appear in the effective Hamiltonian (effective masses, g-factors) for four distinct TMDC materials. We discussed the effects of external magnetic and electric fields, pointing out that the former leads to the splitting of the energy levels in different valleys, while the latter induces a Bychkov-Rashba SOC, which, however, appears to be rather small. We have used the effective Hamiltonian to calculate the spectrum of circular QDs, finding that all bound states are both spin and valley split. Our results suggest that at large magnetic field QDs in TMDCs can be used as spin and valley filters, but that this effect may depend on material-specific details. Finally, we have discussed the possible types of qubits that QDs in TMDC materials can host. We have found that Kramers pairs around B z = 0 appear to be the most realistic candidates.
The effective one-band model and the material parameters that we obtained for different TMDCs will hopefully be helpful in other fields as well, e.g., for studying plasmonic excitations [62].
Note added After the submission of this work another manuscript appeared on the arXiv and has been subsequently published [63] on the spin-splitting in the conduction band of monolayer TMDCs.

V. ACKNOWLEDGMENTS
We acknowledge discussions with Lin Wang. A. K. and G. B. acknowledge funding from DFG under programs SFB767, SPP1285, FOR912 and from the European Union through Marie Curie ITN S 3 NANO. V. Z. acknowledges support from the Marie Curie project CARBOTRON.
Appendix A: Seven-band model

Introduction
Our aim is to derive a low-energy effective Hamiltonian valid close to the K (K ′ ) point of the BZ, which describes the band dispersion, the effects of intrinsic SOC, and the SOC induced by an external electric field (Bychkov-Rashba effect). To this end we will consider the SOC in the atomic approximation, apply k · p perturbation theory, and take into account the effect of an external electric field perturbatively. We consider a seven-band model (without spin) which contains every band from the third band below the VB (which we call VB-3) up to the second band above the CB (denoted by CB+2 henceforth), i.e., we take the basis vb − 2, vb − 1, vb, cb, cb + 1, cb + 2} denotes the band and the lower index µ indicates the pertinent irreducible representation of the point group C 3h , which is the pertinent symmetry group for the unperturbed basis functions at the K point of the BZ. The spinful symmetry basis functions are represented by |Ψ b µ , s = |Ψ b µ ⊗ |s , where s = {↑, ↓} denotes the spin degree of freedom. Note, that the basis states can be separated into two groups. The first group contains those states whose orbital part is symmetric with respect to the mirror operation σ h : , s }; the second group contains antisymmetric states: A ′′ , s }.

Intrinsic spin-orbit coupling at the K (K ′ ) point of the Brillouin zone
The intrinsic SOC is treated in the atomic approximation, whereby the SOC is given by the Hamiltonian [43] Here V (r) is the spherically symmetric atomic potential,L is the angular momentum operator andŜ = (s x , s y , s z ) is a vector of spin Pauli matrices s x , s y , s z (with eigenvalues ±1). One can rewrite the productL·Ŝ asL·Ŝ =L z s z +L + s − +L − s + , whereL ± =L x ±iL y and s ± = 1 2 (s x ± is y ). The task is then to calculate the matrix elements of (A1) in the basis introduced in Sec. A 1 at the K ( K ′ ) point of the BZ. To this end one can make use of the symmetries of the band-edge wave functions. For instance, the diagonal matrix elements are proportional to s z , this is because theL z is symmetric with respect to σ h whereasL ± is antisymmetric. Conversely, most of the off-diagonal matrix elements will be proportional to s ± , reflecting the fact that they are related to matrix elements having different symmetry with respect to σ h . The only exception is the off-diagonal matrix element between |Ψ v−3 E ′ 2 , s and |Ψ c+2 E ′ 1 , s , which connects symmetric states. In addition, one has to consider the transformation properties of the basis functions and angular momentum operators with respect to a rotation by 2π/3. The general result for the K point is shown in Table IV.
Before showing further details of the calculations in subsections A 3 and A 4, some comments are in order here. As long as one considers states close to the K point, the largest energy scale is the band gap and other band-edge energy differences. The next largest energy scale comes from the SOC. As an upper limit of the various diagonal and off-diagonal matrix elements (see Table  IV) one can take the spin-splitting of the VB. The reason is that the main contribution to this band at the K point comes from the metal d orbitals and the metal atoms, being much heavier than the chalcogenides, are expected to dominate the SOC (with the possible exception of the CB). This is smaller than the typical interband energies for the MoX 2 materials and therefore the different bands are only weakly hybridized by the SOC. For the heavier WX 2 compounds the VB spin-splitting is 425−460 meV, indicating that some matrix elements may not be small any more with respect to band-edge energy differences. One is therefore tempted to perform first a diagonalization of the SOC Hamiltonian (see Table IV), to obtain the eigenstates |Ψ b µ,µ ′ , s which will be some linear combination of the original basis states |Ψ b µ , s , and then perform the k · p expansion and the perturbation calculation for the external electric field using this new basis. Diagonalization of the Hamiltonian (IV) is possible if one neglects the matrix elements ∆ v−3,c+1 , ∆ v−3,c+2 and ∆ v−2,c+2 between remote bands. The eigenstates are linear combinations of a symmetric and an antisymmetric basis vector. However, the subsequent calculations in Secs. A 3 and A 4 as well as the final Löwdin partitioning are more tractable if we do not make this diagonalization and stay with the original basis states throughout the calculations. The two approaches give the same results in the leading order of the ratio of the various SOC matrix elements and band-edge energy differences. For MoX 2 compounds the approach outlined below is adequate, for the heavier WX 2 materials it still gives reasonable results, but the numerical estimates for, e.g., the effective g-factor might have to be revised, once experimental and theoretical consensus is reached regarding the magnitude of the band gap and SOC band splittings. The SOC Hamiltonian at K ′ can be obtained by making the following substitutions: These relations follow from the fact the orbital wave functions at K and K ′ are connected by time-reversal symmetry, i.e., |Ψ b µ (K) =K 0 |Ψ b µ ′ (K ′ ) , whereK 0 denotes complex conjugation. Consider, as an Here we have made use ofK 0Lz = −L zK0 . Relations for the matrix elements involving the operatorsL ± can be obtained by noting thatK 0L± = −L ∓K0 and therefore 3. k · p matrix elements at the K (K ′ ) points The Hamiltonian H k·p = 1 2 me (q +p− + q −p+ ) has nonzero matrix elements only between states |Ψ b µ , s and |Ψ b ′ µ ′ , s which are either both symmetric or antisymmetric with respect to the mirror operation σ h . For the discussion in the main text we only need the matrix elements between symmetric states. These matrix elements, which are diagonal in the spin-space, have already been obtained in Ref. 26, but for convenience they are replicated in Table V. We note that in addition top ± , another operator due to SOC appears in the calculation of the k·p matrix elements [43,64], but it can be neglected. The diagonal elements in Table V are the band-edge energies. The matrix elements at the K ′ point can be obtained with the substitutions γ i → γ * i and q ± → −q ∓ . This follows from As mentioned in Ref. 26, concrete values for the γ i parameters can be obtained either from fitting the band dispersion or using the Kohn-Sham orbitals to evaluate directly the matrix elements Ψ b µ |p ± |Ψ b ′ µ ′ . The latter can be done, e.g., with the help of castep code (see Appendix D for computational details). To estimate the effective valley and spin g-factor (Sec. B 1) and the Bychkov-Rashba SOC parameter (Sec. B 4) we will need the value of γ 3 , for which the two approaches give similar results.
External magnetic field The effects of an external magnetic field in the k · p formalism can be obtained by using the Kohn-Luttinger prescription [43], which amounts to replacing the numbers q x , q y in the above formulas with the operatorŝ q = 1 i ∇ + e A, where A is the vector potential and e > 0 is the magnitude of the electron charge. Note that due to this replacementq + andq − become non-commuting operators and their order has to be preserved when one folds down the above multi-band Hamiltonian to obtain a low-energy effective Hamiltonian. Using the Landau gauge to describe a homogeneous, perpendicular magnetic field, the commutation relation is [q − ,q + ] = 2eBz .

External electric field
In order to derive the Bychkov-Rashba SOC, we assume that a homogeneous, perpendicular external electric field is present, which can be described by the Hamiltonian U (z) = eE z z. It breaks the mirror symmetry σ h and therefore couples symmetric and antisymmetric basis states, while the matrix elements between states of the same symmetry are zero. The full symmetry at the K point is lowered from C 3h to C 3 , i.e., the three-fold rotational symmetry is not broken. The matrix elements of H K U between the symmetric and antisymmetric states are shown in Table VI.
can be calculated using the band-edge Kohn-Sham orbitals, as in Ref. 47, where this approach was used to estimate the electric-field-induced band gap in silicene (see Appendix D for computational details). Since the Kohn-Sham orbitals are defined only up to an arbitrary phase, from the actual calculations we cannot extract the real and imaginary parts of ζ b,b ′ . The matrix elements at the K ′ point can be obtained by complex-conjugation of the K-point matrix elements.

Electronic effective Hamiltonian H el
In the electronic Hamiltonian H el we have taken into account the fact that in the presence of an external magnetic field the operatorsq + andq − do not commute. To obtain Eq. (1) in the manuscript, one has to use the commutation relation [q − ,q + ] = 2eBz and re-write 2q2 2me as 2q +q− 2me + eBz 2me . One finds in the K valley and in the K ′ valley. The effective mass m τ,s eff is given by 1 2m τ,s In the above formulasγ i = γ i / . The inverse of the effective mass m τ,s eff can be then re-written in terms of m 0 eff and δm eff , as shown below Eq. (1). The difference δm eff in the effective masses comes mainly from the spin-splitting ∆ v and ∆ c+2 of the VB and CB+2, respectively, other diagonal SOC matrix elements being much smaller. We attribute the heavier effective mass at the K point to the ↑ band. This assignment is based on the following. (i) From DFT calculations we know that both the VB and the CB+2 are composed mainly of d x 2 −y 2 and d xy orbitals. Using group theoretical considerations we take a VB Bloch wave function ∼ d x 2 −y 2 − id xy , whereas in the case of the CB+2 the Bloch wave function is ∼ d x 2 −y 2 + id xy . (ii) Taking into account (i) we as- (K) > 0. Regarding (i), we note that since the states at the K point are related to the states at K ′ by time reversal, our choice for the VB Bloch wave function is equivalent to other choices in the literature [5,27] up to a possible re-labeling of the valleys K ↔ K ′ . The sign of ∆ v , as shown below, affects the sign of the effective spin g-factor, therefore it should be possible to deduce it experimentally. (From symmetry considerations [25,26] and FP results [27] we also know that there is a small X-p orbital contribution to the VB and CB+2 as well, but in contrast to the CB, which is discussed in Sec. B 2, this can be neglected in the case of the VB and CB+2 spin-splitting.) The physical meaning of the term [2|γ 3 | 2 /(ε τ,s c − ε τ,s v )] eB z appearing in Eqs. (B4) and (B5) is probably more transparent if one expands it in powers of The higher-order terms in the expansion determine how the coupling of the spin to the magnetic field is modified due to the strong SOC in TMDCs. Keeping the first-order term only one arrives at the Hamiltoniañ H s sp = 1 2 g ⊥ so µ B B z where g so is an out-of-plane effective spin g factor, where m e is the bare electron mass. The value of ∆ c , i.e., the spin splitting coming from the X-p orbitals in the CB (see Sec. B 2) is not known; however, we can safely assume that it is negligible with respect to ∆ v . As explained above, we assume that ∆ v < 0, so we find that g ⊥ so ≈ 8m e |γ 3 | 2 |∆ v |/(E 2 bg ). We note that in the case of bulk semiconductors a similar formula to Eq. (B8) is called Roth's formula [65].
The relevant parameters ∆ v , |γ 3 |, and E bg to calculate g vl and g ⊥ so are shown in Table VII. in MoSe 2 (∆ MoSe2 c ≈ 23 meV) which contains a heavier chalcogenide than in MoS 2 (∆ MoS2 c ≈ 3meV). On the other hand, the above reasoning would suggest that because of the competition between the two SOC terms of different origins, the splitting in WS 2 (∆ WS2 c ≈ 38meV) should be larger than in WSe 2 (∆ W Se2 c ≈ 46meV), which is not the case according to our DFT calculations. This might be related to the larger orbital weight of the M-d orbitals in the relevant bands in the case of WSe 2 . In any case, the detailed understanding of the SOC in the CB requires further study.

Band-edge shift HU
The Hamiltonian H U in Eqs. (B2c) and (B3c), describes the dependence of the band edge on the external electric field. An order-of-magnitude estimate can be obtained by calculating ζ c,v−2 using LDA Kohn-Sham orbitals, generated by the castep code. As one can see from Table VIII, it is a small effect for the electric field values (E z 10 −2 V/Å), where the perturbation theory should be valid, and therefore we neglect it. We note, that as one can see in Eqs. (B2c) and (B3c), the value of H U also depends (indirectly) on E bg . The band gap, according to GW calculations [32,[68][69][70][71], is most likely to be underestimated by our DFT-LDA calculations. On the other hand, ξ c,v−2 is probably overestimated, because screening is neglected in our perturbative Kohn-Shamorbital-based calculations. As a consequence, the values shown in Table VIII  The shift of the band edge is, in principle, spindependent, but as one can see it from Eqs. (B2c) and (B3c), this is a higher-order effect and can be safely neglected.

Bychkov-Rashba Hamiltonian HBR
Finally, we discuss the Bychkov-Rashba Hamiltonian [ Eqs. (B2d) and (B3d)]. It is a sum of several terms, each having the same structure and related to the matrix elements ξ v,c+1 , ξ v−3,v−1 , ξ c+1,v−1 and ξ c,v−2 . Using Löwdin-partitioning, one finds for the most important term at the K point, To make the results more transparent, in the above formula we have neglected the spin-splittings of the CB and CB+1, which are much smaller than the splitting of the VB. The product γ * 3 ξ v,c+1 ∆ * c,c+1 is in general a complex number and therefore the Bychkov-Rashba coupling constant is also complex. By separating the real and imaginary part of λ One can estimate the magnitude of λ BR in the following way. As mentioned in Sec. A 4, one can calculate the magnitude of ζ z v,c+1 and the parameter γ 3 using the band-edge Kohn-Sham orbitals (see Table IX). The band-edge energies ε ↑,↓ c , ε ↓ v and ε ↑,↓ c+1 are known from DFT-LDA band structure calculations; we have collected their values in Table IX. Unfortunately, the off-diagonal SOC matrix element ∆ c,c+1 is not directly given by the DFT calculations. However, information about the weight of the M-d orbitals in each of the bands can be obtained from DFT computations and therefore we can relate this matrix element to ∆ v , because the dominant contribution to the SOC should come from the M-d orbitals. Since the M-d orbital weight in both the CB and the CB+1 band is similar to the one in the VB, we take |∆ c,c+1 | |∆ v |.
Similar procedure can be performed to estimate the terms proportional to the other non-zero ξ b,b ′ matrix elements as well. We have found that the magnitude of these further terms are significantly smaller than that of λ (1) BR , mainly because of the pre-factors which are inversely proportional to the product of band-edge energy differences between remote bands. Therefore, as an orderof-magnitude estimate of the strength of the Bychkov-Rashba SOC, one can just use λ (1) BR . Taking values for |γ 3 | from Table VII and for the other parameters from  Table IX, one finally arrives at the results shown in Table III.
The method outlined here most likely overestimates the real values of the Bychkov-Rashba parameters. In addition to the uncertainties in the values of the SOC matrix elements and the γ i parameters, there are two other sources of error: i) the calculation of ζ z b,b ′ did not take into account screening effects (see Ref. 47) and ii) according to GW calculations, the real band gap is larger then the DFT one, and this affects the energy denominators in the above formulas.