In the following analysis we consider a quantum dot based on GaAs/Ga1−xAl_{x}As
heterostructure, but assume that the electron effective mass remains constant
across the heterogeneous boundary. This is a reasonable approximation for a
potential step of size 100 meV, as it amounts to only a 15% difference in the
effective mass for the Ga1−xAlxAs outer region in comparison to that of the inner
GaAs region [98]. As the the electron densities for the lower energy states are
negligible in the outer region, this difference does not contribute much to our
problem. Further, we assume a much stronger confinement in the axial direction
of the dot relative to that in its radial direction, thus justifying our two dimen-
sional model for the analysis [106]. Under these assumptions, the Hamiltonian
of a two electron cylindrical quantum dot system in a constant magnetic field B
can be written as

(4.1) H = X

i=1,2

1
2m^{∗}_{e}

~

p_{i}+q ~A_{i}2

+V (r_{i}) + gµ_{B}

¯

h S~_{i}·B~

+ q^{2}

κ|~r_{1} −r~_{2}|

where A~_{i} = ^{B}_{2} (0, r_{i},0) and

V (r_{i}) = 0 ; ri ≤r0

V_{0} ; r_{i} > r_{0} .

Here~r_{i} and ~p_{i} are the conjugate position and momentum of the i^{th} electron, m^{∗}_{e}
is the electron effective mass, q is the electron charge, g is the electron g-factor,
µ_{B} is the Bohr magneton, S~_{i} is the spin angular momentum of the i^{th} electron,
and κ= 4π, in SI units. Here, the value of quantum dot permittivity,, is given
an average value of 13.1_{0} throughout the material. In the operator form Eq.(4.1)
becomes

(4.2) Hˆ =

2

X

i=1

− ¯h^{2}
2m^{∗}_{e}

∂^{2}

∂r_{i}^{2} + 1
r_{i}

∂

∂r_{i} + 1
r_{i}^{2}

∂^{2}

∂φ^{2}_{i}

+ ω_{C}
2

Lˆ_{z}_{i} +1

2m^{∗}_{e}ω_{C}
2

2

r_{i}^{2}
+V (r_{i}) + ¯hω_{L}

2 σˆ_{z}^{(i)}

+ q^{2}

κ|r~_{1}−r~_{2}|
where ˆL_{z}_{i} = −i¯h_{∂φ}^{∂}

i, the z-component of the angular momentum operator of
the i^{th} electron, ω_{C} = _{m}^{qB}∗

e, the cyclotron frequency, and ω_{L} = _{m}^{qB}

0, the Larmor frequency.

By adding and subtracting a harmonic oscillator potential ^{1}_{2}m^{∗}_{e}ω_{0}^{2}r_{i}^{2} for each
electron, we can re-write the above Hamiltonian as

Hˆ = ˆH^{0}(1,2) + ˆH^{0}(1,2) (4.3)
The definitions of ˆH^{0}(1,2) and ˆH^{0}(1,2) are as given below.

Hˆ^{0}(1,2) = ˆH_{ho}(1) + ˆH_{ho}(2) (4.4)
where

(4.5)
Hˆ_{ho}(i) = − ¯h^{2}

2m^{∗}_{e}
∂^{2}

∂r^{2}_{i} + 1
r_{i}

∂

∂r_{i} + 1
r_{i}^{2}

∂^{2}

∂φ^{2}_{i}

+ω_{C}
2

Lˆ_{z}_{i}+ 1

2m^{∗}_{e}ω^{2}r_{i}^{2}
with ω =

q

ω_{0}^{2}+ ^{ω}_{2}^{C}2

. Similarly,

Hˆ^{0}(1,2) = ˆW(1,2) + ˆZ(1,2) + ˆC(1,2) (4.6)

where

Wˆ (1,2) =

2

X

i=1

V (r_{i})−1

2m^{∗}_{e}ω_{0}^{2}r^{2}_{i}

, the residue potential term

Zˆ(1,2) = hω¯ L

2

2

X

i=1

ˆ

σ^{(i)}_{z} , the Zeeman term and
Cˆ(1,2) = q^{2}

