Lagrange mesh and exact diagonalization for numerical study of semiconductor quantum dot systems with application in singlet-triplet qubits

3.7 MB·English

by  Tuukka Hiltunen

Lagrange mesh and exact diagonalization for numerical study of semiconductor
quantum dot systems with application in singlet-triplet qubits
Tuukka Hiltunen, Juha Ritala, Oona Kupiainen, Topi Siro, and Ari Harju
Department of Applied Physics and Helsinki Institute of Physics,
Aalto University School of Science, P.O. Box 14100, 00076 Aalto, Finland
(Dated: January 29, 2013)
Wepresentahighlyflexiblecomputationalschemeforstudyingcorrelatedelectronsconfinedbyan
3 arbitraryexternalpotentialintwo-dimensionalsemiconductorquantumdots. Themethodstartsby
1 aLagrangemeshcalculationforthesingle-particlestates,followedbythecalculationoftheCoulomb
0 interaction matrix elements between these, and combining both in theexact diagonalization of the
2 many-bodyHamiltonian. Weapplythemethodinsimulation ofdoublequantumdotsinglet-triplet
qubits. We simulate the full quantum control and dynamics of one singlet-triplet qubit. We also
n
use our method to provide an exact diagonalization based first-principles model for studying two
a
singlet-triplet qubitsand their capacitative coupling via thelong-distance Coulomb interaction.
J

