Double quantum dots defined in bilayer graphene

TL;DR

This paper investigates double quantum dots in bilayer graphene using atomistic and continuum models, revealing how electron wave functions and interdot coupling influence the quantum states and energy spectrum.

Contribution

It introduces a detailed analysis of electron wave functions and interdot coupling effects in bilayer graphene quantum dots, combining atomistic and continuum approaches.

Findings

Wave functions have opposite parities on sublattices in both layers.

Ground state wave functions switch from bonding to antibonding with interdot distance.

Two-electron states show symmetric or antisymmetric spatial configurations.

Abstract

Artificial molecular states of double quantum dots defined in bilayer graphene are studied with the atomistic tight-binding and its low-energy continuum approximation. We indicate that the extended electron wave functions have opposite parities on each of the sublattices at both graphene layers and that the ground-state wave function components change from bonding to antibonding with the interdot distance. In the weak coupling limit -- the most relevant for the quantum dots defined electrostatically -- the signatures of the interdot coupling include -- for the two-electron ground state -- formation of states with symmetric or antisymmetric spatial wave functions split by the exchange energy. In the high energy part of the spectrum the states with both electrons in the same dot are found with the splitting of energy levels corresponding to simultaneous tunneling of the electron pair from one dot to the other.

Figures

Tables

Table 1. Table 1: The symmetries of the two-electron states, including the spatial Ψ s subscript Ψ 𝑠 \Psi_{s} , the valley Ψ v subscript Ψ 𝑣 \Psi_{v} and the spin Ψ σ subscript Ψ 𝜎 \Psi_{\sigma} factors where the total wave function Ψ ​ ( 1 , 2 ) = Ψ s ​ Ψ v ​ Ψ σ Ψ 1 2 subscript Ψ 𝑠 subscript Ψ 𝑣 subscript Ψ 𝜎 \Psi(1,2)=\Psi_{s}\Psi_{v}\Psi_{\sigma} . The table corresponds to Fig. 7 (d), where the first column orders the state in a growing energy order. Levels 7-8 and 9-10 are degenerate. The color of the lines corresponds to the z 𝑧 z component of the spin. The list corresponds to the dominant contribution to the CI wave functions.

no.	${\Psi }_s$	${\Psi }_v$	${\Psi }_{\sigma }$
16	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K\left(2\right)$	$\uparrow \left(1\right)\uparrow \left(2\right)$
15	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K\left(2\right)$	$\uparrow \left(1\right)\downarrow \left(2\right)+\downarrow \left(1\right)\uparrow \left(2\right)$
14	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)-\uparrow \left(1\right)\downarrow \left(2\right)$
13	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\downarrow \left(2\right)$
12	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)+K^{\prime }\left(1\right)K\left(2\right)$	$\uparrow \left(1\right)\uparrow \left(2\right)$
11	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)-K^{\prime }\left(1\right)K\left(2\right)$	$\uparrow \left(1\right)\uparrow \left(2\right)$
9-10	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)+K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)+\uparrow \left(1\right)\downarrow \left(2\right)$
9-10	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)-K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)-\uparrow \left(1\right)\downarrow \left(2\right)$
7-8	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)-K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)+\uparrow \left(1\right)\downarrow \left(2\right)$
7-8	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)+K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)-\uparrow \left(1\right)\downarrow \left(2\right)$
6	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)+K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\downarrow \left(2\right)$
5	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K\left(1\right)K^{\prime }\left(2\right)-K^{\prime }\left(1\right)K\left(2\right)$	$\downarrow \left(1\right)\downarrow \left(2\right)$
4	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K^{\prime }\left(1\right)K^{\prime }\left(2\right)$	$\uparrow \left(1\right)\uparrow \left(2\right)$
3	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K^{\prime }\left(1\right)K^{\prime }\left(2\right)$	$\uparrow \left(1\right)\downarrow \left(2\right)+\downarrow \left(1\right)\uparrow \left(2\right)$
2	$l\left(1\right)r\left(2\right)+r\left(1\right)l\left(2\right)$	$K^{\prime }\left(1\right)K^{\prime }\left(2\right)$	$\downarrow \left(1\right)\uparrow \left(2\right)-\uparrow \left(1\right)\downarrow \left(2\right)$
1	$l\left(1\right)r\left(2\right)-r\left(1\right)l\left(2\right)$	$K^{\prime }\left(1\right)K^{\prime }\left(2\right)$	$\downarrow \left(1\right)\downarrow \left(2\right)$

Equations

$$V_{QD}(x,y)=-V\exp (-((|x|-d{)}^2+y^2)/R^2),$$

