Spin-valley system in a gated MoS$_2$-monolayer quantum dot

The aim of presented research is to design a nanodevice based on a gate-defined quantum dot within a MoS$_2$ monolayer in which we confine a single electron. By applying control voltages to the device gates we modulate the confinement potential and force intervalley transitions. The present Rashba spin-orbit coupling additionally allows for spin operations. Moreover, both effects enable the spin-valley SWAP. The device structure is modeled realistically, taking into account feasible dot-forming potential and electric field that controls the Rasha coupling. Therefore, by performing reliable numerical simulations, we show how by electrically controlling the state of the electron in the device, we can obtain single- and two-qubit (thus universal) gates in a spin-valley two-qubit system. Through simulations we investigate possibility of implementation of two qubits \textit{locally}, based on single electron, with an intriguing feature that two-qubit gates are easier to realize than single ones.


I. INTRODUCTION
Two-dimensional crystals consisting of single layers of atoms are modern materials that can be used for implementation of quantum computation. 2D monolayers of transition metal dichalcogenides (TMDCs), e.g. MoS 2 , seem to be better candidates than graphene because of their wide band gaps and strong electrically induced spinorbit coupling of the Rashba type 1,2 . By considering the valley degree of freedom of an electron together with its spin we extend our ability to define a qubit into: spin, valley 3 and hybrid spin-valley qubit 4 . However, the most interesting is the definition based on spin and valley of a single electron as a two-qubit system 5,6 .
The area of application of monolayer materials for construction of electronic nanodevices is currently under strong development 1,7-12 . Methods for building devices based on gated TMDC monolayers or nanotubes become increasingly advanced [13][14][15][16][17] , opening the possibility of utilizing the spin and valley index of electrons controlled therein. In particular, it is shown by recent results with electrostatic quantum dots (QDs) with a gated MoS 2 -nanoribbon-QD measured by a single electron transport 18,19 , or tunable TMDC spintronic devices, where spin or valley polarized currents emerge in TMDC monolayer proximitized by nearby ferromagnetic [20][21][22][23][24] . Transistor structures with TMDC monolayer forming active area in tunnel FETs are being developed 25 , also with vertical TMDCs heterostructures 26,27 . The more intriguing lateral, in-plane TMDCs heterojunctions are also constructed, enabling interesting 1D physics at interfaces 28 , or leading to improved FETs switching characteristics [29][30][31] .
Inspired by this, we examined the possibility of realization of a nanodevice based on a MoS 2 monolayer, capable of creating a two-qubit system defined on spin and valley degrees of freedom of a confined electron. For this purpose, we have built a realistic model of the nanodevice and perform numerical simulations that prove its capa-bilities. Thanks to the use of appropriately modulated local control voltages, the system is all-electrically controlled and does not require using photons or external microwaves, thus significantly improving its scalability.

II. MODEL
In this section we will go through the device model. The potential in the entire nanodevice, controlled by the gate voltages, is calculated by solving the Poisson equation, while the electron states in the flake are described with the tight-binding formalism. Let's start with the nanostructure overview.

A. Device structure
The proposed nanodevice structure is presented in Fig. 1. On a strongly doped silicon substrate we place a 20-nm-thick layer of SiO 2 . Then we place two electrodes The proposed nanodevice structure containing a MoS2 monolayered flake deposited on a SiO2 layer, and separated by hBN from the control gates responsible for creating the QD confinement potential. which serve as a source (S) and a drain (D). Directly on them we deposit a MoS 2 -monolayer (hexagonally shaped) flake of 16-nm-diameter. The monolayer is then covered with a 5-nm thick insulating layer of hexagonal boron nitride (hBN) with a large bandgap 32 , forming a tunnel barrier. Finally on top of the sandwiched structure we lay down four 15-nm-wide control gates (G 1..4 ), placed symmetrically around the central square-like gap of size 20 × 20 nm. The gate layout presented here is quite similar to the one proposed by us recently 3 , but with a larger 20 nm clearance between opposite gates, which may ease their deposition.
Source, drain and the gates layout are clearly presented in Fig. 1. Voltages applied to these gates (relative to the substrate) are used to create confinement in the flake. To calculate realistic electrostatic potential φ(r) we solve the Poisson equation taking into account voltages V 1..4 applied to control gates G 1.. 4 and to the highly doped substrate V 0 = 0, together with space-dependent permittivity of different materials in the device 3,33 . Resulting potential in the area between SiO 2 and hBN layers, where the flake is sandwiched, is presented in Fig. 2. We de- plete the electron gas until a single electron remains in the formed dot confinement potential.

