Effective theory of monolayer TMDC double quantum dots

Monolayer Transition Metal Dichalcogenides (TMDCs) are promising candidates for quantum technologies, such as quantum dots, because they are truly two-dimensional semiconductors with a direct band gap. In this work, we analyse theoretically the behaviour of a double quantum dot (DQD) system created in the conduction band of these materials, with two electrons in the (1,1) charge configuration. Motivated by recent experimental progress, we consider several scenarios, including different spin-orbit splittings in the two dots and including the case when the valley degeneracy is lifted due to an insulating ferromagnetic substrate. Finally, we discuss in which cases it is possible to reduce the low energy subspace to the lowest Kramers pairs. We find that in this case the low energy model is formally identical to the Heisenberg exchange Hamiltonian.


I. INTRODUCTION
Monolayers of transition metal dichalcogenides (TMDCs) are a class of 2D materials with very interesting electronic and optical properties [1,2]. They are atomically thin semiconductors with a direct band gap and two-fold degenerate valleys in the Brillouin zone [3,4]. Early studies discovered that their intrinsic spin-orbit interaction splits the spin states in the valence band and that it is possible to optically manipulate the spin and valley degrees of freedom (DOF) [5,6]. Further theoretical works suggested the idea that the conduction band should be spin-split as well [7,8] which found recent experimental confirmation [9].
The possibility to construct purely two-dimensional, electrostatically defined quantum dots (QDs) is one of the reasons that makes monolayer TMDCs so attractive, at least from a fundamental point of view. Compared to III-V semiconductors, such as GaAs [10,11], TMDCs have several isotopes with vanishing nuclear spin, thus lacking hyperfine interactions with the electronic spin. Moreover, TMDCs comprise an additional valley pseudospin. Although these two features are common to several other systems used for QDs, such as Si/SiGe quantum wells [12][13][14][15], graphene [16][17][18][19][20] and carbon nanotubes (CNTs) [21][22][23][24], TMDCs are special because they exhibit very strong spin-orbit coupling (SOC). Theoretical investigation of QDs in TMDCs started with the magnetic field dependence of the single-electron spectrum [7] and it now includes studies of valley hybridisation [25], flake QDs of triangular and hexagonal shape [26], the valley Zeeman effect [27], optical control of a spin-valley qubit in nanostructures [28], spin-degenerate regimes for small QD radius in an external magnetic field [29], a model of spin relaxation [30], electric control of a spin-valley qubit [31] and a model of valley qubit [32]. On the experimental side, gating monolayer TMDCs is not straightforward [34,36] and low material quality has hindered for a long time the experimental study of the intrinsic properties * alessandro.david@uni-konstanz.de of these materials. However, recently there has been a significant progress in the fabrication process of nanostructures in TMDCs. This has enabled the creation of single QDs on monolayer [33,37] or trilayer TMDCs [36], double QD experiments with tunable coupling strength between the dots [35,37] and the observation of gatecontrolled Coulomb blockade effect [33][34][35][36][37].
In this paper we study double quantum dots (DQDs) taking into account the spin-orbit coupling and both spin and valley DOF. We note that, although single [23] and double [38][39][40] quantum dots in CNTs have been studied, where the low energy theory shares some similarities with TMDCs, little is known about TMDC DQDs and the role of the exchange interaction. For simplicity we restrict each dot to the lowest orbital and we consider the system filled with two electrons. Starting from the model introduced in Ref. [41], we add the spin-orbit interaction and we find a low energy effective Hamiltonian for the case where each dot is occupied by one electron. We also investigate the situations where the spin-orbit splitting is different for each dot and where we include the effects of a magnetic field. Furthermore, for large spin-orbit splittings it is convenient to focus on a smaller subset of states formed by the two lowest Kramers pairs of the system. In this smaller subspace we find that in most cases the interaction is formally identical to the Heisenberg exchange interaction used to perform a CNOT gate between spin qubits [42]. This supports the idea that the lowest Kramers pair can serve as a qubit for TMDC, as was suggested in Ref. [7]. This paper is organised as follows. In Sec. II we introduce the basic terms appearing in our Hamiltonian, we present a naming conventions for certain useful projection operators and we briefly evoke the results for the null spin-orbit coupling case. Then, in Sec. III we present the effective Hamiltonians for the case where each QD is occupied by one electron, for three different situations: when the spin-splitting is equal in both dots, when the spin-orbit splitting is different and when the TMDC is deposited on an insulating ferromagnetic substrate or placed in an external magnetic field. Afterwards, we explore in Sec. IV under which conditions it is possible to focus on smaller effective Hamiltonians for these three arXiv:1802.09789v1 [cond-mat.mes-hall] 27 Feb 2018 situations. Finally, in Sec. V we present our conclusions.

