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

TL;DR

This paper proposes a realistic gate-defined MoS$_2$ monolayer quantum dot device that enables electrical control of spin and valley states for quantum computing, demonstrating potential for universal qubit gates.

Contribution

It introduces a detailed model and numerical simulations of a MoS$_2$ quantum dot device capable of implementing single- and two-qubit gates via spin-valley control.

Findings

Electrical control enables single- and two-qubit gates.

Two-qubit gates are easier to realize than single-qubit gates.

The device design is feasible with realistic potentials and electric fields.

Abstract

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.

Equations

$$H=$$

$$+$$

$$\gamma =\left(\begin{array}{ccc}0.1& 5& 17\\ 5& 4& 1\\ 17& 1& 4\end{array}\right)\times 10^{-3}\mathrm{nm}.$$

$$H_R=\sum _{⟨ij⟩\sigma {\sigma }^{\prime }\alpha \beta }ı{\gamma }_R^{\alpha \beta }({\hat{\mathbf{e}}}_{ij}\times \hat{\mathbf{z}})\cdot {\mathbf{s}}_{\sigma {\sigma }^{\prime }}{c}_{i\alpha \sigma }^{†}{c}_{j\beta {\sigma }^{\prime }},$$

$$H_Z=\sum _{i\sigma {\sigma }^{\prime }\alpha \beta }{\gamma }_Z\mathbf{B}\cdot {\mathbf{s}}_{\sigma {\sigma }^{\prime }}{\delta }_{\alpha \beta }{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{i\alpha {\sigma }^{\prime }},$$

$${\theta }_B=\frac{e}{\hslash }{\int }_{{\mathbf{r}}_i}^{{\mathbf{r}}_j}\mathbf{A}(\mathbf{r})\cdot d\mathbf{r}.$$

$${\theta }_B=\frac{e}{\hslash }\left(A_y({\mathbf{r}}_i)+A_y({\mathbf{r}}_j)\right)(y_j-y_i)/2.$$

$$\mathbf{\Psi }(\mathbf{r},t)=\sum _n{c}_n(t){\mathbf{\psi }}_n(\mathbf{r})e^{-\frac{ı}{\hslash }{E}_{n}t},$$

$$ı\hslash \frac{\partial }{\partial t}\mathbf{\Psi }(\mathbf{r},t)=H(\mathbf{r},t)\mathbf{\Psi }(\mathbf{r},t),$$

$$H(\mathbf{r},t)=H(\mathbf{r})+\delta \phi (\mathbf{r},t)+\delta H_R(\mathbf{r},t).$$

$${\dot{c}}_m(t)=-\frac{ı}{\hslash }\sum _n{c}_n(t){\delta }_{mn}(t)e^{\frac{ı}{\hslash }(E_m-E_n)t}.$$

$$\tilde{\mathbf{\Psi }}(\mathbf{k},t)={\int }_F{d}^{2}r\mathbf{\Psi }(\mathbf{r},t)e^{-ı\mathbf{kr}},$$

$$\mathcal{K}(t)=\frac{3a}{4\pi }{\int }_{{\tilde{F}}_{1/3}}{d}^{2}k\tilde{\rho }(\mathbf{k},t)k_x$$

$$s(t)={\int }_F{d}^{2}r{\mathbf{\Psi }}^{†}(\mathbf{r},t){\sigma }_z\otimes 1_3\mathbf{\Psi }(\mathbf{r},t),$$

$$\left(\begin{array}{cc}{H}_{6\times 6}& H_{6\times 16}\\ H_{6\times 16}^{†}& H_{16\times 16}\end{array}\right)\left(\begin{array}{c}{\psi }_6\\ {\psi }_{16}\end{array}\right)=\mathcal{E}\left(\begin{array}{c}{\psi }_6\\ {\psi }_{16}\end{array}\right).$$

$$H_{6\times 6}^{\prime }=H_{6\times 6}+H_{6\times 16}{\left(\mathcal{E}-H_{16\times 16}\right)}^{-1}{H}_{6\times 16}^{†},$$

$$a_0{1}_2+|e|E_z{\gamma }_{\alpha ,\beta }\left(\begin{array}{cc}0& e^{-i\eta }\\ e^{i\eta }& 0\end{array}\right).$$

$${\gamma }_R=|e|E_z\gamma ,\gamma =\left(\begin{array}{ccc}0.09& 5.62& 17.4\\ 5.62& 3.59& 1.26\\ 17.4& 1.26& 4.97\end{array}\right)\times 10^{-3}\mathrm{nm}.$$

$$\gamma =\left(\begin{array}{ccc}0.08& 4.58& 16.9\\ 4.58& 4.65& 0.71\\ 16.9& 0.71& 3.10\end{array}\right)\times 10^{-3}\mathrm{nm},$$

$$C_{\mathcal{K}}NO{T}_s=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),{\bar{C}}_{\mathcal{K}}NO{T}_s=\left(\begin{array}{cccc}0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right).$$