$${\hat{H}}_{1e}=\sum _{i,\sigma }{W}_i^{\sigma }{\hat{c}}_{i,\sigma }^{†}{\hat{c}}_{i,\sigma }+\sum _{ij,\sigma {\sigma }^{\prime }}\left({\hat{c}}_{i,\sigma }^{†}{\hat{c}}_{j,{\sigma }^{\prime }}{t}_{ij}+h.c\right),$$

$$t_{ij}\to t_{ij}{e}^{i\frac{e}{\hslash }{\int }_{{\mathbf{r}}_{\mathbf{i}}}^{{\mathbf{r}}_{\mathbf{j}}}\mathbf{A}\cdot \mathbf{dl}},$$

$$W^{\sigma }(x,y)=V_{QD}(x,y)+\frac{W_B}{2}{\tau }_z+\frac{1}{2}{\mu }_{B}g{\sigma }_z{B}_z,$$

$$\Psi (\mathbf{r})=({\varphi }_A(\mathbf{r}),{\varphi }_B(\mathbf{r}),{\varphi }_{B^{\prime }}(\mathbf{r}),{\varphi }_{A^{\prime }}(\mathbf{r}){)}^T,$$

$${\hat{H}}_D=\left(\begin{array}{cccc}{V}_{QD}& \pi & t_{\perp }& 0\\ {\pi }^{†}& V_{QD}& 0& 0\\ t_{\perp }& 0& V_{QD}& {\pi }^{†}\\ 0& 0& \pi & V_{QD}\end{array}\right)+\frac{W_B}{2}{\tau }_z,$$

$${\tau }_z=\left(\begin{array}{cc}{\mathbb{I}}_2& 0\\ 0& -{\mathbb{I}}_2\end{array}\right),$$

$${\hat{U}}_P=\left(\begin{array}{cc}-{\sigma }_z& 0\\ 0& -{\sigma }_z\end{array}\right)\hat{P},$$

$${\hat{J}}_z={\hat{L}}_z{\mathbb{I}}_4+\frac{\hslash }{2}{\tau }_z-\frac{\hslash }{2}\left(\begin{array}{cc}{\sigma }_z& 0\\ 0& -{\sigma }_z\end{array}\right),$$

$${\hat{H}}_W=W_P\hslash v_F{a}_D\left(\begin{array}{cc}-{\sigma }_z& 0\\ 0& -{\sigma }_z\end{array}\right){\nabla }^2,$$

$$\Psi =\sum _{i=1}^M{d}_i\left({\psi }_{i1}({\mathbf{x}}_1)\otimes {\psi }_{i2}({\mathbf{x}}_2)-{\psi }_{i2}({\mathbf{x}}_1)\otimes {\psi }_{i1}({\mathbf{x}}_2)\right)$$

$${\hat{H}}_{2e}({\mathbf{x}}_1,{\mathbf{x}}_2)=\sum _{i=1}^2{\hat{H}}_{1e}({\mathbf{x}}_{\mathbf{i}})+\frac{\kappa }{|{\mathbf{r}}_{12}|},$$

$${\hat{H}}_{2e}=\sum _{ij}{\hat{d}}_i^{†}⟨{\psi }_i|\hat{H}|{\psi }_j⟩{\hat{d}}_j+\frac{1}{2}\sum _{ijkl}{\hat{d}}_i^{†}{\hat{d}}_j^{†}{\hat{d}}_k{\hat{d}}_l{V}_{ijkl},$$

$$V_{ijkl}=\kappa ⟨{\psi }_i({\mathbf{x}}_1){\psi }_j({\mathbf{x}}_2)\frac{1}{|{\mathbf{r}}_{12}|}|{\psi }_k({\mathbf{x}}_1){\psi }_l({\mathbf{x}}_2)⟩,$$

$$V_{ijkl}$$

$$V_{QD}(x,y)=$$

$$-V\exp (-((|x|-d{)}^2+y^2)/R^2)\cdot \left(1-\varphi \frac{x}{d}\right),$$

Full text

I Introduction

Studies of artificial molecules formed by carrier orbitals extended over double quantum dots (DQDs) focus on carrier tunneling and interactions vqds ; vqds1 ; ab1 ; ab2 ; ab3 ; paulisb . A particular attention is payed to the spin-related phenomena in the context of the quantum information processing ldp using exchange interaction due to the interdot coupling burkard ; paulisb . The coupling with nuclear spins decoh1 ; decoh2 ; decoh3 ; decoh4 which limits the electron spin coherence times motivated studies of DQDs with holes as spin carriers douho ; bohu as well as on systems based on silicon sili and carbon, including nanotubes laird and graphene tr . In graphene the carrier storage is hampered by the Klein tunneling klein but for nanoribbons the lateral confinement opens the transport gap ribbons-gap that makes the carrier storage possible dqdsingraphene . Quantum dots in ribbons are influenced by edge effects and disorder nrb5 . An alternative medium is the bilayer graphene review for which the perpendicular electric field opens the band gap gap in the energy spectrum and allows for carrier confinement by lateral fields biqd0 ; biqd1 ; biqd2 ; biqd3 ; biqd4 ; bizebr .