κ|~r_{1}−~r_{2}|, the Coulomb term.

The eigenvalues and eigenfunctions of ˆH_{ho}are well known [26] and are as given
below.

Hˆ_{ho} ϕ_{n,m}(r, φ) = En,m ϕ_{n,m}(r, φ) (4.7)
wheren = 0,1,2. . .and m= 0,±1,±2. . .are the radial and azimuthal quantum
numbers respectively, and

En,m = ¯hω(2n+|m|+1) + ¯hω_{C}m

2 (4.8)

ϕn,m(r, φ) = 1
a^{1+|m|}_{B} |m|!

r(|m|+n) ! πn! exp

−r^{2}

2a^{2}_{B} +imφ

r^{|m|}1F1[−n,|m|+1, r^{2}
a^{2}_{B}]

(4.9)
where a_{B} =q

¯ h

m^{∗}_{e}ω is the effective Bohr radius of the quantum dot.

Now, in a linear variational theory [107], we consider a trial wavefunction expressed as

|ψi=c_{1}|χ_{1}i+c_{2}|χ_{2}i+· · ·+c_{i}|χ_{i}i+· · ·+c_{d}|χ_{d}i (4.10)
Here d is the dimension of the basis set, coefficients c_{1}, c_{2}, . . . , c_{d} are variational
parameters and|χ_{i}iare the orthonormalized two-electron states. If s_{1} and s_{2} are
the spin quantum numbers of individual electrons (each having a value ^{1}_{2}), and if
s is the total spin quantum number of the two-electron system, then by angular
momentum addition rule s can take only two values viz. s = 0 corresponding to

|s_{1}−s_{2}|, ands= 1 corresponding to|s_{1}+s_{2}|. Fors= 0 state, thez-component of
the total spin operator, ˆS_{z}, can have only one value for its quantum number viz.

ms= 0 (singlet), and for the s= 1 state, it can take three values viz. ms = 0,±1 (triplet). The orbital part of the singlet state must be symmetric and that of the triplet states must be antisymmetric. These can be easily constructed out of one- electron wavefunctions by taking Slater permanent and determinant respectively.

For example, if|ϕ_{n}_{1}_{,m}_{1}iand|ϕ_{n}_{2}_{,m}_{2}iare any two distinct one-electron eigen states,
then

χ^{±}_{n}_{1}_{,m}_{1}_{,n}_{2}_{,m}_{2}(1,2)

= |ϕ_{n}_{1}_{,m}_{1}(1)i |ϕ_{n}_{2}_{,m}_{2}(2)i ± |ϕ_{n}_{2}_{,m}_{2}(1)i |ϕ_{n}_{1}_{,m}_{1}(2)i

√2 (4.11)

is a valid orthonormalized symmetric (antisymmetric) two-electron orbital wave-
function constructed out of them, except for the case (n_{1}, m_{1}) = (n_{2}, m_{2}) where,
an additional factor of √

2 must be taken care of. This, multiplied by their appro- priate spin-wavefunction counterpart (|si), forms the required basis wavefunctions in Eq. (4.10), given by

|χ_{i}i=
χ^{±}_{i}

|s = 0(1)i (4.12)

Since h

H,ˆ Lˆ_{z}i

= 0, the total z-component of the angular momentum, M =
m_{1}+m_{2}, must be the same for all terms in the trial wavefuntion. Similarly, since
hH,ˆ Sˆ^{2}

i

=

hH,ˆ Sˆz

i

= 0, every term in the trial wavefunction must have the same
value for the s and m_{s} quantum numbers. Thus, we can do variational analysis
for each combinations of M, sand m_{s} values separately. For each case, the linear
variational analysis reduces to solving an eigen value problem of the type

[H]{C}=E{C}. (4.13)

Here [H] is a square matrix of size d×d with elements H_{ij} =D
χ_{i}

Hˆ
χ_{j}E

, {C}

is a column matrix of size d×1 with unknown coefficients c_{i} of Eq. (4.10) as
elements, and E is their corresponding eigen energy. The matrix [H] can be
written as a sum of matrices [H^{0}] and [H^{0}] corresponding to the operators ˆH^{0}
and ˆH^{0} defined via Eq. (4.3). The elements of [H^{0}] for the case of singlet (+)