B. Monolayer model
The monolayer flake is made of molybdenum disulfide. MoS 2 monolayers are successfully described by several tight-binding (TB) models, with different numbers of orbitals used, including nearest or next-nearest neighbors. Seven 34 or eleven [35][36][37] Mo and S orbitals construct TB basis to reproduce low-energy physics in the entire Brillouin zone, also near the Γ-point. Although the simpler three-band (including three Mo orbitals) TB model 38 FIG. 3. The MoS2 monolayered flake structure: (left) hexagonal flake employed in the device has sides made of NB = 25 Mo atoms, giving the flake side of 8 nm (lattice distance between Mo nodes is 0.319 nm); (right) MoS2 crystal lattice structure formed of hexagonally packed Mo and S atoms arranged in triangular lattices rotated relative to each other by π. The Mo lattice vectors R k determine the hopping directions in our nearest-neighbors TB model. fails around the Γ-point, it correctly represents the orbital composition around the K point near the (both conduction and valence) band edges, where the Bloch states mainly consist of Mo d orbitals 39 . Thus it is good enough to deal with low energy states near the band minimum. However when considering perpendicular electric field, crucial for the Rashba coupling, we have to include also S orbitals localized above and below Mo plane in the MoS 2 structure (see Fig. 3). For calculating Rashba coupling we will utilize 11-band model with three p orbitals for each S atom in dimer. 35 Consequently, we have described the monolayer structure using three Mo orbitals: d z 2 , d xy , d x 2 −y 2 , and the nearest-neighbors hoppings 38 : The potential energy of the electrostatic confinement at the i-th lattice site: ϕ i = −|e|φ(x i , y i ) together with the on-site energies α enter the on-diagonal matrix elements (α numbers the orbitals). The off-diagonal electron hopping element from the β Mo orbital localized in the j-th lattice site to the α orbital localized in the i-th site is denoted by t αβ ≡ t αβ R k(i,j) . It depends on the hopping direction (between ij neighbor pair) described by the nearest neighbor vectors R k for the molybdenum (Mo) lattice, which are defined as in Fig. 3. They form two non-equivalent families: R 1 , R 3 , R 5 and R 2 , R 4 , R 6 with the nearest sulphur (S) neighbor on the left or right side, as shown in Fig. 3. This symmetry constraint reflects on the reciprocal lattice where in the corners of the first (hexagonal) Brillouin zone, the K points form two non-equivalent families: K and K .
Opposite hoppings are mutually transposed: t αβ (R k ) = t βα (−R k ). Their explicit forms, together with the on-site energies α , can be found in [3 and 38].

C. Rashba coupling
Electric field perpendicular to the monolayer surface breaks the reflection σ h symmetry and modifies the onsite energies of atoms in three MoS 2 sublayers. This leads to externally, electrically controlled spin-orbit interaction (SOI) of the Rashba type. The Rashba coupling can be also introduced to layered TMDCs by a structure asymmetry from ferromagnetic substrate leading to the proximity effect 20 .
The idea to calculate the electrically induced Rashba spin-orbit coupling strength is to take the tight-binding model with atomic spin-orbit coupling (introduced by λ αβ in Eq. 1) including also p orbitals of the sulfur top and bottom sublayers to which we apply on-site potentials V t and V b 34,[40][41][42] . The difference between them results from external electric field: . While d = 0.32 nm is the monolayer thickness (sulfur sublayers distance).
Given such extended tigh-binding model (we take 11orbital model of Ridolfi et al. [35]), we perform downfolding, using the Löwdin partitioning technique, to our Mo-orbitals model and obtain the Rashba coupling γ R within this 3-band base. Further details of the calculation are attached in the Appendix. Resulting coupling matrix elements are proportional to external electric field: γ αβ R (x, y) = |e|E z (x, y) γ αβ , with the explicit form The tight-binding Rashba Hamiltonian 42-46 (in Eq. 1) with characteristic spin-and orientation-dependent hopping between two nearest neighbor bonds is: with the Pauli-matrices vector s = (σ x , σ y , σ z ), and the (unit) versorê ij ≡ (e x ij , e y ij ) pointing along the bond connecting sites i and j. An obvious propertyê ji = −ê ij ensures H R hermiticity. The expanded hopping expression is (ê ij ×ẑ)· s σσ = e y ij s x σσ − e x ij s y σσ , where s x = σ x , s y = σ y , and γ R = |e|E z γ with the above defined γ.