In this work we consider the bilayer graphene review and formation of extended orbitals within the DQDs. In the low-energy continuum approach Hamiltonian eigenstates possess a definite total angular momentum biqd0 with a different orbital angular momentum for each of the four sublattices described by the wave function components. As we discuss below for a lower – point symmetry of the DQD system – the Hamiltonian eigenstates components correspond to opposite spatial parities for each of the sublattices. Thus wave function can be bonding on one sublattice and antibonding on the other. The effect induces a complex dependence of the spectrum on the DQD distance that is similar to the antibonding heavy hole ground state found for vertical self-assembled quantum dots ab1 ; ab2 ; ab3 ; rev1 .

In the low energy part of the two-electron spectrum for the weak coupling case - the typical one for the electrostatic quantum dots – the electrons are localized in separate QDs and the ground state at zero magnetic field $B$ is nearly 16-times degenerate with respect to the valley and the spin. For nonzero magnetic field the energy levels correspond to wave function that can be approximately described by a product of the spatial, valley and spin components. Only the symmetry of the spatial part against the electron interchange influences the energy of the states, and we find that the energy levels shift in pairs with $B$ that correspond to the opposite symmetry split by the exchange integral. In the high energy part of the spectrum the states corresponding to both electrons in the same dot are found. These energy levels also shift in $B$ in pairs which is a manifestation of a collective two-electron intedot tunneling that forms bonding and antibonding two-electron orbitals. For stronger interdot coupling the DQD spectrum resembles the one of a single quantum dot with a ground state triplet zarenia2e .

This paper is organized as follows. In the Theory Section (II) we describe the model system, the atomistic and continuum approaches, as well as the configuration interaction method for the electron pair. In Section III we discuss formation of the extended orbitals and their nature as artificial molecular states and next we discuss the spatial, spin and orbital symmetries of the two-electron states. Summary and conclusions are given in Section IV.

II Theory

II.1 Model structure

We consider a flake in form of a stretched hexagon (see Fig. 1) with the bilayers in Bernal stacking review and two quantum dots defined by external potential. We consider two sizes of the flake. The smaller flake has a length of $L=64.6$ nm and height $h=27.6$ nm, with the horizontal edge length of $L^{\prime}=48.8$ nm, and the total number of carbon atoms $119928$. The larger flake is characterized by $L=89.46$ nm, $L^{\prime}=63.9$ nm, $h=44.51$ nm, and contains $261000$ atoms. Two different sizes of the flake for fixed external potential defining the quantum dots are introduced for the discussion of the effects of the coupling between the dot-confined states and the edge of the flake.

We consider a flake with an armchair boundary for which no edge-localized energy levels appear in the spectrum at the neutrality point HexTriang ; bizebr ; zeb2 . The quantum dots for electrons of the conduction band are defined by external potential which is taken in form

[TABLE]

where the origin ($x=y=0$) is placed in the center of the flake, $V$ is the quantum dots depth, $R$ is the effective radius and $2d$ is distance between the centers of the QDs. In this work, we fix $R$ at $4$ nm, for this value the wave functions localized in the dots vanish off the edges of the larger flake - for which the valley mixing effect of the armchair edge is negligible. For the smaller flake the effect of the edge for the confined states is still present and we use this fact below for description of the valley mixing effects on the single- and the two-electron spectra.

II.2 Tight-binding model

For most of the calculations we use the atomistic tight-binding Hamiltonian review given by

[TABLE]

where $\hat{c}_{i,\sigma}^{\dagger}$ ( $\hat{c}_{i,\sigma}$) creates (annihilates) an electron with spin $\sigma$ at ion $i$ and $t_{ij}$ is the hopping parameter. At zero magnetic field the hopping parameter takes the value of $t_{ij}=-2.6$ eV for the in-plane nearest neighbor atoms and $t_{ij}=0.3$ eV for the vertical dimers between the layers review . In order to account for the magnetic field perpendicular to the layers we include the Peierls phase into the hopping parameters,

[TABLE]

where $\mathbf{B}=\mathbf{\nabla}\times\mathbf{A}$ with $\mathbf{A}=(0,B_{z}x,0)$. The potential $W^{\sigma}_{i}$ in Eq. (2) has the form

[TABLE]

where the first term is given by Eq. (1) and the second one is the potential difference between the layers with $\tau_{z}=\pm 1$, with $+$ ($-$) for the upper (lower) layer. This difference opens the energy gap in the bilayer flake by the asymmetry introduced between the layers review . In our calculations we set $W_{B}=300$ meV which is equivalent to applying the vertical electric field $F\approx 0.9$ V/nm. For this value we obtain the gap of about $E_{g}\approx 141$ meV. The last term in Eq. (4) introduces the Zeeman interaction, where $g=2$ is the Landé factor and $\sigma_{z}$ stands for the Pauli matrix.