$$C_{s}NO{T}_{\mathcal{K}}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\end{array}\right),{\bar{C}}_{s}NO{T}_{\mathcal{K}}=\left(\begin{array}{cccc}0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& 0\\ 0& 0& 0& 1\end{array}\right).$$

$$SWAP=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right),cSWAP=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right).$$

Full text

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. MoS2, seem to be better candidates than graphene because of their wide band gaps and strong electrically induced spin-orbit coupling of the Rashba typeWang et al. (2012); Kormányos et al. (2014a). By considering the valley degree of freedom of an electron together with its spin we extend our ability to define a qubit into: spin, valleyPawłowski et al. (2018) and hybrid spin-valley qubitSzéchenyi et al. (2018). However, the most interesting is the definition based on spin and valley of a single electron as a two-qubit systemRohling and Burkard (2012); Wu et al. (2016).

The area of application of monolayer materials for construction of electronic nanodevices is currently under strong development Papadopoulos et al. (2019); Shang et al. (2019); Paul et al. (2019); Wang et al. (2012); Pisoni et al. (2019); Ghiasi et al. (2019); Huang et al. (2018). Methods for building devices based on gated TMDC monolayers or nanotubes become increasingly advancedKim et al. (2019a); Sahoo et al. (2019); Wang et al. (2018); Kim et al. (2019b); Reinhardt et al. (2019), 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 MoS2-nanoribbon-QD measured by a single electron transportKotekar-Patil et al. (2019); Zhang et al. (2017), or tunable TMDC spintronic devices, where spin or valley polarized currents emerge in TMDC monolayer proximitized by nearby ferromagneticCortés et al. (2019); Zollner and Fabian (2019); Yuan et al. (2018); Li et al. (2014); Majidi and Asgari (2014). Transistor structures with TMDC monolayer forming active area in tunnel FETs are being developedIlatikhameneh et al. (2015), also with vertical TMDCs heterostructures Georgiou et al. (2012); Wu et al. (2019). The more intriguing lateral, in-plane TMDCs heterojunctions are also constructed, enabling interesting 1D physics at interfacesÁvalos-Ovando et al. (2019), or leading to improved FETs switching characteristics Choukroun et al. (2018); Marian et al. (2017); Iannaccone et al. (2018).

Inspired by this, we examined the possibility of realization of a nanodevice based on a MoS2 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 capabilities. 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.