D. External magnetic field
To include electron interaction with a perpendicular magnetic field in the monolayer model, we should add to the Hamiltonian a standard Zeeman term (in Eq. 1): with a magnetic field B. For γ Z = geµ B 2 we arrive at the standard Zeeman energy gµ B 2 s · B. To address also orbital effects related to magnetic field we apply the so-called Peierls substitution 47 . We multiply the hopping matrix, by the additional factor t ij →t ij = t ij exp (ıθ B ) in the Hamiltonian (1). Now the vector potential enters Eq. (1) via the Peierls phase θ B , calculated as the path integral between neighbor nodes: A is a vector potential induced by the B field. We use the Landau gauge, with the vector potential A(r) = [0, dxB z (x, y), 0] T for the perpendicular magnetic field B(r) = [0, 0, B z (r)] T . This leads to the phase: The most important result of applying a magnetic field is a splitting introduced between levels with opposite spin and valley index. Interestingly there are two types of Zeeman splittings: standard, spin type, and Zeeman valley splitting, both presented in Fig. 6 in Section IV. Each of them possesses other Landé factor. This will enable us to separately address each transition between four basis states in the spin-valley two-qubit space.

III. CALCULATION METHOD
Let's have a look at the stationary and time-dependent calculation methodology. Firstly, we solve the eigenproblem for the stationary Hamiltonian (1): H(r)ψ m (r) = E m ψ m (r), and obtain M eigenstates. For our hexagonal flake of size N B = 25 we have 1801 sites × 6 giving M = 10806 eigenstates. For the Hamiltonian matrix eigenproblem we utilize the fast and efficient FEAST routine 48 . The obtained eigenstates are represented by 6-dimensional vectors ψ m (r) = (ψ σα m (r)) , with σ = 1, 2 and α = 1, 2, 3. They belong to te state space H spin 2 ⊗H oribital 3 , with the spin and the 3-dimensional Moorbitals space. To identify them, at first we need an electron density calculated as ρ(r) = |ψ m (r)| 2 to determine if given state is localized at the flake edge forming the so-called edge state, or is confined within the quantumdot. Secondly, we need to identify the state quantum numbers-valley and spin indices. To do this we utilize similar formulas as in Eqs. (12) and (13), here adapted to a stationary state ψ m (r).
During the time-dependent calculations we will be working in the previously found eigenstates base. Therefore the full time-dependent wave function is represented as a linear combination of N basis states ψ n : together with time-dependent amplitudes c n (t) and phase factors of the corresponding eigenvalues E n . We assume a basis of N = 200 < M lowest eigenstates from the conduction band (represented by yellow bullets in the lower inset in Fig. 5). The time evolution is governed by the time-dependent Schrödinger equation: with the time-dependent Hamiltonian being a sum of the stationary part (Eq. 1) and a time-dependent contribution to both the potential energy and the Rashba coupling: The full time-dependent potential energy ϕ(r, t) = ϕ(r)+ δϕ(r, t) contains variable part δϕ(r, t), generated by modulation of the gate voltages. Whole is calculated as ϕ(r, t) = −|e|φ(r, t), with the potential φ(r, t) obtained by solving the Poisson equation for the variable density ρ(r, t) at every time step. Note that the charge density originates from the actual wave-function, thus the Schrödinger and Poisson equations are solved in a selfconsistent way. Similarly, the time-dependent part of the Rashba coupling δH R (r, t) is induced by the variable part of the electric field: (9) with the same formula. Insertion of (7) to the Schrödinger equation (8) gives a system of equations for time-derivatives of the expansion coefficients at subsequent moments of time: The actual matrix elements δ mn (t) = ψ m |δϕ(r, t) + δH R (r, t)|ψ n need to be calculated at every time step due to changes in the potential and the electric field. Then, by using it, we solve the system (10) iteratively using a predictor-corrector method, with explicit "leapfrog" and implicit Crank-Nicolson scheme, obtaining the next time step of the system evolution.
For the electron wave function Ψ(r, t) we calculate the Fourier transform: on the flake surface area F , with 2D-wave vector k ≡ (k x , k y ). The Fourier transform naturally has periodic structure in the reciprocal space, therefore we can limit the k-area toF : k x,y ∈ − 2π a , 2π a , encompassing the (first) Brillouin Zone (BZ). KnowingΨ(k, t) we can calculate density in the reciprocal space expressed as: ρ(k, t) = |Ψ(k, t)| 2 . The k-density calculated for the |K ↓ state from Fig. 5 is presented in Fig. 10(0). We also mark the BZ along with points of high symmetry: Γ in the center, two types of K(K ) at the corners of the hexagonal zone and M on the edges of the hexagon. The coordinates of the high-symmetry points are: Γ = (0, 0), one of K = π a ( 4 3 , 0) and one of M = π a (1, 1 √ 3 ), with the lattice constant a = 0.319 nm. We can clearly see in Fig. 10(0) that density peaks are localized in the neighborhood of the K points (while not next to K ), confirming that in the |K↓ state exactly K valley is occupied. Now, the valley index K is calculated as: on the reciprocal space areaF 1/3 defined as two opposite π/3 sectors (withinF area) encompassing exactly one K point and one K point, i.e. |ϑ| ≤ π 6 ∪ |ϑ| ≥ 5π 6 , with azimuthal angle ϑ = atan2(k x , k y ). Because K(K ) point The electron spin value s is calculated as the expectation value of the Pauli z-matrix s z = σ z : also integrated on the flake area F for the actual wave vector Ψ(r, t). Operation ⊗ 1 3 means that during the spin calculations we trace out over the orbitals subspace.