II.3 Continuum approximation

For analysis of the single-electron envelope wave functions, we also consider the low-energy approximation in the continuum model. Near a single valley the electron wave functions are written in form of the four component spinor review

[TABLE]

where $A,B$ ($B^{\prime},A^{\prime}$) represents the $A$ and $B$ sublattices of the upper (lower) layer respectively. The Dirac Hamiltonian review around the $\mathbf{K}$ valley takes the form

[TABLE]

where $\pi=v_{F}\left(p_{x}+ip_{y}\right)$, $(p_{x},p_{y})$ is the momentum operator, $v_{F}=\frac{3ta}{2}=0.84\cdot 10^{6}\frac{m}{s}$ – the Fermi velocity, and $t_{\bot}=0.3$ eV is the interlayer coupling term. The $\tau_{z}$ operator defined as

[TABLE]

that assigns +1 (-1) for upper (lower) layer respectively opens the gap by applying the voltage difference ($W_{B}$) between the layers.

The quantum dots potential (1) is symmetric with respect to the point inversion at the origin. For this potential the Hamiltonian commutes with a generalized parity operator

[TABLE]

where $\hat{P}$ changes the sign of the spatial coordinates of scalar functions $\hat{P}\phi(\mathbf{r})=\pm\phi(\mathbf{-r})$, and $\mathbb{I}_{n}$ is the $n\times n$ identity matrix.

For a single circular quantum dot the total angular momentum operator

[TABLE]

commutes with the Hamiltonian biqd0 , where $\hat{L}_{z}$ is the angular momentum operator. The second term is the layer index operator with $\tau_{z}$ defined by (7). The last term corresponds to the pseudospin review .

In order to find the eigenstates of the Hamiltonian (6) we use the finite element method with triangular elements and cubic shape functions. The fermion doubling fast-varying spurious solutions Wilson are eliminated with the additional Wilson Wilson term introduced to the Hamiltonian (6)

[TABLE]

where $W_{P}=0.015$ is the dimensionless Wilson parameter Wilson and $a_{D}=1.32$ nm is the discretization constant. $W_{P}$ parameter was fine-tuned to remove the spurious solutions from the discussed energy range leaving the actual states nearly unaffected.

II.4 Configuration-interaction method

The wave functions obtained with the tight-binding approach are used in the configuration interaction method tbci ; tbci2 . The atomistic approach when applied to the exact diagonalization method naturally accounts for the intervalley scattering induced by the short range component of the Coulomb interaction vm1 ; vm2 ; vm3 .

For the description of the two-electron states we use the solution of the one-electron eigenproblem (2). Then, we expand the two-electron wave function in the basis of $M$ Slater determinants

[TABLE]

with spin-orbitals $\psi_{i}(\mathbf{x}),\ i=1,2,\dots M$ where $\mathbf{x}=(\mathbf{r},\sigma)$ represents the orbital and spin coordinates. The number of Slater determinants $M$ is ${K}\choose{2}$, where $K$ is the number of dot-localized single-electron spin-orbitals. The two-electron Hamiltonian is

[TABLE]

where $\kappa=e^{2}/(4\pi\epsilon\epsilon_{0})$. The dielectric constant $\epsilon=6$ is taken from Ref. sicc for graphene grown on SiC. Using, the form of wave function given by (11), we arrive at the Hamiltonian of the form

[TABLE]

where $\hat{d}^{\dagger}_{i}$ creates an electron in the $i$-th spin-orbital. The one-electron energy is taken into account with the matrix elements $\langle\psi_{i}|\hat{H}|\psi_{j}\rangle$. The Coulomb matrix element

[TABLE]

for the single-electron wave functions given by linear combinations of atomic orbitals $p_{z}$ give

[TABLE]

For the Coulomb integral we apply the two-center approximation 2c $\langle p_{z}^{a}({\bf r}_{1})p_{z}^{b}({\bf r_{2}})|\frac{1}{|{\mathbf{r_{12}}}|}|p_{z}^{c}({\bf r_{1}})p_{z}^{d}({\bf r_{2}})\rangle=\frac{1}{r_{ab}}\delta_{ac}\delta_{bd}$ for $a\neq b$ and for the single center integral $a=b$ we take 16.522 eV tbci2 .

III Results

III.1 Single-electron energy levels