A. Basic definitions
We give here the definitions for the different terms of the Hamiltonians that describe the various scenarios studied in this work. The operators c ( †) jτ σ annihilate (create) an electron in QD j with valley τ and spin σ. Here j = L (R) refers to the left (right) QD, τ = K (K) indicates the positive (negative) valley and σ =↑ (↓) specifies spin up (spin down).
The on-site Coulomb repulsion between electrons in the same QD is captured by the Hubbard Hamiltonian, where U > 0 is the positive charging energy of the dot and the number operator is defined as A detuning term specifies the energy difference ε between the dots, Electron-hopping from one dot to the other is accounted for by a tunneling term that preserves spin and valley, where the tunneling coefficient t is generally a complex number. The intrinsic spin-orbit coupling is modeled by a simple time-reversal symmetric (T -symmetric) spin-splitting: ∆τ z σ z [7]. Here, τ i (σ i ) is the i-th Pauli matrix acting on the valley (spin) DOF (i = x, y, z), while ∆ is a real, positive or negative, coupling constant. This implies that the Kramers pair of states in the set P = |K ↑ , |K ↓ is shifted by the energy +∆, while the Kramers pair of states in the set N = |K ↓ , |K ↑ is shifted by the energy −∆. We call P the positive Kramers pair and N the negative Kramers pair. For the double dot system, Eq. (5) assumes that the spin-orbit splitting is the same for every dot, which is usually the case for dots created on the same material. In case the spin-orbit splitting is different we use the following generalisation, where ∆ L and ∆ R are the spin-orbit splittings in the left and right QD respectively. We also consider the coupling of spin and valley to an external magnetic field. The corresponding spin Zeeman term is given by where h SL and h SR are two vectors of coupling constants for left and right QD respectively and σ is the vector of Pauli matrices acting on spin. The valley Zeeman term is the valley counterpart of the spin Zeeman but considering only the z-Pauli matrix acting on the valley. This is motivated by Ref. [7] and recent experiments [43,44], where h V Lz and h V Rz describe the valley splittings in the left and right dot. A concise outline of the interactions and the DOF of our DQD model is shown in Fig. 1. We consider the case where there are two electrons in the system. The possible charge configurations are (2, 0), (1, 1) and (0, 2), where (n L , n R ) means that there are n L electrons in the left QD and n R electrons in the right QD. Because of the spin and valley degrees of freedom and due to the Pauli exclusion principle, there is a total number of 28 linearly independent states: 6 (2, 0)-states, 16 (1, 1)states and 6 (0, 2)-states. Throughout this paper, we always assume a small detuning and weak tunneling, i.e. |t| |U ± ε|.