IV. ELECTROSTATIC QUANTUM DOT
By applying voltages V 1..4 = 1500 mV to all of the gates we form the QD potential energy in the flake area, presented in Fig. 4(a). Calculated electronic eigenstates of the Hamiltonian (1) for the entire flake lattice, forms a ladder in Fig. 5 representing subsequent M = 10806 eigenstates ψ m (r). The bullets color is used to mark the dot occupation, namely the brighter the color is, the electron is more localized in the flake center. E.g. the yellow states are strongly confined, while black color marks the edge states with density localized on the flake border. These states are inaccessible to the electron confined in the QD, forming a forbidden energy range, namely, a bandgap. The bandgap divides the QD eigenstates into Subsequent eigenenergies of the flake Hamiltonian (1) for an electron confined in the created quantum dot, marked as color bullets. The bullet-colors describe the QD occupation, with yellow states for carriers strongly confined and localized at the flake center, while black bullets for the edge states with density at the flake border. The latter levels are forbidden for a carrier confined within the dot and define the energy bandgap.
conduction and valence bands. Lets now zoom into CB minimum. The states therein are presented in insets from Fig. 5.
The first four states form two doublets {|−1, ↑ , |1, ↓ } and {|−1, ↓ , |1, ↑ }, spin-orbit split (see the upper inset in Fig. 5). Their electron density is presented in Fig. 4(b). It turns out, that (states from) both the bottom of the conduction band and the top of the valence band are located at the points K and K (we have a direct band gap here), not at the Γ point. These bands form two non-equivalent valleys K and K which can be occupied by qubit carriers. Subspace spanned by the first four states consists of exactly one valley and one spin two-level system, forming together a 4-dimensional Hilbert space H valley 2 ⊗ H spin 2 of spin-valley two-qubit states |K, s .

A. Two-qubit subspace
If we add an external magnetic field, degeneracy in both pairs is lifted, as presented in Fig. 6 for B z = 1 T. Calculated splitting, here for 1 T (for 2 T resulting factors are the same), between first pair: is 242 µeV, which is in agreement with differences 238 µeV between opposite spin states (with different valleys) in the conduction band minimum for a larger dot [49]. The difference of 242 µeV leads to g v − g s 4.18. For the second pair 84. Therefore, obtained effective valley and spin g-factors are: g v = 6.51, similar to DFT calculations giving the value 7.14 [2, and errata : 50]. While the spin splitting factor g s = 2.33, with agreement with the DFT calculations [2] and experimental result [51], both giving value about 2.2. If we now take into account both spin and valley splitting, it turns out that the higher states pair is more split than the lower one, as presented in Fig. 6. This results from emergence of additional valley Zeeman splitting, for which K levels bend upwards and for K downwards (blue arrows). Similarly, spin-down level bends down, while spin-up bends up (red arrows). Therefore, for the higher levels pair the splittings will add, while for the lower one they will subtract. Thanks to such a form of level splittings, all of the transitions between them (6 in total) can be separately addressed.