The single-electron spectrum obtained with the atomistic tight-binding is displayed in Fig. 2 as a function of the depth of the dots for the interdot distance $2d=20$ nm for zero magnetic field. The color of the lines indicates the extent of the electron localization within QDs – the red ones are localized in the QD potentials, and the grey energy levels correspond to the states which are localized outside the dots. Fig. 2 shows that the latter ignore the QD potential. The vertical line in Fig. 2 indicates the potential value $V=250$ meV that is taken for further calculations. Each of the energy levels plotted in Fig. 2 is nearly fourfold degenerate with respect to the spin and valley.

For identical quantum dots the splitting between the two lowest energy levels is a result of the interdot coupling. These energy levels tend to a further degeneracy at large $d$ – see Fig. 3. Typically, in the coupled quantum dots for small $d$ the bonding and antibonding electron orbitals are formed and the splitting energy in III-V systems is a monotonic function of $d$ burkard . This is not the case that we observe in Fig. 3(a) – where crossings of the two lowest energy levels are seen.

The energy spectrum obtained with the continuum approximation of Fig. 3(c,d) is similar to the exact results of Fig. 3(a,b). In Fig. 3(c) we plotted by the color of the line the average value of the total angular momentum (9) calculated for the envelope function of the continuum Hamiltonian. At large $d$, when the interdot tunneling is negligible and for $d=0$, where the potential has a rotational symmetry, $\hat{J}_{z}$ eigenvalue is a good quantum number. The separate components of the envelope wave function correspond then to the angular momentum quantum number biqd0 of $(m,m+1,m,m-1)$ for sublattices A,B,B’ and A’ respectively. The ground state for the circular QD corresponds to $m=1$ and the lowest excited state with $m=0$ in accordance with Ref. biqd0 .

For a general interdot distance $d$ the total angular momentum is no longer quantized. However, each of the components has a strict symmetry with respect to the point inversion through the center between the dots. The total wave function is an eigenstate of $\hat{U}_{P}$ and the eigenvalues $\pm 1$ are marked in Fig. 3(d) with the color of the line.

In order to analyze in a more detail formation of the extended orbitals by the single-dot wave functions we constructed the localized (ionic) orbitals of the left $l$ and right $r$ quantum dots using the tight-binding approach. The wave functions $l$ and $r$ were constructed separately from 1) the two lowest energy states $\phi_{1}$, $\phi_{2}$ and 2) the third and fourth states $\phi_{3}$, $\phi_{4}$. The ionic functions were taken as a superposition of the Hamiltonian eigenstates, $l_{1,2}=(\phi_{1}+\exp(i\alpha)\phi_{2})$, $r_{1,2}=(\phi_{1}+\exp(i(\alpha+\pi))\phi_{2})$, where the phase $\alpha$ was taken to maximize the electron localization at the left or right side of the origin. Then for Fig. 3(b) the extended states were produced as a constructive and destructive interference of the $l$ and $r$ functions $\phi^{1,2}_{c}=\frac{1}{\sqrt{2}}\left(l_{1,2}+r_{1,2}\right)$ and $\phi_{d}^{1,2}=\frac{1}{\sqrt{2}}\left(l_{1,2}-r_{1,2}\right)$, respectively. A similar operation was performed starting from the third and fourth energy level wave functions $\phi_{3}$ and $\phi_{4}$. The mean values of the energy calculated for $\phi_{c}$ and $\phi_{d}$ wave functions were plotted by blue and orange lines in Fig. 3(b), respectively. The red line shows the exact result that are the same as in Fig. 3(a). Note that, there is a one to one correspondence between the $\hat{U}_{P}$ eigenvalue and the superposition sign taken in Fig. 3(b).

For the scalar electron envelope function in a III-V material $\phi_{c}$ and $\phi_{d}$ would correspond simply to the bonding and antibonding orbitals. For multicomponent wave functions of opposite symmetry the situation is more complex. In particular, for holes in the artificial molecules formed by vertical quantum dots ab1 ; ab2 ; ab3 ; rev1 the wave functions need to be described by multicomponent wave function for each of the valence bands that become degenerate at the $\Gamma$ point. In these systems the parity of the heavy hole component is opposite to the light hole component. For large interdot distances the interdot tunneling is carried by the light hole component a bonding state at the light hole component is formed which triggers an antibonding heavy hole ground state ab1 ; ab2 ; ab3 ; rev1 . In III-Vs’ the crossing of the states of opposite parities is observed only for the vertical ab1 ; ab2 ; ab3 and not the lateral couping of the quantum dots rev1 .