B. Naming basis states and projection operators
It is natural to refer to certain operators that will appear as projectors on specific states. We present a convention to name all 28 states of the total Hilbert space and we give the projection operators that will be useful later.
Not all the states in the (1, 1)-subspace are allowed to tunnel to (0, 2) or (2, 0)-states, but only those which are antisymmetric in spin and valley, because of the Pauli exclusion principle [40,45]. Following the antisymmetric nature of the tunneling states, we begin by introducing a basis consisting of states that are symmetric or antisymmetric in both spin and valley: |s V s S (n L ,n R ) , where (n L , n R ) is the charge configuration and s V (s S ) indicates the exchange symmetry of valley (spin) DOF, which can be either a singlet (S) or a triplet (T − , T 0 , T + ). See Table I for the definitions of singlet and triplet states for spin and valley. When (n L , n R ) = (1, 1) we omit the indication of the charge configuration in the subscript. The equal spin-orbit coupling term H ∆ (Eq. (5)) is not diagonal in this basis. In order to work with a basis that makes H ∆ diagonal, in the (1, 1)-subspace we substitute |ST 0 , |T 0 S , |SS and |T 0 T 0 with the following states, where |n ± (|n ± ) are antisymmetric (symmetric) spinvalley states consisting of only positive (subscript +) or negative (subscript −) Kramers pairs. Note that |n ± (|n ± ) are odd (even) under the time-reversal operator T . For these states there is no defined exchange symmetry for spin or valley alone. States analogous to |n ± are defined in the (2, 0) and (0, 2)-subspaces. In Table II, we list all these states grouped by charge configuration and by symmetry under exchange. Projection operators on each one of the antisymmetric (1, 1)-states are easy to obtain (see Appendix A), where τ ji (σ ji ) is the i-th Pauli matrix acting on valley (spin) in QD j = L, R and σ j , τ j are vectors of Pauli matrices acting on spin or valley respectively. We note that the projection operator over the whole antisymmetric subspace of the (1, 1)-sector can be written as [41] (see Appendix A), We finish this section by briefly considering the case where there is no spin-orbit splitting (∆ L = ∆ R = 0). The total Hamiltonian is, then, where H U , H ε and H t are defined in equations (1), (3) and (4) respectively. This model describes DQD systems with fourfold degenerate spin and valley states in each dot (and no spin-orbit interaction). This case has been extensively treated in Ref. [41]. For later reference and readability we report the resulting effective Hamiltonian, the exchange energy J in this case has the usual form, and the definition of P as is given in Eq. (12). This Hamiltonian shifts down the energy of all the (1,1) states in the antisymmetric subspace by −J. It means that the ground state is any state in this 6-dimensional subspace and it has energy −J. The first excited states are those in the orthogonal 10-dimensional symmetric subspace, with energy 0. This is represented in Fig. 2(a).