B. Intervalley coupling
Let's now examine the intervalley coupling strength and its origin. The doubled intervalley coupling 2Λ can be calculated from the difference between the ground state and the 1st excited state in CB with no spin-orbit interaction 52 . Further, we assume that the intervalley coupling is: The structure of electronic states confined within a nanoflake and mostly presence of edge states, that cross the gap, depends on the flake edge type [53][54][55] . Zigzag edges do not mix valleys, but supports edge states. On the other hand, gap states are missing for an armchair edge type, which mixes valleys and induces transitions in graphene-like structures 56 . Same edge-dependent valley mixing was proven for MoS 2 nanoribbons 57 . Here, to skip the edge influence, we take a flake with a zigzag edge, and induce confinement strong enough to decouple the electron from the flake edge.
To check completely the edge influence on the intervalley mixing, we calculate 2Λ as a function of the flake size N B and the confinement depth, controlled by V dc . The results with no confinement are presented in Fig. 7  It turns out that the approximate relation between the coupling and the applied voltage in this range of the flake sizes is 2Λ ∼ √ V . It is similar to the relation between eigenenergies and potential of a harmonic oscillator. Now the coupling is purely confinement-dependent and does not depend on the flake size. This is because the confined electron is decoupled from the edge and its valley mixing is controlled electrically via the confinement potential.
A 10-nm-scale lithographic process required in in Fig. 1 can be difficult to achieve. Our calculations show that scaling up the structure (with preserved gate voltages) will decrease the intervalley coupling amplitude as 1/L, where L is the gate width. This means that after magnifying gates layout by an order of magnitude the coupling still keeps reasonable values. However, such scaling proportionally extends the intervalley transition time. confinement potential in a way that the dot minimum oscillates back and forth in the x-direction. The potential energy landscape at oscillations start, i.e. at t = 0, is presented in Fig. 8(b). While its form at the maximum left and right displacement from the center position, e.g. at t = π 2ω and t = 3π 2ω , is presented in Figs. 8(a) and 8(c) respectively. Additionally to inducing oscillatory movement of the dot position, the potential shape is modulated and becomes narrower at the maximum shift to the left-see Fig. 8(a), while shallower at the maximum right- Fig. 8(c). This enforces oscillatory squeezing of the electron state density, as seen in Fig. 8(g-i).