For the lateral DQD studied here the origin of the non monotonic behaviour of the spacing between the lowest energy levels of Figs. 3(c,d) can be understood by inspection of the wave function. In Fig. 4 we plotted the real part of the wave functions for energy levels and $d$ values marked by the corresponding letters in Fig. 3(d). For the plots at large interdot distances [cf. Fig. 4(c,e,f)] – each of the wave function components depends on the azimuthal angle as $\Re\left[\exp(iL\Phi)\right]$ as for the angular momentum eigenfunctions with the quantum number $L$. In particular, in the lowest-energy state that is twofold degenerate at large $d$ [cf. Fig. 4(e,f)] the total angular momentum quantum number is $m=1$, and the $\hat{U}_{P}$ eigenvalue is +1 in Fig. 4(e) and -1 in Fig. 4(f). As the interdot distance is varied [cf. Fig. 3(d)] the energies of the states change but in a complex manner since two components of the wave functions are even in terms of $\hat{U}_{P}$ operator while the other two are odd, and the two parities have opposite consequences for energies as $d$ is modified. Moreover, the contributions of the components change. For instance: the state with the wave function given in Fig. 4(f) has a dominant A’ component which is antibonding. For lower $d$ the A and B components [Fig. 4(d)] are increased at the expense of A’ one. The black energy level in Fig. 3(d) becomes the ground state level between the interdot distance 8 nm and 16 nm. The varied contributions of the components involve avoided crossings between the energy levels of the same $\hat{U}_{P}$ parity – and the identity of the lines that avoided cross can be traced by the average angular momentum [Fig. 3(c)]. At $d=0$ the sequence of $J_{z}$ eigenvalues in the order of the energy is the same as in the large $d$ limit.

In order to analyze the effect of the molecular coupling on the energy levels and the crossing observed in Fig. 3(a-d) we plotted in Fig. 3(e) the contribution of the in-plane kinetic energy – for the hopping between A-B and A’-B’ sublattices within each layer that are due to the $\pi$ operator in the Hamiltonian given by Eq. (6). Figure 4(d,e,f) indicates that for the ground state the coupling of the single-dot energy levels – corresponding to $m=1$ each, appears for the largest interdot distances at the $B$ component which corresponds to orbital angular momentum quantum number $L=m+1=2$. The state that forms a bonding component there [Fig. 4(e)] with a positive eigenvalue of $\hat{U}_{P}$ operator decreases in the energy [Fig. 3(e)] due to the interdot coupling. When the dots get closer the components with orbital angular momentum $L=m=1$ at A and B’ sublattices form a molecular orbital which is antibonding (bonding) for the $\hat{U}_{P}$ eigenstate with eigenvalue +1 (-1). Due to the inverse contribution to the energy the first crossing at the ground-state energy level is observed near $2d=18$ nm [Fig. 3(a,d)].

The in-plane kinetic energies for the two states cross for even larger distance [Fig. 3(e)] of $2d=21$ nm. For the interlayer hopping between A and B’ sublattices included the crossing of kinetic energies occurs closer [Fig. 3(f)] to the intersection of the energy levels [Fig. 3(d)].

The second crossing of the ground-state energy levels near $2d=8$ nm results from the activation of the interdot tunnel coupling for the A’ component of the angular momentum $L=m-1=0$ [Fig. 4(d-f)]. Concluding, the ground-state energy level crossings for the states of opposite $\hat{U}_{P}$ parity results from the subsequent switching of the interdot tunnel coupling for separate components of the wave function at the sublattices.

The magnetic field dependence of the energy spectrum is given in Fig. 5. At this energy scale the results for the smaller [Fig. 5(a)] and the larger flake [Fig. 5(b)] become distinguishable. For the smaller flake we notice avoided crossings between the energy levels of opposite valleys which result from the valley mixing effect of the armchair edge of the flake [at $B\simeq 0$ and $B\simeq 0.4$ T] The valley mixing opens avoided crossing between the states of the same spin and $\hat{U}^{P}$ eigenvalue, and outside of the avoided crossings the energy levels in Fig. 5(a) can be attributed with the valley and spin quantum numbers. For the larger flake the valley mixing effects are not resolved on this scale.

III.2 Two-electron spectra

III.2.1 Low $d$ limit

The energy spectrum taken for a small interdot distance of $2d=6$ nm is displayed in Fig. 6(a). At $B=0$ the ground-state is threefold degenerate. One of the electrons occupies the $K$ and the other the $K^{\prime}$ valley. The valley degree of freedom allows the electrons to acquire the same spatial distribution with the three possible components of the total spin. From the dominant contributions to the wave function one concludes that the approximate form of the wave function for the three-fold degenerate ground state can be put in a separable form (normalization skipped) $\Psi_{gs}=\psi(1)\psi(2)\left(K(1)K^{\prime}(2)-K^{\prime}(1)K(2)\right)\Psi_{T}$, where $\psi$ is the spatial orbital for the positive $\hat{U}_{P}$ eigenvalue, and $\Psi_{T}$ is one of the spin triplet wave functions for $S_{z}=-\hbar$, 0 or $\hbar$. The triplet ground state at $B\neq 0$ is split only by the Zeeman effect. The first excited state is non degenerate spin singlet that corresponds to electrons in opposite valleys $\Psi_{4}=\psi(1)\psi(2)\left(K(1)K^{\prime}(2)+K^{\prime}(1)K(2)\right)\Psi_{S}$, where $\Psi_{S}=\left(\uparrow_{1}\downarrow_{2}-\uparrow_{2}\downarrow_{1}\right)$. The second excited state at $B=0$ is two-fold degenerate. Both the electrons occupy either the $K^{\prime}$ valley (the energy level that goes down in $B$), $\Psi_{5}=\psi(1)\psi(2)K^{\prime}(1)K^{\prime}(2)\Psi_{S}$, or the $K$ valley (the energy that goes up). Still the spatial orbital is the same for both electrons which is only possible due to the opposite spins of the states – they are both spin singlets.