and triplet (−) states are given by Eq. (4.14). Due to the orthogonality of one-
electron wavefunctions|ϕ_{n,m}i, it may be noticed that [H^{0}] has non-zero elements
only along its diagonal.

H_{ij}^{0} =
D

χ^{±}_{n}_{1}_{,m}_{1}_{,n}_{2}_{,m}_{2}

Hˆho(1) + ˆHho(2)

χ^{±}_{n}_{3}_{,m}_{3}_{,n}_{4}_{,m}_{4}
E

=δ_{n}_{1}_{n}_{3}δ_{m}_{1}_{m}_{3}δ_{n}_{2}_{n}_{4}δ_{m}_{2}_{m}_{4}(En3,m3 +En4,m4)

(4.14)

Similarly, the elements of [H^{0}] are also evaluated. In the case of triplet states,
the result after substituting for |χ^{−}i from Eq. (4.11) becomes

(4.15)
H_{ij}^{0} =

D

ϕn1,m1(1)ϕn2,m2(2)

Hˆ^{0}(1,2)

ϕn3,m3(1)ϕn4,m4(2) E

−D

ϕ_{n}_{1}_{,m}_{1}(1)ϕ_{n}_{2}_{,m}_{2}(2)

Hˆ^{0}(1,2)

ϕ_{n}_{4}_{,m}_{4}(1)ϕ_{n}_{3}_{,m}_{3}(2)E

A little care must be given while evaluating the elements of [H^{0}] for singlet states.

This is because, we can build symmetric orbital states out one-electron eigenstates
even when n_{1} =n_{2} and m_{1} =m_{2}, and those states are different from the general
symmetric states given by Eq. (4.11). For the general states,

(4.16)
H_{ij}^{0} =D

ϕ_{n}_{1}_{,m}_{1}(1)ϕ_{n}_{2}_{,m}_{2}(2)

Hˆ^{0}(1,2)

ϕ_{n}_{3}_{,m}_{3}(1)ϕ_{n}_{4}_{,m}_{4}(2)E
+D

ϕ_{n}_{1}_{,m}_{1}(1)ϕ_{n}_{2}_{,m}_{2}(2)

Hˆ^{0}(1,2)

ϕ_{n}_{4}_{,m}_{4}(1)ϕ_{n}_{3}_{,m}_{3}(2)E

When n_{1} =n_{2} and m_{1} =m_{2} and n_{3} =n_{4} and m_{3} =m_{4}, the elements of [H^{0}] for
singlet states becomes

(4.17)
H_{ij}^{0} =

D

ϕn1,m1(1)ϕn1,m1(2)

Hˆ^{0}(1,2)

ϕn3,m3(1)ϕn3,m3(2) E

The final case is when quantum numbers are equal only on one side, say n_{1} =n_{2}
and m_{1} =m_{2}, the elements of [H^{0}] for singlet states are given by

(4.18)
H_{ij}^{0} =√

2D

ϕ_{n}_{1}_{,m}_{1}(1)ϕ_{n}_{1}_{,m}_{1}(2)

Hˆ^{0}(1,2)

ϕ_{n}_{3}_{,m}_{3}(1)ϕ_{n}_{4}_{,m}_{4}(2)E

All the elements of [H^{0}] were found out through numerical integration in the
Mathematica software. To simplify the case, we can split and write [H^{0}] as a
sum of three matrices corresponding to the terms defined in Eq. (4.6) as

[H^{0}] = [W] + [Z] + [C] (4.19)

Now, the matrix [Z] is diagonal due to the orthonormality of orbital part of the basis wavefunctions. Furthermore, these diagonal elements are also zero for ms = 0 states

|↑↓i±|↓↑i√ 2