A. Spin and valley transitions
The confinement potential modulation introduces two effects to the system. Firstly, the voltage modulation moves the electron confined in the QD potential in an oscillatory way. The electron position oscillations causes that its momentum also oscillates. The velocity defined as the time derivative of the electron expectation position d dt x is presented in Fig. 9 as the red curve, while the electron position as the orange one. Together with the present perpendicular electric field E z felt by the electron (blue curve in Fig. 9) which induces the Rashba SOI, it creates spin-orbit mediated electron spin resonance transitions. The mean electric field E z is almost constant during oscillations, which is related to fact that the perpendicular electric field component E z (x, y) is mostly uniformly distributed on the flake, as seen in Fig. 8(d-f).
Secondly, oscillatory shallowing of the confinement potential leads to electron packet squeezing, visible as os- cillations of the electron packet size σ in Fig. 9 (green curve) and causes intervalley coupling changes. Resonant modulation of the intervalley coupling generates gradual transitions of the electron between the different valley states 3 . Subsequent stages of a transition between the different valleys in reciprocal space are presented in Fig. 10. Initially, the electron density in the reciprocal spacẽ ρ(k x , k y ) is localized in the K point vicinity within the BZ, as showed in Fig. 10(0). The voltage pumping process with the resonant frequency ω, tuned to the energy spacing ω between the states |1, ↑ and |−1, ↑ (or |1, ↓ and | − 1, ↓ ), leads to a gradual change of occupation to K valley, with density flow between valleys visible in the subsequent stages- Fig. 10( π 4 ) and ( π 2 ). After time t = 345 ps the electron occupies the K valley entirely (Fig. 10(π)) and the intervalley transition is completed. The whole process is presented in Fig. 12, where the green curve represents the valley index K evolution during the entire transition.
Besides inter-spin and valley transitions we can simultaneously obtain spin and valley manipulation leading to inter spin-valley transitions or simply the spin-valley SWAP. Similarly here transitions are resonant and we need to tune the modulation frequency ω to the energy spacing between |−1, ↑ and |1, ↓ (or |−1, ↓ and |1, ↑ ). During voltage oscillations both the Rashba SOI mediated spin transitions and the intervalley coupling modulation effects are enabled, thus allowing for simultaneous spin and valley flipping. Both effect are needed: turning off the Rashba SOI, by setting E z (x, y) = 0, turns off the spin-valley SWAP. Simulation results are presented in Fig. 11 with the resonance frequency ω 0 = 0.244 meV and the driving voltage V ac = 100 mV. We observe here simultaneous spin (blue curve) and valley index (violet curve) flips. These transitions are obviously of the Rabi oscillations type. If we diverge from the resonance, the maximum (minimum) value of the valley (spin) index falls down rapidly, entering the region of incomplete transitions. In Fig. 11, green and orange (red and yellow) curves pair present spin and valley index courses for a driving frequency ω = 1.003 ω 0 (ω = 1.006 ω 0 ) beyond the resonance.
In opposite to spin-valley SWAP, the intervalley transitions are not mediated by the spin-orbit coupling, and are unaffected even if we eliminate the electric field in the simulations, simply by setting E z (x, y) = 0 in (Eq. 3). Moreover, the transition time depends on the amplitude of the intervalley coupling modulation. In Fig. 12(a) are presented transitions between both valleys, starting from K valley with the index K = 1. The transition period T deceases as the voltage modulation amplitude V ac increases. Indeed, for presented in Fig. 12 amplitude V ac ranges (50-150 mV), oscillations frequency Ω turns out to be approximately proportional to V ac , and thus to the intervalley coupling modulation amplitude Λ 1 (we assume here that for such a small voltage modulation range the intervalley coupling 2Λ responses linearly-cf. Fig. 7). This is typical for the Rabi oscillations, where near the resonance the Rabi frequency Ω = (ω − ω 0 ) 2 + (Λ 1 / ) 2 depends linearly on the driving amplitude Λ 1 of the intervalley coupling oscillations.
In case of resonance Ω = Λ 1 / . If we calculate the minimum value reached by the K index for out-of-resonance transitions, we obtain resonance curves presented in Fig. 12(b). The full width at half maximum (FWHM) parameter characterizing resonance curves in case of the Rabi oscillations corresponds to the transition duration. The resonance curve (for driving Λ1 2 cos(ωt)) has the form (Λ1/ ) 2 (ω−ω0) 2 +(Λ1/ ) 2 , which gives FWHM equal 2Λ 1 . This agrees with our calculations. E.g. for the green curve, i.e. V ac = 100 mV, transition time (period) T = 2π/Ω is 680 ps, which corresponds to Λ 1 = h/T 0.006 meV and agrees with FWHM equal 0.012 meV. In comparison, for the violet curve FWHM is just over two times wider than for orange one.  Simply, by tuning the modulation frequency we can select and switch-on desired transition. If we apply a proper magnetic field value (we assume B z = 1 T), the Zeeman splitting together with the spin-orbit induced splitting result in different frequencies among transitions (and allows to separately address each of them). If we now sweep the driving frequency over a range covering all the six transitions, assuming that the system can be initially in four different basis states we observe in Fig. 13(b) all of transitions at their own frequencies. They are highlighted by colors corresponding to colors of the arrows from the scheme in Fig. 13(a). Interestingly we also observe some minor fractional resonances for lower frequencies ω/2. Fortunately, they do not overlap with other peaks and transitions are not disturbed by each other. The driving amplitude applied in presented simulations is V ac = 100 mV, with an exception for SWAPs where V ac = 150 mV.
Let us now translate obtained transitions to the language of qubit operations. Starting from the blue transition from Fig. 13, for pumping at 2.063 meV, we obtain a spin-flip only if K = −1, i.e. K valley is occupied. In case of the valley index K = 1, we would not observe any operation on spin for such a driving frequency. This means that we get the spin NOT quantum operation controlled by the valley qubit. We denote it simply by C K N OT s . For the violet transition (2.333 meV) we get the opposite CNOT operation with the spin qubit flipped if K = 1, denoted byC K N OT s . On the other hand, for the green transition we get complementary operation with the valley index being rotated only if the spin is oriented down. In this case acquiring spin-controlled NOT quantum operation on the valley qubit, analogously denoted as C s N OT K . The purple transition is performed for the opposite, spin-up, thus denoted byC s N OT K . Let us note that the CNOT gates are essential in creating the universal set of quantum gates. Any (multi-qubit) quantum operation can be approximated by a sequence of gates from a set consisting CNOT gate and some singlequbit operation 58 , e.g. the R π/8 gate. If we stop our transitions earlier, we can get various rotation gates. In particular, limiting the operation time in Fig. 14 to 1/4 of the full valley (spin) flip, i.e. 85 ps (230 ps), we realize the R π/8 rotation acting on the valley (or spin) qubit.
Beside the both CNOT operations with spin or valley serving as the control qubit, while the other one being the target qubit, we can create previously mentioned SWAPs. By taking the red transition we get spin and valley states swapped, i.e. | − 1, ↑ ↔ |1, ↓ . The complementary operation, induced by the yellow transition, interchanges the two remaining states: |−1, ↓ ↔ |1, ↑ . We denote it by cSWAP . All the mentioned two-qubit operations are represented in the spin-valley two-qubit subspace of |K, s states by 4 × 4 unitary matrices. Their explicit form can be found in the appendix.