III. RESULTS FOR (1,1)-SUBSPACE
A. Symmetric spin-orbit splitting (∆L = ∆R = ∆) We first consider the case where the spin-orbit coupling is equal for both QDs: ∆ L = ∆ R = ∆. We call this symmetric spin-orbit splitting. The total Hamiltonian is where H ∆ is defined in Eq. (5). This Hamiltonian describes DQDs created in TMDCs [3][4][5][6][7]9], but it may also be used to study DQDs in CNTs [23, 38-40, 46, 47]. Symmetric spin-orbit coupling H ∆ shifts those states formed by two elements of the negative (positive) Kramers pair by −2∆ (+2∆), while it leaves unchanged those states formed by one element of the negative Kramers pair and one element of the positive Kramers pair (see Table II). In order to identify the (1, 1)-sector as our low energy subspace (LES) we have to guarantee first of all that no (2, 0) or (0, 2)-state is lower in energy than any (1, 1)-state. Looking at Table II we see that this condition is met when When both (9) and (17) are satisfied, it follows that (2, 0) and (0, 2)-states are energetically unfavored and tunneling out from the (1, 1)-sector is strongly suppressed. However, virtual tunneling processes must be taken into account. It is important to notice that antisymmetric (1, 1)-states can tunnel only to their (0, 2) and (2, 0) counterpart states that have the same spin and valley configuration. This is due to the spin-and valley-preserving nature of the tunneling term of Eq. (4): there is no transition from a negative Kramers pair to a positive Kramers pair and viceversa. Therefore, all the energy differences between initial and final states in a virtual tunneling process do not depend on ∆ and the exchange interaction that emerges does not change from the case of zero spinorbit splitting reported in Sec. II C. Exchange interaction and symmetric spin-orbit coupling act independently of each other. We obtain the effective Hamiltonian: where J is the same constant defined in (15), P as is given in (12), It is simple to evaluate the consequences of Eq. (18) since it is diagonal in the states of Table II. Assuming now that ∆ > 0, the ground state is |n − because it is the only state that is shifted down by both −J and −2∆, thus the ground state space is one dimensional, see Fig. 2(b),(c). The dimensionality of the ground state space changes from one to six only when ∆ is exactly equal to zero (see Fig. 2(a)) and this degeneracy is lifted linearly in ∆. As shown in Fig. 2(b),(c), there are five groups of excited states with various degrees of degeneracy and whose relative distances in energy depend on the sizes of ∆ and J. For later comparison, we observe that in this case each antisymmetric state is separated from its symmetric counterpart by the same amount, the exchange energy J. In the case ∆ < 0, the results are analogous, but the ground state is |n + .
Regarding the experimental situation, the few available experiments indicate U ≈ 2 meV for QDs created in WSe 2 [33] and MoS 2 [35][36][37] with QD radii around 100 nm. On the other hand, the theory predicts 2∆ to be about 3 meV for MoS 2 and larger for other compounds [7,48]. Therefore, at the moment, condition (17) seems not to be satisfied for TMDCs, but further experiments with smaller QD sizes may lead to higher charging energies. However, when Eq. (17) is not satisfied, it is possible to focus on a smaller subspace inside the (1, 1)-sector, as explained in Sec. IV.  (1)) and detuning (Eq. (3)) for the different charge configurations. The three columns on the right hand side show which states are shifted by −2∆, 0 and +2∆ by the action of the symmetric spin-orbit coupling H∆ defined in Eq. (5). The states are grouped as (1, 1)-states (antisymmetric and symmetric), (2, 0)-states and (0, 2)-states (only antisymmetric).