8
PACSnumbers: 73.22.-f,81.07.Ta
2

]
l I. INTRODUCTION particle eigenstates are computed using the Lagrange
al meshmethodandarethenusedtocreatethemany-body
h Thedevelopmentofexperimentalmethodshasenabled basis in the calculations. The use of the Lagrange mesh
- method as a source of single-particle states allows us to
s the fabrication of “artificial atoms” with a controlled
e number of electrons, ranging from a few to a few hun- study very flexibly various forms of the confinement po-
m tential. For example, the relatively complex two-DQD
dred, confined in a tunable external potential inside a
. semiconductor.1–3Thesequantumdots(QD’s)havebeen case can be handled with ease using this method. Still,
t
a proposedasapossiblerealizationforthequbitofaquan- the one-body basis needed for good accuracy is much
m tum computer4,5. more compact in the Lagrange mesh than, e.g., in the
finite-difference formulation.
- A framework for using two-electron spin eigenstates
d as qubits was proposed by Levy in 20026. The two- This paper is organized as follows. In Section II, the
n
electrondoublequantumdot(DQD)spinstateshavenat- theoretical model used in our computations is briefly
o
c ural protectionagainstthe decoherence by the hyperfine discussed. In Section III, we introduce the Lagrange
[ interaction and allow a scalable architecture for quan- meshmethod for many-body problems in quantumdots.
tum computation7. The universal set of quantum gates The computational results are shown in Section IV. The
1
for two spinsinglet-tripletDQD qubits has been demon- single-particlestates andtheir convergencearediscussed
v
7 strated experimentally. These gates include one qubit inSectionIVA.InSectionIVB,wemodelthefullquan-
1 rotationsgeneratedbytheexchangeinteraction8andsta- tum control and dynamics of a singlet-triplet qubit. In
5 bilizedhyperfinemagneticfieldgradients9,andtwoqubit Section IV C, we use the Lagrange mesh method to cre-
6 operations using long distance capacitative coupling by ate a realistic first-principles ED model for studying the
1. the Coulomb10 interaction. interplay and entanglement of two singlet-triplet qubits.
0 Increatingtheinter-qubitgatesandoperations,quan-
3 tum entanglement is essential11. The aforementioned
1
capacitative coupling is one possible method to create II. MODEL
:
v entangled states and implement two-qubit operations
i in singlet-triplet qubits, the other possibility being ex-
X We model a lateral GaAs quantum dot system with
change based methods6,12. In the capacitative dipole-
the two-dimensional Hamiltonian
r
dipole coupling, the entanglement is achieved by differ-
a
ing charge densities in the singlet and triplet states that N (p +eA(r ))2 e2
resultindifferentCoulombrepulsionbetweenthequbits. Hˆ = j j +V(r )+V (r ) + ,
j Z j
2m 4πǫr
This conditioning can be used to create entangled states j=1(cid:20) ∗ (cid:21) j<k jk
X X
and to implement the two-qubit gates required for uni-
versal quantum computing7,10,13–15. where V (r ) = g µ B(r ) S is the Zeeman term
Z j ∗ B j j
·
Althoughother methods, like the variationalquantum with the effective GaAs g-factor g = 0.44. A is
∗
Monte Carlo16 andthe density functional theory17, have the magnetic vector potential, and m −0.067m and
∗ e
≈
shown to give reasonably accurate results, exact diago- ǫ 12.7ǫ are the effective electronmass and permittiv-
0
≈
nalization is still the most reliable technique for small ity in GaAs, respectively. In numerical work, it is con-
particle numbers. In this paper, we use the Lagrange venient to switch into effective atomic units by setting
mesh method18 and exact diagonalization to simulate m =e=~=1/4πǫ=1. In these units, energy is given
∗
one and two singlet-triplet qubit systems. The one- by Ha∗ ≈11.30 meV and length in a∗0 ≈10.03 nm.
2
In our computations, the external potential V(r) for Many different Lagrange meshes, mostly based on or-
quantum dot systems consists of several parabolic wells. thogonal polynomials or trigonometric functions, have
Aconfinementpotentialofnparabolicwellscanbewrit- been proposed19 both for finite intervals and the infinite
ten as intervals (0, ) and ( , ). The meshes can also be
∞ −∞ ∞
modified to distribute the mesh points optimally for a
1
V(r)= 2m∗ω021mjinn{|r−rj|2}, (1) particular system.20
≤ ≤ One of the most simple Lagrange meshes is the sinc
where r are the locations of the minima of the mesh.19 It is defined over the interval ( , ), but de-
j 1 j n
parabo{licw}e≤lls≤, andω is the confinement strength. The signed to treat fairly well localized wav−e∞fun∞ctions. The
0
kinks caused by the min-function are smoothed in the mesh points distributed uniformly around the origin are
potential.
x =a, a N 1, N 1 +1, ... , N 1 , (6)
a ∈ − 2− − 2− 2−
(cid:8) (cid:9)
III. METHOD and all the weights in the Gauss quadrature are λ = 1.
a
The Lagrange-sinc functions are
The Lagrange mesh method18 is a very efficient
method for solving the Schr¨odinger equation, and it has
sin[π(x a)]
thesimplicityofafinite-differencemeshcalculation,since L (x)=sinc(x a)= − .
a
− π(x a)
no integrations need to be performed. It also does not
−
suffer from the same limitations regarding to the con-
The matrix elements of the derivatives ∂ and ∂2 be-
finement potential as for example using the analytical x x
tween two sinc functions can be calculated analytically,
Fock-Darwin basis.
resulting in
We will use this technique to solve the eigenstates of
the one particle Hamiltonian,
∞
p2 (∂x)a′a = dxLa′(x)∂xLa(x)
Hˆ = +V(r), (2) Z−∞
2m 0 , a =a
∗ ′
omitting the magnetic vector potential here. The eigen- =
slattaitoens adroenethbeyntuhseedexaascatbdaiasgisonfoarlitzhaetiomnanteyc-hbnoidqyuec,aalcnud-  (−a1′)−a′a−a , a′ 6=a,
the Zeeman term and Coulomb interaction are included and 
in it.
π2 , a =a
− 3 ′
∂2 =
A. One-particle problem (cid:0) x(cid:1)a′a  −2((−a′1)aa′)−2a , a′ 6=a.
−
A set of N Lagrange functions Lk defined over an in- The potential energymatrix elements can be calculated
terval (a,b) is associated with N mesh points x (a,b)
k analyticallyforsomepotentials,but itturnsoutthatfor
∈
and a corresponding Gauss quadrature
the smooth potentials, these can be accurately approxi-
mated using the Gauss quadrature of Eq. (3) as
b N
dxf(x) λ f(x ). (3)
k k
Za ≈kX=1 Va′a ≈V(a)δa′a.
The Lagrange functions are infinitely differentiable real
Strictly speaking, this approximation breaks the varia-
functions, which are orthonormal,
tional principle. The validity of the Gauss quadrature
approximationis discussed in the Appendix.
b
dxL (x)L (x)=δ , (4) Generalized to two dimensions and an area
i j ij
Za ( L,L) ( L,L), the Lagrange-sinc functions are
−2 2 × −2 2
and satisfy the Lagrange conditions given by
Li(xj)=λi−1/2δij. (5) La(r)= NL sinc NL(x−xax) sinc NL(y−yay) , (7)
From the conditions of Eqs. (4) and (5), it follows that (cid:2) (cid:3) (cid:2) (cid:3)
where the N N mesh points are scaled to
the Gauss quadrature is exact for any product of two ×
Lagrange functions: ra = NLa, ax,ay ∈ −N2−1, −N2−1 +1, ... , N2−1 ,
b N (cid:8) (cid:9)
dxL (x)L (x)=δ = λ L (x )L (x ). having grid spacing h=L/N and weights λa =h2. The
i j ij k i k j k
Za k=1 matrixelementsoftheHamiltonian(2)betweenthebasis
X
3
functions in Eq. (7), in effective atomic units, are where the expansion coefficients α multiply the inter-
action matrix elements vabcd between the sinc basis
3πh22 +V(ra) , a′x=ax,a′y =ay functions. To calculate these, we start with the two-
dimensional Fourier transform of 1/r , namely
wheHrea′xVa′y(arx)ayis=theex0hht((e22−−((r11aan))′y′xaaa−−′y′xl−a−ayxpaa))yxo22tent,,,iaaaal′x′x′x. 6==6=Daaaixxxa,,,gaaao′y′y′yn=6=6=alaaaiyyyza,tion whereFθ(cid:20)isr11t2h(cid:21)e(akn)g==lebZZe0Rt∞2wdederrn112r2eZ−0ria2k1πn2·rd1d2kθ1e.2iUkrs1i2ncgost(hθ+eπJ)a,cobi-
12
of this Hamiltonian matrix gives the one-particle eigen-
Anger identity of Bessel functions,
statesand-energies. Theaccuracyoftheresultsobtained
can be tested by varying the number of mesh functions
N and the side length of the simulation square L. eizcosφ = ∞ inJ (z)einφ,
n
n=
X−∞
B. Many particles
leads to
After the one-particle eigenstates are obtained, these
can be used as the single-particle basis for solving the 1 (k)= ∞dr 2πdθ ∞ ( i)nJ (kr )einθ
12 n 12
F r −
eigenstates of the interacting many-body system by ex- (cid:20) 12(cid:21) Z0 Z0 n=
X−∞
act diagonalization. The N-particle Hamiltonian can be
∞ 2π
written in the second quantization formalism as =2π dr12J0(kr12)= .
k
Z0
1
Hˆ = εjaˆ†jaˆj + 2 Vijklaˆ†iaˆ†jaˆlaˆk The potential 1/r can now be written as the inverse
j i,j,k,l 12
X X
Fourier transform of 1 as:
+ Vi,jaˆ†iaˆj, (8) F r12
h i
i,j
X
1 1 1 2π
where εj are the energyeigenvaluesof the single-particle = −1 = d2k eik·r12
Hamiltonian, r12 F ◦F (cid:20)r12(cid:21) (2π)2 ZR2 k
(cid:0)1 (cid:1) 1
V = ψ ψ Vˆ(r,r ) ψ ψ (9) = dkxdky eikx(x2−x1)eiky(y2−y1). (10)
ijkl i j ′ k l 2π R2 k
Z
D (cid:12) (cid:12) E
are the matrix elements of t(cid:12)he Coulo(cid:12)mb two-body inter-
action Vˆ(r,r ) = 1/r r(cid:12) in the s(cid:12)ingle-particle basis, With the identity of Eq. (10), the integrations over dif-
′ ′
| − | ferent coordinates factorize in the interaction matrix el-
andV = ψ Vˆ(r) ψ ,whereVˆ(r)containstheZee-
i,j i j ement:
man interaDction(cid:12) and(cid:12)addEitional external potentials that
(cid:12) (cid:12)
are not included(cid:12) in E(cid:12)q. (2). N 1
The interaction matrix elements of Eq.(9) can be cal- vabcd= 2πL R2 dkk
culated as follows. Let ψ be the single-particle eigen- Z
i
functions expanded in the sinc basis of Eq. (7), ∞ dx sinc(x a )sinc(x c )eikxx1
1 1 x 1 x
× − −
Z−∞
ψi(r)= αiaLa(r) . ∞ dy sinc(y a )sinc(y c )eikyy1
a × 1 1− y 1− y
X Z−∞
The interaction matrix elements are then ∞ dx2sinc(x2 bx)sinc(x2 dx)e−ikxx2
1 × − −
V = dr dr Ψ (r )Ψ (r ) Ψ (r )Ψ (r ) Z−∞
ijkl ZR2 1ZR2 2 ∗i 1 ∗j 2 r12 k 1 l 2 ∞ dy2sinc(y2 by)sinc(y2 dy)e−ikyy2.(11)
× − −
= αia∗αjb∗αkcαld dr1 dr2 Z−∞
a,b,c,d ZR2 ZR2
X The sinc functions can be replaced by their integral rep-
1
La(r1)Lb(r2) Lc(r1)Ld(r2) resentation
× r
12
= αia∗αjb∗αkcαldvabcd , sinc(x)= 1 πdteixt,
a,Xb,c,d 2π Z−π
4
andtheintegralsoverxandycoordinatesareoftheform the threads with eachthreadcorrespondingto a value of
d. The threads then loop over the indices a,b and c,
I (k)= ∞ dxsinc(x a)sinc(x b)eikx and in the end the results of all threads in the block are
ab
− − summed with a parallel prefix sum algorithm to obtain
isigZn(−k∞)( 1)a b eika eikb , k 2π, a=b the final result.
2π(a−b) − − − | |≤ 6 InEq.(12),vabcd doesnotdependonthestateindices
 (cid:0) (cid:1) i,j,k,l,anditisbeneficialtocalculateitbeforehandand
= 21πeika(2π−|k|) , |k|≤2π, a=b sctisoiroenitfloinatiangtapboleinitnatrhitehmGePtUic,mtheme soirzye.oIfnthdeoutbablelepfroer-
Bysubs0titutingthisresultintoEq.(11,),|kt|h>eo2rπig.inalfour- saizNe i×s aNppmroexsihmiaste8lNy83.b3ygteigs.abFyotersa, w12hi×ch1fi2tsmienstho, tthhee
6 gigabyte global memory of the state of the art Tesla
dimensional integral over two planes reduces into a two-
cards,suchastheC2070. Wealsoutilizethefaston-chip
dimensional integral over a finite square in k-space,
shared memory by caching the expansion coefficients α
2π 2π 1 before the calculation. The GPU speeds up the matrix
vabcd =2 dkx dky elementcalculationbyafactorofaround13whendouble
k
Z−2π Z−2π precision arithmetic is used.
I (k )I (k )I (k )I (k )
× axcx x aycy y bxdx x bydy y Unfortunately, we have found the 12 12 mesh insuf-
×
2π K(θ) ficient for double quantum dot calculations if a realis-
= dθ dk tic distance between the minima is used. Therefore, the
Z0 Z0 calculation has to be divided so that the whole vabcd
I ( kcos(θ))I ( ksin(θ))
× axcx − aycy − matrix is not calculated at once. We lowered the mem-
×Ibxdx(kcos(θ))Ibydy(ksin(θ)) , (12) oryrequirementbycalculatingfirstViajkxl byfixingtheax
index in Eq. (12). As a consequence, it is sufficient to
where K(θ) = 2π/max( cos(θ), sin(θ)) is the radial
integration limit correspo|nding t|o|the sq|uare. The last calculate the vabcd also using a fixed ax index, and the
memory requirement is dropped to 8N7 bytes, which al-
form can be used in numerical calculations.
lowscalculationwitha 17 17mesh. The V elements
One canseethat inEq.(12), one obtains fivedifferent areobtainedbysummingt×heVax elementsoijvkelra . The
integrals depending on how many of the four functions ijkl x
sum is updated after the calculationofeachVax matrix
Iab have the same indices. In addition, the case with ijkl
tosavememory. Thismodificationofthealgorithmadds
two equal index pairs is naturally split into two cases,
someserialwork,whichslowsdownthecomputation,but
depending on whether the equal indices belong to the
a compromise between memory requirement and speed
same Cartesian component of k. In most cases, some
must be made.
further analytic work can be done to handle the angu-
The sum in Eq. (12) could be further divided to allow
lar integral. For instance, in the case when all the index
largermeshsizesbyfixingmoreindices,butthecomputa-
pairs differ, such that a = c , b = d , a = c and
x x x x y y
6 6 6 tiontime becomesfasta limiting factor. The calculation
b =d , the integrand can be written as a sum of terms
y y
6 of interaction matrix elements for the 24 lowest single-
of the form cos k[mcos(θ)+nsin(θ)] , and the angular
{ } particle states using a mesh size of 17 17 takes almost
part can be integrated analytically, and we are left with
×
two days. The exact diagonalization part is much faster
a one-dimensional numerical integral. In this way, we
than this.
are able to calculate the interaction matrix elements be-
tween the sinc basis functions, andthen for any external
confinement potential, Eq.(10)can be usedto construct
IV. RESULTS
V .
ijkl
It turns out that the calculationof V fromEq.(12)
ijkl
iscomputationallyverytime-consuming,becauseonehas A. Convergence of the single particle states
toloopoverfourindicesonboththeright-andleft-hand
sides of Eq. (12). Luckily, this basis change can be triv- In this section, we compute the single-particle eigen-
ially parallelized and a very efficient scheme can be ob- states of systems consisting of 1, 2, and 4 minima using
tained using graphics processing units (GPUs). theLagrangemeshmethod. Theaccuracyofthemethod
We performed the calculation of the interaction ele- and optimal parameters are discussed. In the following
ments in Eq. (12) with an Nvidia Tesla C2070 graphics sections, we then use the obtained single particle states
processing unit, which was programmed with CUDA21, in the actual many-body computations.
a parallel programmingmodel for Nvidia GPUs. On the First,westudiedtheconvergenceofthemethodinthe
GPU,thecomputationisparallelizedacrosstensofthou- analytically solvable case of just one parabolic well. The
sands of lightweight computational threads, which are confinement strength was ~ω = 4 meV. We computed
0
organized in independent blocks. In our parallelization the 24 first single particle energies with different mesh
scheme,eachblockcomputesoneelementofV . Inside parameters N and L and compared the results with the
ijkl
the block,the sumoverthe index d is parallelizedacross analytical Fock-Darwin eigenenergies. The relative dif-
5
est eigenvalue of R. The effect of the rounding can be
100 seen in Fig. 2.
e
c
en
nergydiffer10−5 eV)1105 δδ==02 meV
e m
ative10−10 V( 5
el
R
0
−50 0 50
10−15 x(nm)
5 10 15 20 25 30 35 40
N
FIG.2. (Coloronline)TheeffectoftheroundingontheDQD-
potential. The potential is shown in the x-axis. The minima
FIG. 1. (Color online) The convergence of the 24 first single
are located at x = ±40 nm. The confinement strength is
particleenergiesinthecaseofoneparabolicdotasafunction ~ω0 = 4 meV. Both non-rounded (δ = 0, red dashed curve)
of N (the grid being N ×N). The relative differences of the
and rounded (δ = 2 meV, black solid curve) potentials are
single-particle energies (computed with the Lagrange mesh
shown.
method)withtheanalyticalFock-Darwinenergiesareshown.
Thethinredcurvesshowthefirst24statescasewithL=280
The current maximum grid size in the computation of
nm(thethickredcurveshowstheaveragerelativedifference).
the V -elements is N = 17 due to the GPU memory
Thethickdashedcurveshowstheaveragerelativedifferences ijkl
in thesmaller grid case with L=200 nm. limitations (larger grids can in principle be computed,
butwiththeexpenseofconsiderablylongercomputations
times). As the accuracy of the method depends non-
ference of the energies can be seen as a function of the trivially on both the simulation area L and the grid size
grid size N in Fig. 1. N, the value of L was optimized.
Fig. 1 shows that given large enough N, the relative We comparedthe obtainedeigenenergieswith thoseof
difference of the energies converges to the order of the alargesystem(N =68andL=300nmorL=320nm)
numerical double precision accuracy in the L = 280 nm andchosethevalueforLthatgavethesmallesterrorwith
case. Theeffectofthesizeofthesimulationareacanalso respect to the more accurate large system. The relative
be seen in the figure. The smaller area case (the black difference of the energies as a function of L in the two
dashed curve, L = 200 nm) shows faster convergence dot case is shown in Fig. 3 and in Fig. 4 in the four-dot
with respect to N. However, the finite simulation area case. The potentials for the two- and four-dot systems
resultsinsomeerroraswell,andthustheconvergencein areillustratedintheinsetsofFigs. 3and4. Thetwo-dot
the L = 200 nm case stops before it reaches the double potential consists of parabolic dots with the distance 80
precision. nm between their minima. In the four-dot system, dots
The main topic of this paper is the simulation of 1 and 2 are 80 nm apart,dots 2 and3 are 120nm apart,
singlet-triplet qubits. We will first study one-qubit dy- and 3 and 4 are 80 nm apart.
namics and then use our model to simulate a system of Figs. 3 and 4 show that with N = 17 the optimal
two singlet-triplet qubits. Next we discuss the conver- value of L is between180 nm and200nm in the two-dot
gence of the method in these systems. caseandbetween260nmand290nminthefour-minima
In the potential in Eq. (1), the derivative of the po- case. Up to this point, the convergenceof the energiesis
tentialisnotcontinuous;themin-functioncausesanedge monotonous. WithtoosmallL,thewavefunction’leaks’
at the interface of two branches. This sharp edge can be out of the simulation area, and with too high L the grid
problematic in the Lagrangemesh method due to the fi- spacing becomes too large. The singularity like dips in
nite number of mesh points. To alleviate this, rounding the relative difference curves probably result from the
of the edges was used in the case of multiple dots. The fact that the errors due to finite L and N have different
rounding was found to speed up the convergence of the signs. At the dip, these errors nearly cancel each other
single particle states. out.
The rounding is achieved by defining a matrix R at
each grid point. R has the different dot potentials in
its diagonal. For example, in the case of four dots at
locations r ...r the diagonal entries are R = V =
1 4 11 1
1m ω2 r r 2 and so on. The non-diagonalentries are
2 ∗ 0| − 1|
constant δ and define the strength of the rounding. The
potentialattheparticulargridpointisgivenasthesmall-
6
to model a singlet-triplet qubit. The confinement
strength was ~ω = 4 meV and the distance between
0
the dots was a= r r =80 nm. Our DQD potential
1 2
| − |
is illustrated in the inset of Fig. 3.
e
nc10−2 The logical basis of a singlet-triplet qubit consists
e
er of the two lowest eigenstates, the singlet state, S =
diff 1 ( ), and the S = 0 triplet state, |T i =
gy √2 | ↑↓i−| ↓↑i z | 0i
er 1 ( + )(the arrowsdenotedirectionofthe elec-
en10−4 √2 |↑↓i |↓↑i
e tron spins). An arbitrary state ψ of the qubit can be
Relativ 80 nm written as | i
θ θ
10−6 |ψi=cos 2 |Si+sin 2 eiφ|T0i. (14)
(cid:18) (cid:19) (cid:18) (cid:19)
100 120 140 160 180 200 220
L(nm) Here, θ [0,π] and φ [0,2π). The state of the qubit
∈ ∈
canthusbevisualizedusingthesurfaceofaBlochsphere
FIG.3. Theconvergenceofthe24firstLagrangemeshsingle with S and T0 atthenorthandsouthpoles,andθand
| i | i
particle energies as a function of the simulation area length φ denoting the angles with respect to z and x-axes.
L in the two dot case. The two-dot potential is illustrated It should be noted that the DQD is not a true two-
in the inset. The rounding is set to δ = 2 meV. The grid level system. There are higher excited states as well,
size is 17 × 17 (N = 17). The relative differences of the and Eq. (14) is just an approximation. However, in our
Lagrange mesh energies with the energies of a larger system simulations, the weighs of the higher states were found
(L=300 nm and N =68, also computed with the Lagrange
to be negligible with practicalparametervalues,and the
meshmethod)areshown. Thethickcurveshowstheaverage
system can be considered as two level in this sense.
relative difference for the 24 states.
Universal quantum control of the qubit requires ro-
tations around at least two different axes in the afore-
100 mentioned Bloch sphere. In DQD singlet-triplet qubits,
rotations around the z-axis are controlled by the ex-
change interaction (the singlet-triplet energy difference)
difference10−2 gdJeont=se.rEaTtTeh0de−baxyEisaSomafnathgdenrerototitcaatfitiioeolndnsignarratodhuienenBdtlto∆hcehBzxsp-bahexetirwseeciesannthtbehnee
y
nerg 80 nm n=Jeˆz+µBg∆Bzeˆx, (15)
ativee10−4 120 nm and the frequency is9
el
R
f = J2+(µ g∆B )2/h. (16)
B z
10−6 80 nm p
Here,g = 0.44istheGaAsgyromagneticratio,µ the
B
−
150 200 250 300 Bohr magneton and h the Planck’s constant.
L(nm) The single particleeigenstatesare computedusing the
Lagrangemeshmethodandtheyarethenusedinthetwo
FIG.4. ThesameasinFig. 3butforthefourdotcase, with
particle ED-calculations. In our model, the z-rotations
thereference system being L=320 nm and N =68.
are created by detuning (a potential energy difference ǫ
between the minima of the dots) the two parabolic dots,
which lifts the degeneracy of the S and T states and
B. Singlet-triplet qubit | i | 0i
results in exchange interaction. The x-rotations are cre-
ated using a local magnetic field gradient that is taken
In this section, we use the Lagrange mesh and ED
into account by the Zeeman-term.
methods to simulate the time evolution of the state of
The detuning potential and the local magnetic field
a singlet-triplet DQD qubit. We demonstrate that, by
are modeled as step functions that are zero far away
applying local electric and magnetic fields in our model,
from the dot minima and have different signs in the two
wecanachievefullquantumcontroloverthe stateofthe
dots. We calculate the the matrix elements V =
qubit and reproduce realistic dynamics of the system in i,j,σ1,σ2
ψ V ψ , where V is either the detuning po-
our simulation. h i,σ1| | j,σ2i
tential or the Zeeman-interaction, in the eigenbasis
Weusedapotentialthatconsistsoftwoparabolicdots,
ψ obtained using the Lagrange mesh method (σ
i,σ
{ }
1 denotes the spin quantum number). The detuning
V(r)= m∗ω0min r r1 2, r r2 2 , (13) and the Zeeman-term are then taken into account in
2 {| − | | − | }
7
the two-body ED through the one-body operator Vˆ =
i,j,σ1,σ2Vi,j,σ1,σ2a†i,σ1aj,σ2. 1
The evolutionof the initial state ψ(0) of the qubit is
P | i
computed by propagation, using
0.8 J≈0
i J≠0
ψ(t+∆t) =exp Hˆ(t)∆t ψ(t) , (17)
| i (cid:18)−~ (cid:19)| i ) 0.6
〉
S
where Hˆ(t) is the (time-dependent) two-body Hamilto- (|
p
nian. The matrix exponent is computed using Lanczos 0.4
method. To study the evolution of the qubit’s state in
theBlochsphere,theanglesθandφin(14)areextracted
0.2
from ψ(t) by using the properties of the two-body spin
| i
operator Sˆ2, i.e. Sˆ2 S =0 and Sˆ2 T =2~2 T .
0 0
Thefirst24single|-piarticlestates|weirecomp|utiedusing 0
0 0.2 0.4 0.6 0.8 1
the Lagrange mesh method. The mesh parameters were
time (ns)
N =17 and L=210nm. The rounding was setto δ =2
meV.TheV -andV -elements(correspondingtoboth
ijkl ij FIG.5. (Coloronline)Thetimeevolutionofthesingletprob-
the detuning and the Zeeman term) were computed for
ability p(|Si) (red curve). The detuning is ǫ=4.3 meV, and
the 24 single-particle states. themagnetic fieldis∆Bz =0.4 T.Thispart isomitted from
We firstdemonstratethe controlofthe qubitina sim- the figure, as p(|Si) is constant during it. The simulation
ple case. In this simulation, the system is initially in the time is 1 nsand thetime step length ∆t=1 ps. Thedashed
singlet state. The dots are detuned so that the differ- cyan curveshows thenon-detunedcase with J ≈0.
ence between their energy minima is ǫ = V V = 4.3
2 1
−
meV. The detuning lifts the degeneracy of the singlet
and triplet states, resulting in an exchange energy of
J 3.748µeV.A magnetic fielddifference of∆B =0.4
z
≈
is then put between the dots and the system is let to
evolvefor1ns. Thesingletandtripletprobabilitieswere
computed by projecting the state of the qubit onto the
S2 operator.
The computed time evolution of the singlet probabil-
ities p(S ) can be seen in Fig. 5. Fig. 6 shows the
| i
evolution of the state of the qubit on the Bloch sphere.
In Fig 5, the detuned singlet probability oscillates be-
tween its maximum 1 and minimum 0.12. The singlet
probabilitynevergoestozeroduetothez-rotationdriven
bytheexchangeenergyJ 3.748µeV.Thefrequencyof
≈
the oscillation is f 2.618 GHz, which is very close to
≈
the value givenby Eq. (16), f 2.625GHz. In the non-
≈
detuned case, the probability oscillates between 1 and 0,
as expected. In this case too, the computed frequency
coincides very well with Eq. (16). FIG.6. (Coloronline)TheevolutionofthestateoftheDQD
Fig. 6showsthatinthedetunedcase,the planeofthe qubitontheBlochsphere. Thequbitisinitiallyinthesinglet
rotation is tilted from the S -T plane, as expected by state (at the north pole). A detuning of ǫ = 4.3 meV and
0
Eq. (15). The state never|reia|cheis T (the south pole) a magnetic field gradient of ∆Bz = 0.4 are turned on and
0
| i the state is let to evolve. The simulation time is 1 ns, and
during the simulation. The non-detuned case oscillates
∆t=1 ps. The thick red curve denotes the trajectory of the
between S and T , passing through the spin localized
| i | 0i qubit’s state on the sphere. The thin cyan curve shows the
states = 1/√2(S + T0 ) and = 1/√2(S non-detuned case, with J ≈ 0. The dashed black curve (the
| ↑↓i | i | i | ↓↑i | i−
T0 ). equator) shows thexy-plane,where p(|Si)=1/2.
| i
We also tried more complicated pulse sequences and
tracked the evolution of the state in the Bloch sphere.
One such is demonstrated in Fig. 7. Here, the detuning Inconclusion,ourmodel,basedonLagrangemeshED,
strength was oscillating, ǫ(t) = ǫ sin(2πft), where ǫ = canbeusedtosimulateGaAssinglet-tripletDQDqubits.
0 0
4.6meVandf =1GHz. Anon-triviallytimedependent We cansimulate the realisticfull quantum controlofthe
J causes the axis of the states’s rotation to change as a qubitstartingfromthefirstprinciples. Next,weproceed
function of time, which leads to quite complicated paths to use the model in studying two-qubit dynamics.
on the Bloch sphere.
8
150
200
100
r1 aA
r
50 A 2 150
m)
n 0 d
y( 100
B
−50 r
3
−100 aB r4 50
−150
−100 0 100
x(nm)
FIG.8. Contourplot of thetwo-DQD potential. Thesystem
consistsoffourQDs,dividedtotwoqubitsAandB.Thecon-
FIG. 7. (Color online) evolution of the state of the DQD finement strength is ~ω0 =4 meV. The qubit-qubit distance
qubit on the Bloch sphere. The qubit is initially in the is d=120 nm, and thedistance of the dots in thequbits are
singlet state (at the north pole). An oscillating detuning, aA=aB =80 nm. The contours are shown in meV.
ǫ(t)=ǫ0sin(2πft), where ǫ0 =4.6 meV and f =1 GHz, and
magneticfield∆Bz =0.4T,arethenapplied. Thesimulation
time is 1 ns, and ∆t=1 ps.
lowersthepotentialenergyinoneofthedotsofthequbit,
the singlet state chargedensity becomes more located in
this dot. However, if the detuning is not too high the
C. Two singlet-triplet qubits
tripletdensityisunaffectedduetotherepulsiveexchange
force in the spatially anti-symmetric triplet state.
The Lagrange mesh allows the study of interplay of
The singlet and triplet states have differing charge
singlet-triplet qubits. For example, the entanglement
densities, and hence the Coulomb repulsion between the
of two singlet-triplet qubits by the long distance dipole-
qubits depends on the states of the qubits. This condi-
dipoleinteractioncanbesimulatedusingthismethod. In
tioning creates an entangled state when the qubits are
thissection,wefirstcomputeandstudythelowesteigen-
evolved under exchange.
states of the two-DQD system, using different detunings
A bipartite state ψ ( and
AB A B A B
for the two qubits. The main topic of this section is the | i ∈ H ⊗ H H H
are the Hilbert spaces of the subsystems A and B) is
entanglementofsinglet-tripletqubits. We show thatour
an entangled state if it cannot be written as a tensor
model can be used to simulate the entangling procedure
product ψ = ψ ψ . In general, if the vec-
demonstrated recently by Shulman et al10. | iAB | iA ⊗ | iB
tor ψ is written in any orthonormal product basis
AB
We model the two-DQD system with an external con- | i
e e ,
i A j B ij
finementpotentialthatistheminimumoffourquadratic {| i ⊗| i }
wells,
ψ = M e e , (19)
AB ij i A j B
| i | i ⊗| i
1
V(r)= 2m∗ω02mj in4{|r−rj|2}. (18) Xi,j
≤ it is an entangled state if and only if the matrix of coef-
OursimulationsystemcanbedividedtoqubitsAandB. ficients, M= M , is not singular.
ij
Aconsistsofthewellsatr1 andr2,withthedotdistance Thedegreeo{fent}anglementcanbedeterminedbysome
aA = r1 r2 =80nm. Similarly,the inter dotdistance entanglement measure. One such measure is the con-
| − |
of the qubit B is aB = r3 r4 = 80 nm. The inter currence. In case of pure states, and two-level systems
| − |
qcounbfiitnedmisteanntcsetriesnggitvhenisb~yωd0==|r42m−eVr3.|T=he12p0ontemn.tiaTlhies (wqhuebrietsλ)1AaannddλB2,acroentchuerreeignecnevCaluisegsiovfenmaastrCix=M2†√Mλ1.λI2t,
illustrated in Fig. 8. The inter-qubit distance and the is easy to see that this simplifies to the formula
confinement strangth are large enough that there is no
tunneling between A and B, so the qubits interact only C =2 det(M). (20)
| |
through the Coulomb repulsion of their electrons. Also,
the qubits interactmainly via the electronsin the dots 2 Concurrence can also be generalized to mixed states11.
and 3, as the inter-dot distances are quite large. Concurrence assumes values between 0 and 1. A non-
The qubits A and B can become entangled due to the zeroC isapropertyofanentangledstate,andthehigher
factthatundertheexchangeinteraction,thechargeden- the value of C, the higher the degree of entanglement.
sitiesofthe S and T statesdiffer. Whenthedetuning The maximally entangled Bell states have C =1.
0
| i | i
9
In our two singlet-triplet qubit system, the Hilbert thequbitsarelettoevolveinthexy-plane,theweightsin
spaces are given as the two lowest eigenstates of the M matrix obtain phase factors proportionalto these
a DQD-system, = = S , T . The energies.
A B 0
M matrix is thuHs obtaineHd by pro{j|eciti|ngi}the four- Let the system be initially in the state with M = 1
ij 2
electron wave function onto the computational basis i,j 1,2 , i.e. ψ = . The system
A B
∀ ∈ { } | i | ↑↓i ⊗| ↑↓i
SS , ST , T S , T T ; M = SS ψ , M = is then let to evolve. If we approximate the projected
0 0 0 0 11 12
{| i | i | i | i} h | i
ST ψ , M = T S ψ and M = T T ψ . Hamiltonian to be diagonal, the time dependence of the
0 21 0 22 0 0
ahnIdt |TshioSuldarebeeignheontsetda|ttehisaotfwthheilfeouthr-ephastratitcels|e|SiS2Soip,e|rSaTto0ir,, Ecoeffi=ciEents aMndijsiosogniv.eInnsaesrtiMngijt(ht)es=e in21eEiqE.ij(t/2~0,)wyiheeldres
0 11 SS
| i
T T is not. Indeed, it is not given as an eigenstate by the formula for the concurrence,
0 0
| i
theLanczositeration. Inordertodotheprojectionsonto
1
the computational basis, the state T0T0 was generated C(t)= 2 2cos(∆ t/~), (22)
| i E
using localized magnetic fields. 2 −
As T = 1 ( ), T T = T T can p
be wr|itt0ein as√2 |↑↓i−|↓↑i | 0 0i | 0iA⊗| 0iB where, ∆E =ESS +ET0T0 −EST0 −ET0S.
In Eq. (22), the parameter ∆ represents the cou-
E
1 1 pling between the qubits. Eq. (22) shows that the
T T = +
0 0 A B A B entanglement indeed arises from the differences in the
| i 2|↑↓i ⊗|↑↓i 2|↑↓i ⊗|↓↑i
1 1 charge densities. If all the computational basis states
+ A B+ A B. (21) have identical charge densities, the Coulomb repulsion
2|↓↑i ⊗|↑↓i 2|↓↑i ⊗|↓↑i
between the two qubits, γ, is the same for all these
Inthis decomposition, T0T0 iswrittenusingthe Sz =0 states. In this case, the energies of the computational
| i
eigenstates. These Sz = 0 eigenstates can be generated basis states are ESS = 2ES +γ, ET0T0 = 2ET0 +γ and
using strong localized magnetic fields in the four dots E =E =E +E +γ. and∆ =0. Theconcur-
ST0 T0S S T0 E
of the two-qubit system. For example, A B rence is thus zero, i.e. there is no entanglement between
| ↑↓i ⊗| ↑↓i
is obtained as the ground state of a system where the the qubits. Inthe detuned case,the states havedifferent
magnetic field is up in the first dot, down in the second, densities and different values of the Coulomb repulsion.
up in the third and down in the fourth (the Zeeman- Hence,∆ =0, andthe concurrenceoscillatesaccording
E
6
termalignesthespinsoftheelectronsalongthemagnetic to Eq. (22).
fields). The single-particle states and the V -elements were
ijkl
Decompositions similar to Eq. (21) can be written for again computed using the Lagrange mesh method. The
theotherthreestatesaswell. Indeed,inthenon-detuned simulation cell area was L = 280 nm, and the mesh size
case,thesingletstatesgivenbysuchdecompositionswere was N =17. The rounding of the edges was set to δ =2
found to be the same eigenstates that Lanczos iteration meV. The 24 first one-particle states were used in the
would find. However, the aforementioned magnetic field four-particleEDcomputations. Asintheonequbitcase,
scheme for creating the S = 0 eigenstates can only be the detuning and local magnetic field matrix elements
z
usedtocreatestatesthathaveidenticaldensityinthetwo V = ψ V ψ were computed in order to
i,j,σ1,σ2 h i,σ1| | j,σ2i
dotsofthequbits. Hence,itisnotwellsuitedforcreating achieve full quantum control over the qubits.
the detuned singlet states. Fortunately if the detuning First, we study the lowest eigenstates of the four-
is in the practical operation regime of DQD-qubits, the electron system. Without detuning (ǫ = V V = 0
A 1 2
−
|T0T0idensityremainssymmetricwithrespecttothetwo and ǫB = V4−V3 = 0) the six first states (with Sz =0)
dots of the qubits. are close to each other in energy.
Thus, the computational basis can be created as fol- The ground state is SS (s = 0), and the next two
| i
lows. Thestates SS , ST and T S aregivenaseigen- statesare ST and T S (s=1). Thenextthreestates
| i | 0i | 0 i | 0i | 0 i
states by Lanczos and they can be identified by their given by Lanczos are what we call the triplet states, su-
spin. The state T T is generatedusingthe decomposi- perpositions of s = 0, s = 1 and s = 2 eigenstates (in
0 0
| i
tion Eq. (21). The wave function can then be projected thefour-particlecase,therecanbeseveralS2 eigenstates
onto this basis, and the concurrence can be computed withgivenquantumnumberssands )Thesethreeeigen-
z
according to Eq. (20). states of S2 are so degenerate that Lanczos mixes them.
In the scheme where the qubits A and B are first The three triplet states share the same energy as T T
0 0
| i
brought to the xy-plane and then let to evolve under (which also is not an S2 eigenstate, but another linear
exchange(usedforexamplebyShulmanetal.10),wecan combination of the triplet states). The electrons of the
deriveasimpleanalyticformulaforthetime dependence firstsixstatesaresymmetricallylocatedinthefourdots,
of the concurrence. one electron in each.
In the absence of magnetic fields, the Hamiltonian of With non-zero detuning, one begins to see differences
the two-qubit system is close to a diagonal one in the in the charge densities of the lowest eigenstates. Fig.
basis SS , ST , T S , T T (this was verified nu- 9 shows the effect of the sign of the detunings on the
0 0 0 0
{| i | i | i | i}
merically). The diagonal entries of the projected Hamil- ground state SS . The charge densities of the lowest
| i
tonian are the energies E , E , E and E . As states(givenbyLanczos)inthedetunedcase,ǫ =ǫ =
SS ST0 T0S T0T0 A B
10
FIG.9. Theeffectofthesignofdetuningonthedensityofthegroundstate|SSi. Inthemiddle: thenon-detuned|SSidensity.
Thenumbersinthemiddleplot refertothedots1-4andAandBdenotethequbits. Thesignsofthedetunings,ǫA =V2−V1
and ǫB = V3−V4, are shown by the axes ǫA and ǫB. In the upper left corner: ǫA = 4.35 meV and ǫB = −4.35 meV (dots 1
and 3 havelow potential). In theupperright corner: ǫA=ǫB =4.35 meV (dots1 and 4 havelow potential). In thelower left
corner: ǫA = ǫB = −4.35 meV (dots 2 and 3 have low potential). In the lower right corner: ǫA = −4.35 meV and ǫB = 4.35
meV (dots 2 and 4 havelow potential).
4.35 meV are shown in Fig. 10. Fig. 10) the dots that are the furthest apart from each
other have the lowest potential. This facilitates the lo-
calization of the singlet state in these dots, as it reduces
the Coulomb repulsion. In the cases when dots 2 and 4
or 1 and 3 are in the low potential, the singlet localizes
onlyinthe furtherawaydots,1and4. Inthecasewhere
dots2and3havethelowdetuning,therearefouridenti-
cal peaks, the singlets cannotlocalize in the neighboring
dots due to the Coulomb repulsion.
Thedipole-dipoleentanglementeffectreliesonthedif-
ferences of the singlet and triplet densities. Thus, to
this end, the optimal detuning configuration should be
as in Fig. 10, the dots furthest away are detuned to low
potential energy. The densities in Fig. 10 show that
the singlets localize to the dots that have lower poten-
tial (dots 1 and 4). The fourth plot represents all the
triplet states and is identical also to the T T density.
0 0
| i
It shows four identical peaks, with exactly one electron
in eachdot. When the detuning is further increased,the
singles localize fully to the low lying dots, and at very
high detunings, the triples start to localize as well.
FIG. 10. The charge density of the lowest eigenstates of the
DQD system. The detunings are ǫA =ǫB =4.3 meV (dots 1 With high detuning (ǫA and ǫB above 5 meV), |SSi,
and4havelowpotential). Theupperleftplotshowsthe|SSi ST , and T S were still lowest in energy. However,
0 0
| i | i
state, upper right the |ST0i state, lower left the |T0Si state, the triplet states were not the next three in this case.
and lower right the density of the triplet states (identical to There are states lower in energy than the triplets, in-
thedensity of |T0T0i). cluding other instances of the states SS , ST and
0
| i | i
T S (i.e. states that are of the form X Y ,
0 A B
| i | i ⊗ | i
There is a difference in the densities depending on where X and Y are s = 0 or s = 1 eigenstates of
which of the dots have the low potential, as can be seen the two|-eliectron|Si2-operator). Next we consider states
inFig. 9. IntheupperrightcornerofFig. 9(andalsoin SS , ST , T S , and T T that have the lowest en-
0 0 0 0
| i | i | i | i