C. Single-qubit gates
Two-qubit gates are easy to implement here, because the both qubits are specified on two degrees of freedom of the same particle, thus defined in the same localization. Therefore, coupling between them emerges naturally and two-qubit operations require a single transition between one of the four electron basis states. On the other hand, to obtain a single-qubit gate, acting on a given qubit within such a subspace must be done independently from the other qubit state. It turns out that joining two opposite CNOTs makes the operation on the target qubit independent from the control one. If we perform simultaneously both valley-controlled spin NOTs, i.e. C K N OT s andC K N OT s we arrive at single spin-NOT quantum gate, denoted as N OT s , independent from the valley degree. Similarly, for simultaneous C s N OT K and C s N OT K we get valley-NOT, N OT K quantum operation.
Indeed, to obtain correct operations on spin or valley separately, we need to pump two transitions at the same time. Luckily, it turns out that such twofold transitions are possible, and to do that we need to simultaneously induce oscillations in perpendicular directions, i.e. V 1 (t) = V dc + V x ac sin(ω x t) and V 2 (t) = V dc + V y ac sin(ω y t), by feeding both G 1 and G 2 gates. In Fig. 14  The both pumping frequencies ω x and ω y (blue and violet transitions in Fig. 14(bottom)) are very slightly different from these for single separate transitions, e.g. for spin-NOT single-qubit gate they change from (2.06, 2.33) to (2.07, 2.35) meV for (ω x , ω y ) respectively. Whereas the voltage oscillation amplitudes pair (V x ac , V y ac ) should be selected in a way that the both transitions in the pair lasts the same time. The ratio between them should be properly tuned, e.g. for valley-NOT single-qubit gate (green and purple transitions in Fig. 14(top)) V x ac /V y ac = 100mV/62mV for (ω x , ω y ) = (1.82, 2.58) meV respectively.
We see an intriguing feature that local-defined on single electron-two-qubit gates are easier to implement than single-qubit. However, the same single-qubit operations can be obtained a bit easier, simply by performing appropriate CNOTs one by one. Unfortunately, in this simpler approach the operation time is twice as long as for simultaneous twofold transitions. Applying series of operations representing the particular operations.

D. Gate fidelities and qubits readout
In the course of transitions from Figs. 14 or 12 we can notice a minor oscillation structure of frequency ω related to a single cycle of pumping induced by the voltages oscillations. This can be viewed as a reference frame rotation with ω frequency in the standard RWA approximation 59 . However, it should be emphasized that our numerical calculations are strict. These small oscillations affect the qubit operations fidelity. Fortunately, their amplitude decreases as we get closer to the basis states (i.e. poles on the Bloch sphere). Moreover, we can reduce them arbitrarily by decreasing the amplitude V ac of the voltage oscillations. This is at the expense of increasing the number of pumping cycles, and thus increases the operation time.
We have performed simulations of the operations for decreasing voltage amplitudes, with simultaneously increasing the gates time. In Fig. 15 there is presented fidelity of various operations as a function of their gate time. One should find an appropriate trade-off between high gate fidelity and low operation time. For example, to obtain an error of the order of 1% (99% fidelity) we find valley-related gates (C s N OT K and N OT K ) duration of about half a nanosecond, and spin-gates about 1 ns. The lowest fidelity is for SWAP operations with a 2 ns optimum for 99%-fidelity. The duration of operations should be much shorter than the coherence time. Estimated 6 coherence times of electron valley and spin degrees of freedom in MoS 2 monolayers, related to the hyperfine interaction decoherence from Mo nuclear spins, are of 100 ns. This is about two orders of magnitude longer than our operation times. The coherence timescale may be slightly overestimated 60 . However, it scales with the system size as √ N I , with N I the number of nuclear spins covered by the wavefunction density. Thus, coherence time scales up linearly with the system length, and can be extended by increasing the confinement width. Each full qubit implementation has to comprise initialization and readout. To do this we can utilize the valley-and spin-Pauli blockade, so far observed in carbon nanotubes 61,62 . The blockade utilizes selection rules, which block electron transport between an adjacent dot with the same valley and spin state. Let us add to the setup a nearby auxiliary dot with an electron in the ground state, i.e. |K 0 , s 0 = | − 1, ↑ . Assuming that valley and spin are conserved during tunneling, the electron carrying our qubits cannot tunnel to the nearby dot if the electron confined there occupies the same spin and valley state 5,63 . The electron is blocked in the same state as its neighbor: | − 1, ↑ (1, 1)| − 1, ↑ and both qubits are initialized. However, when we perform operation on valley or spin qubit, the blockade is lifted and the electron can freely enter the nearby dot: |−1, ↑ (1, 1)|−1, ↑ → (0, 2)|K 1 , s 1 |−1, ↑ with K 1 = 1 or s 1 = ↓. In this way, by extending the system with adjacent reference electron, we can perform both-spin or valley qubits readout.