B. Asymmetric spin-orbit splitting (∆L = ∆R)
Suppose now that the spin-orbit strength is different in the two dots, ∆ L = ∆ R , this is the case of asymmetric spin-orbit splitting. The total Hamiltonian is where H ∆ has been substituted by H ∆ L ,∆ R defined in Eq. (6). There is theoretical evidence that for (rather small) QDs on TMDCs the spin-orbit strength ∆ depends on the radius of the QD [29]. Therefore, the Hamiltonian in Eq. (19) can reflect a situation where the two dots have different sizes. Changing the size of the QD would allow for a smooth and tunable adjustment of ∆. Another situation where ∆ L = ∆ R can be relevant is when the dots would be created in lateral heterojunctions such that each dot is created in different types of TMDC, which intrinsically possess different spin-orbit splittings [48]. Lateral heterojunctions of different TMDCs have already been demonstrated by several groups [49][50][51]. To our knowledge no experiment has been reported yet involving such a DQD. We notice that for heterojunctions the situation might be further complicated by mismatch in the energy band gap, acting as a detuning away from the interface, and by likely non-negligible differences in the Coulomb charging energy U of the two compounds. We point out that H ∆ L ,∆ R is not diagonal in the basis of the states listed in Table II. Each of the antisymmetric states |T + S , |T − S , |ST + , |ST − is mixed with its symmetric counterpart and the off-diagonal elements are ±(∆ L − ∆ R ). This is relevant to determine the ground state, but not for the form of the effective Hamiltonian, because the mixing happens inside the (1, 1)-subspace. Another difference is that, in contrast to Table II, the ±2∆ on the diagonal of the Hamiltonian are replaced by ±(∆ L + ∆ R ) in the (1, 1)-subspace, by ±2∆ L in (2, 0) and by ±2∆ R in (0, 2). When condition (9) is also valid, we can ensure that (1, 1)-states are lower in energy than (2, 0) and (0, 2) states by requiring that In this case we can use the Schrieffer-Wolff transformation [52][53][54] to obtain the effective Hamiltonian. The term H ∆ L ,∆ R has three significant consequences. Firstly, the following matrix element differences now depend on where E |ψ = ψ| H tot |ψ . Secondly, in order to apply perturbation theory to the (1, 1)-subspace, the coupling constants must satisfy the conditions |t| |U ± ε|, |t| |U ±(ε+∆ L −∆ R )| and |t| |U ±(ε−∆ L +∆ R )| simultaneously. In turn, these imply that |∆ L − ∆ R | U . Lastly, the exchange energies relative to |n ± acquire a dependence on the asymmetry of the spin-orbit splittings, c.f. Eq. (15). The exchange energies for states |T + S , |T − S , |ST + and |ST − are not affected. An effective Hamiltonian can be written in a compact form in the following way, where P |n− and P |n+ are the projectors on |n − and |n + defined in equation (11c) and in this case we write the restriction H ∆ L ,∆ R | (1,1) = ∆ L Σ L +∆ R Σ R , with Σ j = τ jz σ jz . By diagonalizing the Hamiltonian (23), we obtain the energy levels shown in Fig. 3. Interestingly, the two levels that appear at energy 0 and energy −J for the symmetric spin-orbit splitting case move to the energies, This is a consequence of the level repulsion due to the off-diagonal ±(∆ L − ∆ R ) matrix elements coming from H ∆ L ,∆ R . However, the degrees of degeneracy remain the same. Looking at Eq. (24) and at Fig. 3, we realise that, for positive ∆ L + ∆ R , the ground state is |n − , unless |∆ L −∆ R | is larger than ∆ L +∆ R , in which case the states |T + S , |T − S , |ST + , |ST − (or their symmetric counterparts) span a 4-fold degenerate ground state space. We want to briefly mention here that a generalisation of the asymmetric spin-orbit coupling which preserves the T -symmetry, including couplings to the in-plane components of the spins, does not affect the exchange interaction. Consider (25) where ∆ j is the vector of the three coupling constants ∆ jx , ∆ jy , ∆ jz for QD j, multiplied by the vector of spin Pauli matrices σ. Eq. (25) resembles the spin Zeeman term (cf. Eq. (7)), but it also preserves T -symmetry because of the τ z that multiplies the spin Pauli matrices. Qualitatively, a term like this should appear in the Hamiltonian of TMDC DQD systems for which it is possible to define an average local tilt and an average local curvature of the TMDC sheet for each dot [55]. However, tilt and curvature add couplings between conduction and valence band as well and one would need to check whether these couplings are small enough to be neglected. This is not the focus of this work and we only want to stress that substituting H ∆ L ,∆ R with H ∆ L ,∆ R ,xyz in Eq. (19) does not modify the exchange interaction given in Eq. (23).

