Since the three QDs and the magnetic field are cylindrically symmetric, the wave function, the spin density, and the paramagnetic current density can be repre-sented as
Ψq(r) = e−ilθφq(r, z) , q ≡ {nlσ} (3.1) nσ(r) = nσ(r, z) =
Nσ
q |φq(r, z)|2, (3.2) jp(r) = − h
m(r, z, εq)
N
q
l|φq(r, z)|2eθ, (3.3)
where n is the principal quantum number, l = 0,±1, ±2, · · ·, is the quantum number of the projection of angular momentum onto the magnetic field axis, i.e., the z-axis and eθ is the azimuthal unit vector. The KS equations are then reduced to a 2D problem in the (r, z) coordinates as
HlσKSφq(r, z) = εqφq(r, z) , (3.4)
where the KS Hamiltonian is now defined by
In 2D setting, the solution domain for (3.4) is again expressed by the same notation as that of 3D, that is, Ω = ΩInAs ∩ ΩGaAs ⊂ R2. We choose the domain ΩGaAs sufficiently large so that the wave function is negligibly small at the boundary of ΩGaAs. By symmetry, the domain Ω can be reduced to Ω = {(r, z) : 0 ≤ r ≤ rmax,−zmax ≤ z ≤ zmax} for sufficiently large rmax > 0 and zmax > 0 as shown in Fig. 1.1.
The explicit formula for the potential Vxc(r) in (2.13) is extremely complex in 3D coordinates. Transforming it to the (r, z) space is prohibitively lengthy and impractical. We use all the original formulas (2.13)-(2.19) for calculating Vxc(r) in the 3D space and then obtain the potential in the (r, z) coordinates, i.e., Vxc(r, z) for (3.4).
Since we are dealing with the hard-wall confinement potential, the interface conditions of the wave function in (3.4) has to be specifically imposed, namely,
⎧⎪
where Γ denotes the interface between two materials, i.e., Γ = ΩInAs∩ ΩGaAs, n is an outward normal unit vector on the boundary of ΩInAs, and Γ− and Γ+are the sets of limiting points to the curve Γ from the interior and the exterior of ΩInAs, respectively. The momentum operator Π is similarly defined for the 2D case. The boundary conditions for (3.4) are
⎧⎪
where W , S, E, and N denotes the west, south, east, and north side boundaries of the domain Ω. Note that on the west side of the boundary the values of the wave function are taken to be the same for satisfying the continuity condition across W . In actual implementation this condition is replaced by taking the values of the two horizontal grid points adjacent to W as the same. Moreover, to avoid numerical over-flow due to the term 1/r in (3.6), we do not define unknowns at the grid points on W .
Note that the potential functions Vext(r) and VB(r) can be directly reduced to the (r, z) space since these functions are independent of the azimuthal coordinate.
In real space approximation, the Hartree potential VH(r) is usually calculated by solving the Poisson equation [32]
∇ · (r)∇VH(r) =− e2 4π 0
n(r). (3.10)
By cylindrical symmetry, this equation can be written as
∂2 method of separating variables and substituting a solution of the form
VH(r, θ, z) = VH(r, z) V (θ) (3.13) Obviously, by setting V (θ) = k where k is an arbitrary constant, a function VHp(r, z) = kVH(r, z) satisfying the 2D Poisson equation is a particular solution of (3.14) in view of a second order nonhomogeneous or-dinary differential equation with respect to θ. The corresponding homogeneous general solution is eikθVHh(r, z) satisfying the Laplace equation
∂2 The general solution of the nonhomogeneous equation (3,14) is therefore of the form
k
eikθVHh(r, z) + VHp(r, z) . (3.17) In 2D setting, we also have similar interface conditions for Poisson’s problem (3.10), namely,
Similarly, the boundary conditions for (3.10) are
⎧⎪
By imposing these boundary conditions to the general solution (3.17), we deduce that the particular solution VHp (r, z) is in fact a general solution of (3.14) and thus of (3.11), i.e., VHh(r, z) = 0.
Note that for atomic systems the far side boundary condition of the Hartree potential is usually taken the values obtained by using efficient multipole expan-sion techniques [33]. For our model problem, the size of the domain is 180× 180 nm2 which in comparison with that of atomic systems is quite large and hence
the zero boundary condition for the potential on the far side of the boundary is numerically feasible, see Chapter. 4 below for numerical evidence on the choice of the size.
We then use the standard finite difference method to approximate our model problem. Since the mass in (2.3) and the Land´e factor of (2.12) are energy depen-dent, the KS equation (3.4) and its interface and boundary conditions (3.8) and (3.9) will result in a system of cubic eigenvalue equations
A0+ λA1+ λ2A2+ λ3A3 x= 0, (3.20)
where the unknown eigenpair (λ, x) is an approximate solution of (εq, φq) for some q. Starting from the Schr¨odinger equation, finite difference discretization, to the coefficient matrices A0, A1, A2, and A3, a detailed derivation of a similar cubic eigenvalue system is given in Appendix (or[24]). Several Jacobi-Davidson methods are proposed and compared in [25] for solving this type of eigenvalue problems.
Analogously, the Poisson equation (3.15) with its interface and boundary con-ditions (3.18) and (3.19) leads to a system of algebraic equations
Ax = b, (3.21)
where now the unknown vector x corresponds to the approximate values of VH(r, z) at the grid points.
We briefly describe our algorithm for the implementation of the model system
in CSDFT as follows:
Algorithm 1. A self-consistent method for the current spin density functional theory.
(1) Set VH = 0, Vxc = 0, and solve (3.20) for φ(0)q (r, z) with σ =↑ and then with σ =↓ by using the cubic Jacobi-Davidson method. Set k = 0. When B = 0, the first three lowest energies correspond to n = 1 and l = 0, 1, 2.
We therefore must solve (3.20) six times. At each time, we only seek for the smallest eigenpair. As for B = 15, the first three lowest energies correspond to n = 1, 2, 3 and l = 0. We thus solve (3.20) two times. At each time, we then seek for the three smallest eigenpairs.
(2) Evaluate the electron densities n↑(r), n↓(r), n(r), and the electron energies Eq(k). If the energies converge within an error tolerance then stop. Otherwise, set k = k + 1.
(3) Solve (3.21) for the Hartree potential VH by using GMRES [34].
(4) Evaluate Vxc via (2.13) and then solve (3.20) for the next iterate φ(k)q (r, z).
Go to (2).
There are several numerical issues deserved to be elaborated due to the special formulation of the present model when compared with the existing models of multielectronic systems of QDs. The most prominent feature of the present model is the nonparabolic dispersion relation used to define the effective mass (2.3), the
Land´e factor (2.12), and the interface condition (3.8). As a result, the set of eigenvalues that interests us is embedded in the interior of the spectrum of the eigenvalue problem (3.20) which is a nonsymmetric system. Moreover, with some feasible magnetic fields, we expect to have degenerate eigenstates due to the two identical smaller dots, i.e., the eigenvalue system is defective. Instead of using deflation scheme in the JD solver [25], we extend the generalized Davidson method of Crouzeix, Philippe, and Sadkane [35] to our cubic JD method that allows us to compute several eigenpairs simultaneously and to have a block implementation of Krylov subspaces and search direction transformation. Our JD algorithm is summarized as follows:
Algorithm 2. A cubic Jacobi-Davidson method.
(1) Choose an arbitrary orthonormal matrix V := [v1,...,vn] and let K be a given integer that limits the dimension of the basis of the subspace. Here n can be taken as 3 for six electrons with spin-up and spin-down.
(2) Compute Wk = AkV and Mk = V∗Wk, for k = 0, 1, 2, 3, where the matrices Ak are given in (3.20).
(3) For j = n, ..., K, do
(3a) Compute the eigenpairs of (θ3M3+ θ2M2+ θM1+ M0)φ = 0 by solving the generalized linear eigenvalue problem
⎡
using the QZ algorithm [36].
(3b) Select the desired eigenpairs (θi, φi) with φi 2 = 1, for i = 1, ..., n.
• Orthonormalize t against V by the modified Gram-Schmidt (MGS) method.
(4) Use n Ritz vectors u1, ..., un to create a new V := M GS(u1, ..., un) and go to Step (2) for restarting.
Note that the classical approach for dealing with the nonlinear matrix equation (3.20) is to transform the equation into a generalized linear eigenvalue system with the matrix dimension of 3 times that of (3.20) and then solve the system by the Lanczos or Arnoldi method. The matrix dimension of (3.20) for the present QDM model is about 290000. The JD method described here instead solves the generalized linear eigenvalue system in Step (3a) in a much smaller subspace V . The matrix dimension of the matrices Mi, i = 0, 1, 2, 3, is about 50 in our numerical implementation. The matrix dimension of the linearized system in Step (3a) is thus about 150. The KS Hamiltonian (2.10) is based on the nonparabolic band structure approximation. If the Hamiltonian is based on the Kane’s original form, the resulting eigenvalue problem will then be of linear form but with the matrix dimension of 4 times that of (3.20). The nonparabolic approximation thus reduces computational efforts significantly at the cost of more delicate nonlinear eigenvalue systems.