For the strong coupling case of Fig. 6(a) the two dots effectively form a single system and the spectrum resembles the one calculated for a single dot in Fig. 6(b). The difference for the circular dot is that the valley unpolarized singlet goes higher in the energy at $B=0$. The avoided crossing with the $KK$ energy level as well as the avoided crossing between the two lower singlets at $B=0$ is due to the valley mixing effects of the armchair boundary condition. The ground state is still the triplet – which was previously found for a single-dot study of Ref. zarenia2e .

III.2.2 Weak coupling

For electrostatic quantum dots the coupling between the dots is usually weak, so this case deserves a closer inspection as the one which is the most likely to be encountered in an experimental situation. For a weak coupling of $2d=20$ nm the states can be separated into subgroups with both electrons in separate dots or with both electrons in the same dot. The latter states appear higher in the energy spectrum and they are characterized by a stronger electron-electron interaction. The right panel of Fig. 7 displays the energy spectrum in a wider energy range, as a function of the external field, and the color of the lines display the average electron-electron interaction energy. In Fig. 7(a,b) we display a zoom of the spectrum at the energy levels with the electrons in the same dots, for the smaller flake Fig. 7(a) and for the larger one Fig. 7(b). The spectrum of Fig. 7(a) is very close to the one displayed in Fig. 6(b) – for a single dot in a smaller flake. In Fig. 7(b) the spectrum is similar, only the avoided crossing for opposite valleys are closed. The closing of the avoided crossing in the two-electron spectra Fig. 7(b,d) (for the larger flake) with respect to Fig. 7(a,c) (for the small flake) is of the same origin as closing the avoided crossing of the single-electron spectra of Fig. 5(b) with respect to Fig. 5(a), i.e. it results from removed coupling of the quantum-dot confined energy levels with the edge. The considered armchair edge of the flake is not equivalent with respect to the A and B, and respectively A’ and B’ sublattices. The result of this non-equivalence is intervalley mixing of the states coupled to the edge. When the coupling is removed for a larger distance between the dots and the edge, the valley mixing disappears and the state of opposite valleys change their order in a crosssing instead of an avoided crossing.

The difference between the single QD [Fig. 6(b)] result and the DQD [Fig. 7(a)] is that for the latter, the energy levels shift in pairs. The electron couple can be stored by the left $l$ or the right $r$ dot, with the factor of the spatial wave function that can be put in one of the forms $\left\{l(1)l(2)\pm r(1)r(2)\right\}/\sqrt{2}$. The tunnel coupling between the dots leads to splitting of the energy levels into bonding and antibonding pairs that are observed in Fig. 7(a).

The ground-state of the two-electron spectrum of Fig. 7 corresponds to separated electrons and is displayed in Fig. 7(c) and (d) for the smaller and for the larger flake respectively. Here, the two electrons occupy different spatial orbitals $l$ or $r$. For that reason the Pauli exclusion does not forbid them to occupy any of the spin-valley combinations, hence the ground-state at $B=0$ is nearly 16 times degenerate – for the larger flake [Fig. 7(d)]. For the smaller flake [Fig. 7(c)] at $B=0$ the energy levels split into four quadruples. In the lowest (highest) quadruple the both electrons levels occupy the lower (higher) energy level: they have the four spin states to occupy hence the number of the nearly degenerate states. In both the lowest and highest quadruple a spin-singlet has a slightly lower energy than the spin triplet.

For the larger flake as well as for the smaller one at higher $B$ the single-electron energy levels shift away from the valley mixing avoided crossing and the two-electron spectrum forms groups depending on the valley configurations. By analysis of the CI components of the wave function we found that the wave functions at high $B$ are approximately separable into product of spatial, valley and spin factors $\Psi(1,2)=\Psi_{s}\Psi_{v}\Psi_{\sigma}$ which are listed in Table 1 for the energy order that is found in Fig. 7(b) for $B=0.3$T. In the group marked by 1 (3) both electrons occupy the $K^{\prime}$ ($K$) valley. Each group is formed by 4 energy levels with the splitting that is due to the Zeeman interaction. Since the valley degree of freedom is frozen, the structure of each of the group is identical with the one found for two-electron quantum dots in III-V materials, with the spin-singlet $S$ and spin-triplet energy level of zero $z$ component of the spin $T_{0}$ split by the exchange energy. Here the splitting due to exchange energy is of the order of 9 $\mu$eV.