C. Spin and valley Zeeman and symmetric spin-orbit coupling
In this section we investigate the influence of spin and valley Zeeman terms on the symmetric spin-orbit splitting case. The total Hamiltonian for this scenario reads where H S and H V are defined by equations (7) and (8) respectively. This is the case depicted in Fig. 1. This model approximates at least two similar but distinct situations. First, when a monolayer TMDC is placed in a non-uniform and weak magnetic field, the spin of the electron in each dot couples to the local field in all three spacial directions through the spin Zeeman interaction H S . On the other hand, we assume that the valley only couples along z, the orthogonal direction to the plane [7]. In contrast to Ref. [7] here we only consider the lowest orbital state in each of the dots but we also take into account the Coulomb charging energy U . The other possible situation that came to our attention is that of monolayer TMDC deposited on a ferromagnetic insulator, e.g. europium oxide (EuO). Such a situation was considered in Ref. [56], where Eq. (7) was used to model the giant and tunable valley splitting. A Rashba term is also included, mixing conduction and valence band, that we do not consider here. Although Ref. [56] does not consider the valley Zeeman term, we note that it is allowed in the effective Hamiltonian of the TMDC because the ferromagnetic substrate breaks time-reversal symmetry. The presence of such magnetic exchange field in TMDC and ferromagnetic semiconductor heterostructures has been recently demonstrated [43,44]. Now both H S and H V are non-diagonal in the basis of states presented in Table II and they mix a large number of states. We do not list all the newly introduced matrix elements but we observe that they are confined inside the charge sectors (0, 2), (1, 1) and (2, 0), as H S and H V do not contain terms of the form c † Lτ σ c Rτ σ (or h.c.). The only coupling between (1, 1)-states and (2, 0) or (0, 2)-states remains tunneling. Nevertheless, we have to ensure that (1, 1)-states are lower in energy than other states. In Appendix B we report the discussion of which conditions must be fulfilled in order for this to be the case. In addition, to apply the Schrieffer-Wolff transformation, Eq. (9) and the conditions where = V, S, must also be satisfied. The matrix element differences of states |T − S , |T + S , |ST − , |ST + with their (2, 0) and (0, 2) counterparts acquire a dependence on h V Lz − h V Rz or h SLz − h SRz . Defining the exchange energies the effective Hamiltonian can be written as where J is the standard exchange energy defined in Eq. (15) and 1) .

IV. RESULTS FOR NEGATIVE KRAMERS PAIRS SUBSPACE
In this section we explain in which cases it is possible to restrict the LES to the tensor product of the two lowest Kramers pairs in energy and we present the effective Hamiltonians in this smaller subspace. In general, this happens when the spin-orbit strength is quite large, thus the results of this section are relevant for TMDCs. From now on we assume ∆ L,R > 0, so that for each dot the lowest Kramers pair is N , the negative one (see Sec. II A). Similar results can be obtained in case the spin-orbit splittings are negative and the lowest Kramers pair is P, the positive one.
A. Symmetric spin-orbit splitting (∆L = ∆R = ∆) In the case of symmetric spin-orbit splitting (∆ L = ∆ R = ∆, see H tot in Eq. (16)) there is no matrix element that allows a transition between N and P. Assuming Eq. (9) and 2∆ > J (which is usually the case for TMDCs at low detuning) we see from Fig. 2(c) that the subspace spanned by the states of N ×N (also called N ×N -sector) is the LES, where, . Now we present the effective Hamiltonian when the system is restricted to this LES.
We can easily identify the negative Kramers pair N as a spin-1/2 DOF with so that the effective basis N × N can be seen as the basis of all the states of two spins, In this basis, the spin-orbit coupling always assumes the value of −2∆, which we ignore in the following since it is just an energy shift. It is well known that the effective Hamiltonian for two electrons with only the spin DOF in two QDs with only tunneling, detuning and Hubbard potential is given by [11] H eff = −JS L · S R , with the usual exchange energy J multiplied by the projector on the singlet state |S , the only antisymmetric state for two spins (S j is the vector of spin operators acting on QD j). Analogously, using the following operators, the Hamiltonian of the system restricted to basis N × N is where S ji = 1 2 σ ji is the spin operator proportional to the new Pauli operator σ i , i = x, y, z, acting on QD j = L, R. The ground state space is one-dimensional as for the full (1, 1)-subspace. Indeed, the ground state is the same, Very similar considerations are valid when we relax the constraint of equal spin-orbit splittings in the dots (∆ L = ∆ R , see H tot in (19)). To reduce the LES to the subspace spanned by N × N we need see Fig. 3 and Eq. (24). Of course, to apply perturbation theory, the conditions |t| |U ± (ε + ∆ L − ∆ R )| and |t| |U ± (ε − ∆ L + ∆ R )| must be valid here too. Following the same steps that led to Eq. (34), we arrive to the effective Hamiltonian, where J |n− is given in Eq. (22). Indeed, from Fig. 3 we see that the energy levels for the N ×N -sector are similar between symmetric and asymmetric spin-orbit splitting, only the exchange energy J is altered.
To better understand what is the effect of ∆ L − ∆ R in J |n− with respect to J, we may write Since (∆ L − ∆ R )/(U ± ε) 1, we can write The effect of an almost symmetric spin-orbit coupling on the exchange energy is that of fine tuning around the value of J by a small positive or negative quantity, depending on the sign of ε and ∆ L − ∆ R . Moreover, as ε → 0, J |n− → J. The above results remain valid for the generalisation to H ∆ L ,∆ R ,xyz (Eq. (25)).