II.1 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 SiO2. Then we place two electrodes which serve as a source ($S$) and a drain ($D$). Directly on them we deposit a MoS2-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 bandgapLaturia et al. (2018), 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\times 20$ nm. The gate layout presented here is quite similar to the one proposed by us recentlyPawłowski et al. (2018), 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 $\phi(\mathbf{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 devicePawłowski et al. (2018, 2016). Resulting potential in the area between SiO2 and hBN layers, where the flake is sandwiched, is presented in Fig. 2.

We deplete the electron gas until a single electron remains in the formed dot confinement potential.

II.2 Monolayer model

The monolayer flake is made of molybdenum disulfide.

MoS2 monolayers are successfully described by several tight-binding (TB) models, with different numbers of orbitals used, including nearest or next-nearest neighbors. SevenRostami et al. (2013) or elevenRidolfi et al. (2015); Cappelluti et al. (2013); Fang et al. (2015) Mo and S orbitals construct TB basis to reproduce low-energy physics in the entire Brillouin zone, also near the $\mathrm{\Gamma}$-point. Although the simpler three-band (including three Mo orbitals) TB modelLiu et al. (2013) fails around the $\mathrm{\Gamma}$-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$ orbitalsKadantsev and Hawrylak (2012). 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 MoS2 structure (see Fig. 3). For calculating Rashba coupling we will utilize 11-band model with three $p$ orbitals for each S atom in dimer.Ridolfi et al. (2015)

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 hoppingsLiu et al. (2013):

[TABLE]

The potential energy of the electrostatic confinement at the $i$-th lattice site: $\varphi_{i}=-|e|\phi(x_{i},y_{i})$ together with the on-site energies $\epsilon_{\alpha}$ enter the on-diagonal matrix elements ($\alpha$ numbers the orbitals).

The off-diagonal electron hopping element from the $\beta$ Mo orbital localized in the $j$-th lattice site to the $\alpha$ orbital localized in the $i$-th site is denoted by $t_{\alpha\beta}\equiv t_{\alpha\beta}\!\left(R_{k\left(i,j\right)}\right)$. It depends on the hopping direction (between $\langle ij\rangle$ 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_{\alpha\beta}(R_{k})=t_{\beta\alpha}(-R_{k})$. Their explicit forms, together with the on-site energies $\epsilon_{\alpha}$, can be found in [Pawłowski et al., 2018; Liu et al., 2013].

II.3 Rashba coupling

Electric field perpendicular to the monolayer surface breaks the reflection $\sigma_{h}$ symmetry and modifies the on-site energies of atoms in three MoS2 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 Cortés et al. (2019).

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 $\lambda_{\alpha\beta}$ in Eq. II.2) including also $p$ orbitals of the sulfur top and bottom sublayers to which we apply on-site potentials $V_{\mathrm{t}}$ and $V_{\mathrm{b}}$Ochoa and Roldán (2013); Rostami et al. (2013); Petersen and Hedegård (2000); Konschuh et al. (2010). The difference between them results from external electric field: $V_{\mathrm{b}}(x,y)=E_{z}(x,y,0)d/2$, $V_{\mathrm{t}}(x,y)=-V_{\mathrm{b}}(x,y)$. While $d=0.32$ nm is the monolayer thickness (sulfur sublayers distance).

Given such extended tigh-binding model (we take 11-orbital model of Ridolfi et al. [Ridolfi et al., 2015]), we perform downfolding, using the Löwdin partitioning technique, to our Mo-orbitals model and obtain the Rashba coupling $\gamma_{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: $\gamma^{\alpha\beta}_{R}(x,y)=|e|E_{z}(x,y)\,\gamma_{\alpha\beta}$, with the explicit form

[TABLE]

The tight-binding Rashba HamiltonianKlinovaja and Loss (2013a, b); Konschuh et al. (2010); Ezawa (2014); Kane and Mele (2005) (in Eq. II.2) with characteristic spin- and orientation-dependent hopping between two nearest neighbor bonds is:

[TABLE]

with the Pauli-matrices vector $\mathbf{s}=(\sigma_{x},\sigma_{y},\sigma_{z})$, and the (unit) versor $\hat{\mathbf{e}}_{ij}\equiv(e^{x}_{ij},e^{y}_{ij})$ pointing along the bond connecting sites $i$ and $j$. An obvious property $\hat{\mathbf{e}}_{ji}=-\hat{\mathbf{e}}_{ij}$ ensures $H_{R}$ hermiticity. The expanded hopping expression is $(\hat{\mathbf{e}}_{ij}\!\times\hat{\mathbf{z}})\!\cdot\mathbf{s}_{\sigma\sigma^{\prime}}=e^{y}_{ij}s^{x}_{\sigma\sigma^{\prime}}-e^{x}_{ij}s^{y}_{\sigma\sigma^{\prime}}$, where $s^{x}=\sigma_{x}$, $s^{y}=\sigma_{y}$, and $\gamma_{R}=|e|E_{z}\gamma$ with the above defined $\gamma$.

II.4 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. II.2):

[TABLE]

with a magnetic field $\mathbf{B}$. For $\gamma_{Z}=\frac{g_{e}\mu_{B}}{2}$ we arrive at the standard Zeeman energy $\frac{g\mu_{B}}{2}\mathbf{s}\cdot\mathbf{B}$.

To address also orbital effects related to magnetic field we apply the so-called Peierls substitutionHofstadter (1976). We multiply the hopping matrix, by the additional factor $t_{ij}\rightarrow\tilde{t}_{ij}=t_{ij}\exp\left(\imath\theta_{B}\right)$ in the Hamiltonian (II.2). Now the vector potential enters Eq. (II.2) via the Peierls phase $\theta_{B}$, calculated as the path integral between neighbor nodes:

[TABLE]

$\mathbf{A}$ is a vector potential induced by the $\mathbf{B}$ field. We use the Landau gauge, with the vector potential $\mathbf{A}(\mathbf{r})=[0,\int dxB_{z}(x,y),0]^{T}$ for the perpendicular magnetic field $\mathbf{B}(\mathbf{r})=[0,0,B_{z}(\mathbf{r})]^{T}$. This leads to the phase:

[TABLE]

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 (II.2): $H(\mathbf{r})\bm{\psi}_{m}(\mathbf{r})=E_{m}\bm{\psi}_{m}(\mathbf{r})$, and obtain $M$ eigenstates. For our hexagonal flake of size $N_{B}=25$ we have $1801\,\mathrm{sites}\,\times 6$ giving $M=10806$ eigenstates. For the Hamiltonian matrix eigenproblem we utilize the fast and efficient FEAST routinePolizzi (2009). The obtained eigenstates are represented by 6-dimensional vectors $\bm{\psi}_{m}(\mathbf{r})=(\psi^{\sigma\alpha}_{m}(\mathbf{r}))^{\intercal}$, with $\sigma=1,2$ and $\alpha=1,2,3$. They belong to te state space ${\cal H}_{2}^{\mathrm{spin}}\otimes{\cal H}_{3}^{\mathrm{oribital}}$, with the spin and the 3-dimensional Mo-orbitals space. To identify them, at first we need an electron density calculated as $\rho(\mathbf{r})=|\bm{\psi}_{m}(\mathbf{r})|^{2}$ to determine if given state is localized at the flake edge forming the so-called edge state, or is confined within the quantum-dot. 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 $\bm{\psi}_{m}(\mathbf{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 $\bm{\psi}_{n}$:

[TABLE]

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:

[TABLE]

with the time-dependent Hamiltonian being a sum of the stationary part (Eq. II.2) and a time-dependent contribution to both the potential energy and the Rashba coupling:

[TABLE]

The full time-dependent potential energy $\varphi(\mathbf{r},t)=\varphi(\mathbf{r})+\delta\varphi(\mathbf{r},t)$ contains variable part $\delta\varphi(\mathbf{r},t)$, generated by modulation of the gate voltages. Whole is calculated as $\varphi(\mathbf{r},t)=-|e|\phi(\mathbf{r},t)$, with the potential $\phi(\mathbf{r},t)$ obtained by solving the Poisson equation for the variable density $\rho(\mathbf{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 self-consistent way. Similarly, the time-dependent part of the Rashba coupling $\delta H_{R}(\mathbf{r},t)$ is induced by the variable part of the electric field: $E_{z}(\mathbf{r},t)=E(\mathbf{r})+\delta E_{z}(\mathbf{r},t)$. That is, $E(\mathbf{r})$ induces (3), while $\delta E_{z}(\mathbf{r},t)$ enters $\delta H_{R}(\mathbf{r},t)$ in (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:

[TABLE]

The actual matrix elements $\delta_{mn}(t)=\langle\bm{\psi}_{m}|\delta\varphi(\mathbf{r},t)+\delta H_{R}(\mathbf{r},t)|\bm{\psi}_{n}\rangle$ 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 $\bm{\Psi}(\mathbf{r},t)$ we calculate the Fourier transform:

[TABLE]

on the flake surface area $F$, with 2D-wave vector $\mathbf{k}\equiv(k_{x},k_{y})$. The Fourier transform naturally has periodic structure in the reciprocal space, therefore we can limit the $k$-area to $\tilde{F}$: $k_{x,y}\in\left[-\frac{2\pi}{a},\frac{2\pi}{a}\right]$, encompassing the (first) Brillouin Zone (BZ). Knowing $\tilde{\bm{\Psi}}(\mathbf{k},t)$ we can calculate density in the reciprocal space expressed as: $\tilde{\rho}(\mathbf{k},t)=|\tilde{\bm{\Psi}}(\mathbf{k},t)|^{2}$. The $k$-density calculated for the $|K\!\!\downarrow\rangle$ state from Fig. 5 is presented in Fig. 10([math]). We also mark the BZ along with points of high symmetry: $\Gamma$ 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: $\Gamma=(0,0)$, one of $K=\frac{\pi}{a}(\frac{4}{3},0)$ and one of $M=\frac{\pi}{a}(1,\frac{1}{\sqrt{3}})$, with the lattice constant $a=0.319$ nm. We can clearly see in Fig. 10([math]) that density peaks are localized in the neighborhood of the $K$ points (while not next to $K^{\prime}$), confirming that in the $|K\!\!\downarrow\rangle$ state exactly $K$ valley is occupied.

Now, the valley index $\mathcal{K}$ is calculated as:

[TABLE]

on the reciprocal space area $\tilde{F}_{1/3}$ defined as two opposite $\pi/3$ sectors (within $\tilde{F}$ area) encompassing exactly one $K$ point and one $K^{\prime}$ point, i.e. $|\vartheta|\leq\frac{\pi}{6}\cup|\vartheta|\geq\frac{5\pi}{6}$, with azimuthal angle $\vartheta=\mathrm{atan2}(k_{x},k_{y})$. Because $K$($K^{\prime}$) point in $\tilde{F}_{1/3}$ has coordinates $1$($-1$)$\times(\frac{4\pi}{3a},0)$, the valley index $\mathcal{K}\in\left[-1,1\right]$. $\mathcal{K}=1$ represents the $K$ valley, whereas $\mathcal{K}=-1$ the $K^{\prime}$ valley. For example the state $|K\!\!\downarrow\rangle$ is thus represented by $|\mathcal{K},s\rangle=|1,\downarrow\rangle$ ket.

The electron spin value $s$ is calculated as the expectation value of the Pauli $z$-matrix $s^{z}=\sigma_{z}$:

[TABLE]

also integrated on the flake area $F$ for the actual wave vector $\bm{\Psi}(\mathbf{r},t)$. Operation $\otimes\,\mathbf{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 (II.2) for the entire flake lattice, forms a ladder in Fig. 5 representing subsequent $M=10806$ eigenstates $\bm{\psi}_{m}(\mathbf{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 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,\uparrow\rangle,|1,\downarrow\rangle\}$ and $\{|\!-\!1,\downarrow\rangle,|1,\uparrow\rangle\}$, 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 $\Gamma$ 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 ${\cal H}_{2}^{\mathrm{valley}}\otimes{\cal H}_{2}^{\mathrm{spin}}$ of spin-valley two-qubit states $|\mathcal{K},s\rangle$.

IV.1 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: $\mathcal{E}_{K^{\prime}\uparrow}=\frac{1}{2}(g_{s}-g_{v})\mu_{B}B$ (split down) and $\mathcal{E}_{K\downarrow}=\frac{1}{2}(-g_{s}+g_{v})\mu_{B}B$ (split up) is $242$ $\mu$eV, which is in agreement with differences $238$ $\mu$eV between opposite spin states (with different valleys) in the conduction band minimum for a larger dot [Dias et al., 2016]. The difference of $242$ $\mu$eV leads to $g_{v}-g_{s}\simeq 4.18$. For the second pair $\mathcal{E}_{K^{\prime}\downarrow}=\frac{1}{2}(-g_{s}-g_{v})\mu_{B}B$ (split down) and $\mathcal{E}_{K\uparrow}=\frac{1}{2}(g_{s}+g_{v})\mu_{B}B$ (split up) we have $512$ $\mu$eV, thus $g_{v}+g_{s}\simeq 8.84$. Therefore, obtained effective valley and spin g-factors are: $g_{v}=6.51$, similar to DFT calculations giving the value 7.14 [Kormányos et al., 2014a, and errata: Kormányos et al., 2014b]. While the spin splitting factor $g_{s}=2.33$, with agreement with the DFT calculations [Kormányos et al., 2014a] and experimental result [Marinov et al., 2017], 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^{\prime}$ 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.

IV.2 Intervalley coupling

Let’s now examine the intervalley coupling strength and its origin. The doubled intervalley coupling $2\Lambda$ can be calculated from the difference between the ground state and the 1st excited state in CB with no spin-orbit interactionLiu et al. (2014). Further, we assume that the intervalley coupling is: $\Lambda(t)=\Lambda_{0}+\frac{\Lambda_{1}}{2}\cos(\omega t)$, with a modulation amplitude $\Lambda_{1}$.

The structure of electronic states confined within a nanoflake and mostly presence of edge states, that cross the gap, depends on the flake edge typeBrey and Fertig (2006); Zarenia et al. (2011); Szafran et al. (2018a). 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 structuresSzafran et al. (2018b). Same edge-dependent valley mixing was proven for MoS2 nanoribbonsRostami et al. (2016). 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\Lambda$ 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, by dark-blue bullets. While there is significant coupling for very small flakes with $N_{B}\sim 10$, it suddenly decreases with the flake size, reaching two orders smaller value for $N_{B}\sim 15$. Then the coupling slightly increase, but for larger flakes with $N_{B}>25$ it generally does not exceed several $\mu$eV. However if we add the confinement potential by applying $V_{dc}$, the intervalley coupling increases with the potential depth, which is clearly visible for $N_{B}>25$. Subsequent voltage values, marked by the brighter bullets, are $V_{dc}$ = $-500,-1000,-1500,-2000,-2500$ mV. It turns out that the approximate relation between the coupling and the applied voltage in this range of the flake sizes is $2\Lambda\sim\sqrt{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.

V Spin-valley two qubit system

Let’s now switch to the time-dependent calculations and examine the process of inducing transitions within the defined spin-valley subspace. By applying additional oscillating voltage to the single control gate $G_{1}$: $V_{1}(t)=V_{dc}+V_{ac}\sin(\omega t)$, $V_{ac}=-100$ mV together with $V_{2}=V_{3}=V_{4}=V_{dc}=-2500$ mV, we modulate the 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=\frac{\pi}{2\omega}$ and $t=\frac{3\pi}{2\omega}$, 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).

V.1 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 $\frac{d}{dt}\langle x\rangle$ is presented in Fig. 9 as the red curve, while the electron position as the orange one. Together with the present perpendicular electric field $\langle E_{z}\rangle$ 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 $\langle E_{z}\rangle$ 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 oscillations of the electron packet size $\sigma$ 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 statesPawłowski et al. (2018).

Subsequent stages of a transition between the different valleys in reciprocal space are presented in Fig. 10.

Initially, the electron density in the reciprocal space $\tilde{\rho}(k_{x},k_{y})$ is localized in the $K$ point vicinity within the BZ, as showed in Fig. 10([math]). The voltage pumping process with the resonant frequency $\omega$, tuned to the energy spacing $\hbar\omega$ between the states $|1,\uparrow\rangle$ and $|\!-\!1,\uparrow\rangle$ (or $|1,\downarrow\rangle$ and $|\!-\!1,\downarrow\rangle$), leads to a gradual change of occupation to $K^{\prime}$ valley, with density flow between valleys visible in the subsequent stages—Fig. 10($\frac{\pi}{4}$) and ($\frac{\pi}{2}$). After time $t=345$ ps the electron occupies the $K^{\prime}$ valley entirely (Fig. 10($\pi$)) and the intervalley transition is completed. The whole process is presented in Fig. 12, where the green curve represents the valley index $\mathcal{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 $\omega$ to the energy spacing between $|\!-\!1,\uparrow\rangle$ and $|1,\downarrow\rangle$ (or $|\!-\!1,\downarrow\rangle$ and $|1,\uparrow\rangle$). 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 $\hbar\omega_{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 $\omega=1.003\,\omega_{0}$ ($\omega=1.006\,\omega_{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 $\mathcal{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 $\Omega$ turns out to be approximately proportional to $V_{ac}$, and thus to the intervalley coupling modulation amplitude $\Lambda_{1}$ (we assume here that for such a small voltage modulation range the intervalley coupling $2\Lambda$ responses linearly—cf. Fig. 7). This is typical for the Rabi oscillations, where near the resonance the Rabi frequency $\Omega=\sqrt{(\omega-\omega_{0})^{2}+(\Lambda_{1}/\hbar)^{2}}$ depends linearly on the driving amplitude $\Lambda_{1}$ of the intervalley coupling oscillations. In case of resonance $\Omega=\Lambda_{1}/\hbar$.

If we calculate the minimum value reached by the $\mathcal{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 $\frac{\Lambda_{1}}{2}\cos(\omega t)$) has the form $\frac{(\Lambda_{1}/\hbar)^{2}}{(\omega-\omega_{0})^{2}+(\Lambda_{1}/\hbar)^{2}}$, which gives FWHM equal $2\Lambda_{1}$. This agrees with our calculations. E.g. for the green curve, i.e. $V_{ac}=100$ mV, transition time (period) $T=2\pi/\Omega$ is $680$ ps, which corresponds to $\Lambda_{1}=h/T\simeq 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.

V.2 Two-qubit gates

In Fig. 13(a) there are presented four basis states spanning the two-qubit subspace together with all six transitions—each with a different resonant frequency. The voltage modulations induce electron momentum oscillations and intervalley coupling amplitude modulations which enables us to obtain spin operations (blue and violet arrow pair), intervalley transitions (green and purple), or spin-valley swapping (red and yellow). This twofold control of the electron state allows to fully operate within the defined spin-valley two-qubit subspace. 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 $\omega/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 $\mathcal{K}=-1$, i.e. $K^{\prime}$ valley is occupied. In case of the valley index $\mathcal{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_{\mathcal{K}}NOT_{s}$. For the violet transition ($2.333$ meV) we get the opposite CNOT operation with the spin qubit flipped if $\mathcal{K}=1$, denoted by $\bar{C}_{\mathcal{K}}NOT_{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}NOT_{\mathcal{K}}$. The purple transition is performed for the opposite, spin-up, thus denoted by $\bar{C}_{s}NOT_{\mathcal{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 single-qubit operationShi (2003), e.g. the $R_{\pi/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_{\pi/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,\uparrow\rangle\leftrightarrow|1,\downarrow\rangle$. The complementary operation, induced by the yellow transition, interchanges the two remaining states: $|\!-\!1,\downarrow\rangle\leftrightarrow|1,\uparrow\rangle$. We denote it by $cSW\!AP$. All the mentioned two-qubit operations are represented in the spin-valley two-qubit subspace of $|\mathcal{K},s\rangle$ states by $4\times 4$ unitary matrices. Their explicit form can be found in the appendix.

V.3 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_{\mathcal{K}}NOT_{s}$ and $\bar{C}_{\mathcal{K}}NOT_{s}$ we arrive at single spin-NOT quantum gate, denoted as $NOT_{s}$, independent from the valley degree. Similarly, for simultaneous $C_{s}NOT_{\mathcal{K}}$ and $\bar{C}_{s}NOT_{\mathcal{K}}$ we get valley-NOT, $NOT_{\mathcal{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_{ac}^{x}\sin(\omega_{x}t)$ and $V_{2}(t)=V_{dc}+V_{ac}^{y}\sin(\omega_{y}t)$, by feeding both $G_{1}$ and $G_{2}$ gates. In Fig. 14 there are presented twofold transitions which are composed of two intervalley transitions for both spin orientations making up the valley-NOT operation (top: $\mathcal{K}$-qubit NOT gate), and two spin transitions for both valley occupations ($K$ and $K^{\prime}$) forming $s$-qubit NOT gate (bottom).

The both pumping frequencies $\omega_{x}$ and $\omega_{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 $\hbar(\omega_{x},\omega_{y})$ respectively. Whereas the voltage oscillation amplitudes pair $(V_{ac}^{x},V_{ac}^{y})$ 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_{ac}^{x}/V_{ac}^{y}=100\mathrm{mV}/62\mathrm{mV}$ for $\hbar(\omega_{x},\omega_{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 $C_{\mathcal{K}}NOT_{s}\cdot\bar{C}_{\mathcal{K}}NOT_{s}=\bar{C}_{\mathcal{K}}NOT_{s}\cdot C_{\mathcal{K}}NOT_{s}=\mathbf{1}_{2}\otimes NOT_{s}$ gives spin-NOT. Similarly, $C_{s}NOT_{\mathcal{K}}\cdot\bar{C}_{s}NOT_{\mathcal{K}}=\bar{C}_{s}NOT_{\mathcal{K}}\cdot C_{s}NOT_{\mathcal{K}}=NOT_{\mathcal{K}}\otimes\mathbf{1}_{2}$ results in valley-NOT. Both relations can be easily verified by multiplying matrices (included in the appendix) representing the particular operations.

V.4 Gate fidelities and qubits readout

In the course of transitions from Figs. 14 or 12 we can notice a minor oscillation structure of frequency $\hbar\omega$ related to a single cycle of pumping induced by the voltages oscillations. This can be viewed as a reference frame rotation with $\omega$ frequency in the standard RWA approximationZeuch et al. (2018). 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}NOT_{\mathcal{K}}$ and $NOT_{\mathcal{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. EstimatedWu et al. (2016) coherence times of electron valley and spin degrees of freedom in MoS2 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 overestimatedRivera et al. (2016). However, it scales with the system size as $\sqrt{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 nanotubesPei et al. (2012); Laird et al. (2013). 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. $|\mathcal{K}_{0},s_{0}\rangle=|\!-\!1,\uparrow\rangle$. 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 stateRohling and Burkard (2012); Osika and Szafran (2017). The electron is blocked in the same state as its neighbor: $|\!-\!1,\uparrow\rangle(1,1)|\!-\!1,\uparrow\rangle$ 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,\uparrow\rangle(1,1)|\!-\!1,\uparrow\rangle\rightarrow(0,2)|\mathcal{K}_{1},s_{1}\rangle|\!-\!1,\uparrow\rangle$ with $\mathcal{K}_{1}=1$ or $s_{1}=\,\downarrow$. 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 tight-binding 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.

Acknowledgements.

Author would like to thank Grzegorz Skowron and Paweł 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.

Appendix A Rashba coupling parameters

To calculate the Rashba coupling parameter matrix we utilize the 11-band model from [Ridolfi et al., 2015] with five $d$-orbitals in the Mo atom $d_{z^{2}}$, $d_{xy}$, $d_{x^{2}-y^{2}}$, $d_{xz}$, $d_{yz}$, and six $p$-orbitals for the S atoms, three for top ($t$) and three for bottom ($b$) layers: $p^{t,b}_{x}$, $p^{t,b}_{y}$, $p^{t,b}_{z}$. We add the atomic spin-orbit interaction of the form taken from [Roldán et al., 2014; Kośmider et al., 2013] with intrinsic parameters $\lambda_{\mathrm{Mo}}=0.086$ and $\lambda_{\mathrm{S}}=0.052$ eV.Kośmider et al. (2013) After calculating appropriate Slater-Koster elements we obtain the full Hamiltonian $22\times 22$ (together witch spin) [see Appendix B in Ridolfi et al., 2015], with additional onsite potentials $V_{t,b}$ for the top/bottom layers, serving as parameters. Now, it can be expressed in an infinite layer-form as a function of the momentum $H(k_{x},k_{y})$, simply by substituting each hopping $t_{ij}\rightarrow t_{ij}\exp\left(\imath\mathbf{k}\cdot\mathbf{R}_{ij}\right)$.

We calculate numeric value of $H(k_{x},k_{y})$ for different $V_{t,b}=\mp E_{z}d/2$ around the CB minimum, i.e. for $k=K=(\frac{4\pi}{3a},0)$ (or $k=K^{\prime}=(\frac{2\pi}{3a},\frac{2\pi}{\sqrt{3}a})$) and energy level $\mathcal{E}=2.22$ eV. By utilizing the Löwdin partitioning technique Löwdin (1962); Bochevarov and Sherrill (2006); Jin and Song (2011) we downfold it to a reduced Hamiltonian $6\times 6$ within our 3-band modelLiu et al. (2013); Pavlović and Peeters (2015), used in the simulations. The Schroödinger equation for the full $22\times 22$ block Hamiltonian is

[TABLE]

We perform downfolding by eliminating $\psi_{16}$ and arrive at representation $H^{\prime}_{6\times 6}\,\psi_{6}=\mathcal{E}\,\psi_{6}$ where

[TABLE]

which is equivalent to (14). Afterwards calculating (15), we finally obtain a numeric value of $H^{\prime}_{6\times 6}$ as a function of $E_{z}$. The vector $\psi_{6}$ is represented in the basis of orbitals $(d_{z^{2}},d_{xy},d_{x^{2}-y^{2}})\otimes(\uparrow,\downarrow)$, while $\psi_{16}=(d_{xz},d_{yz},p^{t}_{x},p^{t}_{y},p^{t}_{z},p^{b}_{x},p^{b}_{y},p^{b}_{z})\otimes(\uparrow,\downarrow)$.

After the procedure of downfolding to the 3-band model we obtain a $3\times 3$ matrix representing the Rashba Hamiltonian, indexed by the orbital numbers $\alpha$ and $\beta$, made up of $2\times 2$ spin blocks. Each of this block has the form $a_{0}\mathbf{1}_{2}+a\sigma_{x}+b\sigma_{y}$, which we write down as

[TABLE]

As expected, resulting specific numeric value $e|E_{z}\gamma_{\alpha\beta}$ of each block $(\alpha,\beta)$ is proportional to the external electric field $E_{z}$. Moreover, it is multiplied by the matrix with phase factors $e^{\pm\imath\eta}$ of the form equivalent to $e^{y}_{ij}\sigma_{x}-e^{x}_{ij}\sigma_{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 $\eta$ was disordered.

In this way we obtain the parameter matrix $\gamma^{\alpha\beta}_{R}=|e|E_{z}\gamma_{\alpha\beta}$, where

[TABLE]

We get the same results for the two remaining $K$ points. However, values obtained around $K^{\prime}$ point are slightly different

[TABLE]

Finally, we take to the model an approximate—average value of the $\gamma$ as written down in (2), remembering, however, of slightly different values between valleys. The obtained value corresponds to the Rashba coupling amplitude from [Kormányos et al., 2014a], where $\lambda_{R}=3.3\times 10^{-4}\,E_{z}$ (eV nm), for $E_{z}$ expressed in the $\mathrm{V}/\mathrm{nm}$ units. E.g. for $1$ V/nm and $k$ of order $\frac{4\pi}{3a}$ we get $\lambda_{R}k\sim 4.3$ meV ($a=0.319$ nm). While in our calculations, for a similar electric field, the $\gamma$ 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 [Pearce and Burkard, 2017] (or [Kormányos et al., 2014a]), calculated within the k.p model for $50$-nm-size (or $40$ nm) MoS2 QD. However here we have more than 10 times smaller dot, thus to obtain similar orbital effects relatively to the magnetic scale $\sqrt{\frac{\hbar}{eB}}$ (size of the Landau ground state), equal $\sim 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\times 10^{4}$ T, shows that Landau levels posses complicated self-similar structure, called the Hofstadter butterflyHofstadter (1976); Goldman (2009); Hunt et al. (2013); Wang et al. (2015). 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 $|\mathcal{K},s\rangle=|\mathcal{K}\rangle\otimes|s\rangle$. The first pair constitutes CNOT operations where the valley is the control qubit and the spin is the target one:

[TABLE]

While in the second pair, the spin controlled valley-NOT operations are:

[TABLE]

Finally, the SWAP and its complementary operations are given by:

[TABLE]

Note that performing jointly both SWAP operations is equivalent to the NOT operation on the both qubits $\sigma_{x}\otimes\sigma_{x}$, i.e. $SW\!AP\cdot cSW\!AP=cSW\!AP\cdot SW\!AP=NOT_{\mathcal{K}}\otimes NOT_{s}$.