. But for ms = ±1 states (|↑↑i,|↓↓i), it give rises to
a constant value, ±¯hωL for all the diagonal elements. Similarly, many of the
elements in [W] are zeros due to the othogonality condition when m_{1} 6= m_{3} or
m_{2} 6= m_{4}. Again, because of the azimuthal and exchange symmetries of the
Wˆ (1,2) term, all the non-zero [W] elements require only one–variable numerical
integrations. The evaluation of elements in the [C] matrix involves calculation of
terms like

C_{12,34}=D

ϕ_{n}_{1}_{,m}_{1}(1)ϕ_{n}_{2}_{,m}_{2}(2)

Cˆ(1,2)

ϕ_{n}_{3}_{,m}_{3}(1)ϕ_{n}_{4}_{,m}_{4}(2)E

. (4.20)

If we substitute ϕ_{n,m}(r, φ) = ρ_{n,m}(r)^{e}^{√}^{imφ}

2π in Eq. (4.20) we get,

(4.21)
C_{12,34}= q^{2}

2πκ Z ∞

r1=0

r_{1}dr_{1}ρ^{∗}_{n}_{1}_{,m}_{1}(r_{1})ρ_{n}_{3}_{,m}_{3}(r_{1})
Z ∞

r2=0

r_{2}dr_{2}ρ^{∗}_{n}

2,m2(r_{2})ρ_{n}_{4}_{,m}_{4}(r_{2})I(r_{1}, r_{2})
where

I(r_{1}, r_{2}) =
Z 2π

φ=0

dφ_{1}
Z 2π

φ=0

dφ_{2} exp[i(−m_{1}+m_{3})φ_{1}+i(−m_{2}+m_{4})φ_{2}]

|~r_{1}−~r_{2}| .
Using the multipole expansion of _{|~}_{r} ^{1}

1−~r2| and simplifying further,I(r_{1}, r_{2}) becomes
(4.22)
I(r1, r2) = 4π^{2} 1

r_{>}

∞

X

l=0

(r_{<}

r_{>})^{l}(l−m_{1}+m_{3})!

(l+m_{1}−m_{3})!(P_{l}^{m}^{1}^{−m}^{3}(0))^{2}

where,r_{>}=r_{1} andr_{<}=r_{2}whenr_{1} > r_{2}and vice versa, andP_{l}^{m} is the associated
Legendre function. Substituting this expression forI(r_{1}, r_{2}) in Eq. (4.21), we get

C_{12,34}= 2πq^{2}
κ

∞

X

l=|m1−m3|

(l−m_{1}+m_{3})!

(l+m1 −m3)!(P_{l}^{m}^{1}^{−m}^{3}(0))^{2}
Z ∞

r1=0

r_{1}dr_{1}ρ^{∗}_{n}_{1}_{,m}_{1}(r_{1})ρ_{n}_{3}_{,m}_{3}(r_{1})

"

Z r1

r2=0

dr_{2}ρ^{∗}_{n}_{2}_{,m}_{2}(r_{2})ρ_{n}_{4}_{,m}_{4}(r_{2})
r2

r_{1}
l+1

+ Z ∞

r2=r1

dr_{2}ρ^{∗}_{n}_{2}_{,m}_{2}(r_{2})ρ_{n}_{4}_{,m}_{4}(r_{2})
r_{1}

r_{2}
l#

(4.23) Thus we have replaced the original four-variable integration given in Eq. (4.20) with an infinite sum of two-variable integrations. But since the result is pretty accurate with only first few (say 30) terms of Eq. (4.23), this new formula helps us to achieve fast computation of the elements in [C].

Once all the elements of [H] are evaluated, we proceed to solve the eigen
value problem of Eq. (4.13) which will then give us the energy spectrum and
corresponding electron wavefunctions. We have done this by using the subrou-
tine Eigensystem available in the mathematica software. This process is done
separately for each value of quantum numbers (M, s) and it is then repeated for
values of B ranging from 0 to 30 Tesla. The electron density, η(r), and pair cor-
relation function, f_{pc}(~r), are then evaluated using these resultant wavefunctions,
following the definitions mentioned in reference [106] as given below

η(~r) =

2

X

i=1

hδ(~r−r~_{i})i
fpc(~r) =hδ(~r−r~1+r~2)i

(4.24)