C. Spin and valley Zeeman and symmetric spin-orbit coupling
We now consider the total Hamiltonian of Eq. (26) restricted to the N × N -sector. The results in this section are valid when the spin and valley Zeeman coupling constants are weak compared to the spin-orbit strength: |h Sji |, |h V jz | 2∆, j = L, R, i = x, y, z. We show directly and discuss the effective Hamiltonian for the subspace spanned by N × N as a 4 × 4 matrix. The columns are associated, from left to right, to the rotated basis states |n − , |n − , |T + T − and |T − T + , where J is the standard exchange energy as in Eq. (15). We use the notation with = S, V , j = L, R and A ± and A are given by Since |h Sjz | 2∆, we can use the approximation A + ≈ A − ≈ A, thus the contribution of A to Eq. (40) is an energy shift that can be ignored. Therefore, in this case we can also write where H S | N ×N and H V | N ×N are the restrictions of H S and H V to the N × N -sector.
Given the generality of H ∆ L ,∆ R ,xyz (Eq. (25)) as a Tsymmetric term and the fact that H ∆ L ,∆ R ,xyz and H S only differ in the presence of τ z that multiplies the spin Pauli matrices, we can assert that the new matrix elements (H z S , H z V , δh z S , δh z V , A ± and A) that appear in Eq. (40) as compared to Eq. (34), are allowed only by breaking the time-reversal symmetry.