In the central group of energy levels – marked by $2$ in Fig. 7(b) the electrons occupy opposite valleys - hence a nearly constant dependence on $B$, which here is only due to the Zeeman interaction. Each of the $S_{z}=0$ states is two-fold degenerate (energy levels 7-10 in Table 1) – only the symmetry of the spatial wave function with respect to the interchange of electrons influences the energy of the state via interdot tunneling effect, and there are two spin and valley factors for both symmetric and antisymmetric states. For each pair of energy levels of group (2) at high $B$ the one that is lower in the energy corresponds to a symmetric spatial wave function. The splitting energy is nearly the same as the one between the red energy levels of group 1 and 3. For all the energy levels that shift in pairs in Fig. 7(d) the splitting is due to a difference of the expectation value of the electron-electron interaction energy calculated for a symmetric and antisymmetric spatial wave function, i.e. to the exchange integral burkard .

IV Asymmetric DQD

For the symmetric system of quantum dots we described above, the splitting of the energy levels is uniquely due to the interdot tunnel coupling. In an experimental situation a deviation from the ideal symmetry is inevitable. For a generalization we consider the case of an asymmetric quantum dots. The asymmetry is introduced with confinement potential that includes an in-plane electric field

[TABLE]

with $V$ and $R$ kept unchanged. The value of the electric field is controlled by a dimensionless $\phi$ parameter, and $d$ in the denominator of the expression in brackets keeps the potential of the left and right dot in a fixed offset as the interdot distance $2d$ is varied. The in-plane field in an experimental situation corresponds to a bias between the dots that is introduced to induce the flow of the current.

The single-electron spectrum presented in Fig. 8(a) indicates that the in-plane electric field splits the degeneracy of the energy levels at $B=0$. For large interdot distance [Fig. 8(b)] the single-electron states are localized either in the left or the right quantum dots. In presence of the in-plane field the angular momentum is no longer quantized, but the quantum mechanical expectation values are preserved in the weak coupling limit with respect to the symmetric case [cf. Fig. 8(b) and Fig. 3(c)]. Moreover, the parity symmetry of the external potential is broken by the in-plane electric field. The Hamiltonian eigenstates of the $\hat{U}_{P}$ operator, which opens avoided crossing between the corresponding energy levels as the interdot distance is varied [Fig. 8(b)]. Nevertheless, for a small interdot distance the tunnel coupling prevails over the interdot asymmetry and the results for both the energy levels and the angular momenta are similar to the ideal case [cf. Fig. 8(b) and Fig. 3(c)].

Figure 9 shows the spectrum in the external magnetic field. The asymmetry of the confinement potential increases the energy spacing between the degenerate doublets at $B=0$ displacing the crossing of the $K$ and $K^{\prime}$ valley energy levels to higher values of the magnetic field (compare with Fig. 5).

For the low-energy states of the two-electron system at $d\simeq 20$ nm the electrons occupy separate dots. The energy level in one of the dots is shifted up and the other down, so no change to the ground-state energy level structure is observed unless the offset of the quantum dot potentials exceeds 60 meV, where the states with the electrons in the same dot occupy the same dot with the energy spectrum of the form presented in Fig. 6(b).

V Summary and conclusions

We have considered formation of extended orbitals in bilayer graphene quantum dots using an atomistic tight-binding and the continuum approach. The various angular momentum in each of the four sublattices for the total angular momentum eigenstates of a circular quantum dot evolve into various spatial parities for each of the wave function components for the double quantum dot system. The symmetry amounts in mixed bonding and antibonding character of wave function on separate sublattices and a complex dependence of the energy spectrum on the interdot distance. Both the flake large enough to be considered infinite for the dot-localized states and a smaller flake for which the valley mixing effects of the boundary could be described.

For the electron pair we used the configuration interaction method based on the atomistic single-electron wave functions. We have reproduced the limit result of the spectrum for a single quantum dot when the DQDs are close to one another and identified the effects of the interdot tunnel coupling between the dots for the more realistic weaker interdot tunnel coupling. The effects include – the splitting of energy levels in the high energy part of the spectrum where the electron pair forms bonding and antibonding two-electron orbitals, and near the ground-state – the exchange energy that splits the symmetric and antisymmetric pairs of energy levels that shift parallel in the external magnetic field.

Acknowledgments

This work was supported by the National Science Centre according to decision DEC-2013/11/B/ST3/03837 and by the Flemish Science Foundation (FWO-VL). Calculations were performed in the PL-Grid Infrastructure.