VI. SUMMARY
The emerging branch of electronics utilizing the valley degree of freedom, called valleytronics in analogy to spintronics, introduces new intriguing methods for defining qubits. Nanodevices with gated monolayer QDs currently become more reliably fabricated. To advantage this, we investigate the possibility of realization of the spin-valley two-qubit system defined on a single electron, that is confined in a QD, controlled by voltages applied to the device local gates. The proposed nanodevice is modeled after structures that were experimentally realized.
A realistically calculated QD confinement potential and electric field via the Poisson equation together with the exact form of the Rashba coupling within the tightbinding monolayer model, leads to reliable modeling of both the intervalley coupling and the Rashba SOI. We solve the time-dependent Schrödinger equation with such variable confinement, and track the transitions by calculating actual values of the spin and valley index. We also analyze the edge influence on the intervalley coupling with the increasing flake size and confinement depth, concluding that proposed nanodevice will work even if we enlarge the flake and move its edges away.
As a result of the performed simulations, we show feasibility of electrically controlling both the electron spin and valley degrees of freedom, simultaneously. By applying an appropriate magnetic field we get such spin and valley Zeeman splittings, that all of the six possible transitions within the spin-valley subspace can be separately addressed. These transitions are interpreted as a variety of two-qubit gates (i.e. CNOTs and SWAPs), and properly combined, they give single-qubit NOT gates.
Encoding two qubits locally on two degrees of freedom of a single electron reverse difficulty in such a way that two-qubit gates are easer to implement than a single qubit. The latter, however, can be achieved in one go or as two consecutive transitions. By examining the exact course of transitions, we can also estimate fidelity of the implemented gates.
Finally, we remark that to implement fully scalable system, we also need to control interaction between valley (and spin) indexes of nearby electrons in the register. This requires adding additional electrons to the system and research interactions among them.

ACKNOWLEDGMENTS
Author would like to thank Grzegorz Skowron and Pawe l Potasz for invaluable discussions. This work has been supported by National Science Centre, under Grant No. 2016/20/S/ST3/00141. This research was supported in part by PL-Grid Infrastructure.
As expected, resulting specific numeric value e|E z γ αβ of each block (α, β) is proportional to the external electric field E z . Moreover, it is multiplied by the matrix with phase factors e ±ıη of the form equivalent to e y ij σ x − e x ij σ y in the H R Hamiltonian from (Eg. 3). The phase depends on the hopping direction, however in our calculations it is undetermined and resulting η was disordered.
In this way we obtain the parameter matrix γ αβ R = |e|E z γ αβ , where γ R = |e|E z γ, γ = Finally, we take to the model an approximate-average value of the γ as written down in (2), remembering, however, of slightly different values between valleys. The obtained value corresponds to the Rashba coupling amplitude from [2], where λ R = 3.3 × 10 −4 E z (eV nm), for E z expressed in the V/nm units. E.g. for 1 V/nm and k of order 4π 3a we get λ R k ∼ 4.3 meV (a = 0.319 nm). While in our calculations, for a similar electric field, the γ is reaching comparable values of a few meVs.

Appendix B: energy spectrum in magnetic field
Applying an external magnetic field introduces a splitting of levels with opposite spin (and also valley). If we now gradually increase the magnetic field we obtain typical energy levels structure in a quantum dot with a magnetic field, presented in Fig 16. Similar structures are shown in [70] (or [2]), calculated within the k.p model for 50-nm-size (or 40 nm) MoS 2 QD. However here we have more than 10 times smaller dot, thus to obtain similar orbital effects relatively to the magnetic scale eB (size of the Landau ground state), equal ∼ 25.66 nm for 1 T, we need to increase B field more than 100 times. To keep the results in Fig. 16 comparable, we omit the Zeeman term here. Presented results are calculated for N = 15 and V dc = −1500 mV, forming the QD confinement.
If we further increase the magnetic field, the energy levels will start to attract to each other and form characteristic Landau levels. Calculations for the whole range of artificially high magnetic fields 0 < B < B 0 , B 0 = 9.4 × 10 4 T, shows that Landau levels posses complicated self-similar structure, called the Hofstadter butterfly 47,71-73 . It is presented in Fig. 17. Complex regularities are also manifested in colors that represent the QD occupation. Same here we skip the Zeeman energy term.
Appendix C: two-qubit gate matrices Here we present the explicit forms of matrices representing the all two-qubit operations obtained in the simulations. They act in the 4-dimensional Hilbert space of the two-qubit spin-valley states |K, s = |K ⊗ |s . The first pair constitutes CNOT operations where the valley is the control qubit and the spin is the target one: Note that performing jointly both SWAP operations is equivalent to the NOT operation on the both qubits σ x ⊗ σ x , i.e. SWAP · cSWAP = cSWAP · SWAP = N OT K ⊗ N OT s .