V. CONCLUSIONS
In this work we have studied and presented the influence of a spin-orbit coupling on the low energy properties of a DQD system with spin and valley DOF in the (1, 1) charge configuration. In our analysis we have also explored the possibility of a different spin-orbit splitting in each dot and we have included a T -symmetry breaking magnetic field. In addition, we have discussed under which conditions the LES corresponding to the (1, 1) charge configuration, which is 16-dimensional, can be further restricted to a 4-dimensional subspace (the N × Nsector). We found that an equal spin-orbit splitting in each dot has no effects on the induced exchange interaction with respect to the case without spin-orbit coupling. On the other hand, asymmetric spin-orbit splitting, spin Zeeman and valley Zeeman modify the exchange coupling constants of three different pairs of antisymmetric states respectively. The modification of the exchange energies is similar for all these three pairs of states and it depends on the asymmetry of the interaction between left and right QD.
We also found that TMDCs satisfy the conditions to restrict the LES to the N × N -sector, where the effective Hamiltonian for the symmetric spin-orbit splitting case is formally identical to the Heisenberg exchange interaction between two spin-only qubits in valley non-degenerate materials. This renders the Kramers pair an ideal implementation of a qubit in TMDC, as was suggested in Ref. [7]. If the τ x operation can be effectively implemented (theoretical proposals to achieve this include the use of impurities [31,57] or the use of oscillating confinement potentials [32]), a recipe for a CNOT gate with these states is readily available from the original Loss and DiVincenzo proposal for spin-only qubits [42]. Moving to the asymmetric spin-orbit splitting we found that only the exchange energy is affected, while the form of the exchange Hamiltonian remains unchanged. The spin-orbit coupling asymmetry offers a way to tune the exchange energy other than the detuning of the dots. Finally, for the spin and valley Zeeman case, new couplings appeared in the reduced effective Hamiltonian, which originate from the breaking of the time-reversal symmetry.
with |τ 1 σ 1 ; τ 2 σ 2 = c † Lτ1σ1 c † Rτ2σ2 |0 . The operators that project on these states are shown in Eq. (11). We can follow two approaches to obtain them, one is more intuitive and gives compact results, assembling projectors on larger parts of the Hilbert space, but in practice it works only when the states we are considering have a structure for which we already know the correct basic projectors. The other one is rather formal but general, however it does not give the projectors in a compact form.
In order to show that equations (11) are really the projectors we are looking for, we will present here the intuitive way to derive them, using symmetries in the structure of states (A1). Take τ z , this operator has eigenvalue +1 when the state is in valley |K and eigenvalue −1 when the state is in valley |K . Thus, (1 + τ jz )/2 is the projector on all those states which have an electron in the j-th QD in valley |K and (1 − τ jz )/2 projects on states where the electron in QD j has valley |K . Analogous considerations hold for the spin operator σ z and the spin states |↑ , |↓ . In other words, we wrote down projectors on states which have a certain valley or possess a certain spin in a specific dot. Now focus on τ Lz τ Rz , this string of operators has eigenvalue +1 when both electrons are in the same valley (|KK or |KK ) and eigenvalue −1 when the valleys are different (|KK or |KK ). Then, operators (1 ± τ Lz τ Rz )/2 project on states whose valleys are the same (+) or are opposite (−). Similar considerations are valid for their spin counterparts.
The states with the simpler structure are |T + S , |T − S , |ST + and |ST − . They all have one of the properties fixed for both dots (either positive valley, negative valley, spin up or spin down, respectively). Consider |T + S (Eq. (A1a)), for both dots the valley is |K , but the spins are in a singlet state. The combination projects on all the states with positive valley in both dots.
To complete the expression we multiply by the projector on the spin singlet, given by (1 − σ L · σ R )/4. Finally, the complete projector is as in (11a). The other very similar expressions in (11a) and (11b) follow the same derivation with the appropriate changes of signs and spin/valley operators. States |n ± of equations (A1e) and (A1f) have a more complicated structure. First of all, they are both composed of states with opposite valley and opposite spin in the dots. The operator projects on the subspace spanned by {|K ↑; K ↓ , |K ↓; K ↑ , |K ↓; K ↑ , |K ↑; K ↓ }.
To distinguish between |n + and |n − , we note that for |n + the product of the eigenvalues of spin and valley inside each dot is positive, while for |n − it is negative (this is why they get an energy shift of ±2∆ by the spin-orbit coupling H ∆ ). We use τ Lz σ Lz to distinguish them and the projectors become positive sign for the subspace spanned by {|K ↑; K ↓ , |K ↓; K ↑ } and negative sign for the subspace spanned by {|K ↓; K ↑ , |K ↑; K ↓ }. Finally, the operator that identifies the correct superposition is τ Lx σ Lx τ Rx σ Rx , the same for both states. The total projectors for |n ± are as in Eq. (11c). Formally, the above results can also be obtained by noticing that the projectors we are looking for are spinand-valley operators (i.e. acting on the (1, 1)-subspace) and every spin-and-valley operator P can be expanded in a linear combination of Pauli operator strings, P = i,j,k,l=0,x,y,z P ijkl τ Li σ Lj τ Rk σ Rl , where P ijkl are complex coefficients. If |ψ is the state we want to project on, the projector is P = |ψ ψ| and we can use the inner product provided by the trace operation to obtain the expansion coefficients, P ijkl = Tr [|ψ ψ| τ Li σ Lj τ Rk σ Rl ] .
Again, recollecting a compact expression from this set of coefficients requires some work (even when most of them are zeros) and does not give any physical insight. For completeness, we briefly discuss P as , the projector on the whole antisymmetric subspace of (1, 1)-states of Eq. (12). Although this operator could be obtained as the sum of the 6 projectors of Eq. (11) on the single, linearly independent antisymmetric states, we can write it down following another procedure, derived in Ref. [41]. It can be expressed as a combination of projectors on the singlet (S, antisymmetric) and on the triplet (T , symmetric) subspaces for spin and valley, P as = P S spin P T valley + P T spin P S valley .