Numerical simulation of three dimensional pyramid quantum dot
Tsung-Min Hwang
a, Wen-Wei Lin
b, Wei-Cheng Wang
b, Weichung Wang
c,*aDepartment of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan
bDepartment of Mathematics, National Tsing Hua University, Hsinchu 300, Taiwan
cDepartment of Applied Mathematics, National University of Kaohsiung, Kaohsiung 811, Taiwan Received 30 May 2003; received in revised form 27 October 2003; accepted 28 October 2003
Abstract
We present a simple and efficient numerical method for the simulation of the three-dimensional pyramid quantum dot heterostructure. The pyramid-shaped quantum dot is placed in a computational box with uniform mesh in Cartesian coordinates. The corresponding Schr€odinger equation is discretized using the finite volume method and the interface conditions are incorporated into the discretization scheme without explicitly enforcing them. The resulting matrix eigenvalue problem is then solved using a Jacobi–Davidson based method. Both linear and non-linear eigenvalue problems are simulated. The scheme is 2nd order accurate and converges extremely fast. The superior performance is a combined effect of the uniform spacing of the grids and the nice structure of the resulting matrices. We have successfully simulated a variety of test problems, including a quintic polynomial eigenvalue problem with more than 32 million variables.
2003 Elsevier Inc. All rights reserved.
Keywords: Finite volume method; Heterostructure; Large scale polynomial eigenvalue problem; Semiconductor pyramid quantum dot;
Schr€odinger equation
1. Introduction
The purpose of this paper is to compute all the relevant energy states (eigenvalues) and the corre- sponding wave functions (eigenvectors) of a three-dimension (3D) semiconductor pyramidal quantum dot (QD) for electrons in the conduction band. As shown in Fig. 1, the dot is embedded in a cuboid matrix of different material. A typical example is an InAs pyramid QD embedded in a cuboid GaAs matrix.
www.elsevier.com/locate/jcp
*Corresponding author. Tel.: +886-7-5919521; fax: +886-7-5919344.
E-mail addresses: min@math.ntnu.edu.tw (T.-M. Hwang), wwlin@am.nthu.edu.tw (W.-W. Lin), wangwc@math.nthu.edu.tw (W.-C. Wang), wwang@nuk.edu.tw (W. Wang).
0021-9991/$ - see front matter 2003 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2003.10.026
The pyramid-shaped QDs are standard product of modern semiconductor manufacturing procedures.
These nano-scale QDs have been of great interest over the past few years and have stimulated numerous research activities concerning their physical properties [9,18,22,29] and applications [6,15,16,21,28].
The governing equation for this problem is the Schr€odinger equation
r h2 2mðr; kÞru
þ V ðrÞu ¼ ku; ð1Þ
where h is the reduced Plank constant, k is the unknown eigenvalue and uðrÞ the corresponding eigen- function. Here the effective electron mass mðr; kÞ and the confinement potential V ðrÞ are discontinuous across the heterojunction.
The dependence of mðr; kÞ on k can be derived from the eight-band k p analysis and the effective mass theory. This is the simplest possible model taking into account the spin–orbit split-off coupling.
First of all, starting from the Hamiltonian of the electron with a nearly periodic potential VpðrÞ which varies on a small scale (the lattice constant), we apply the L€owdin perturbation expansion up to second order to get an 8 8 Hamiltonian He that acts on the envelope functions for electrons in the lowest conduction band and the top valence band. The corresponding potential term VðrÞ is the local average of VpðrÞ. Both V ðrÞ and the envelope functions vary on a scale much larger than the lattice constant.
Secondly, instead of solving the eigenvalue problem for the Hamiltonian He, the effective mass theory proposes to project the 8 8 Hamiltonian onto the conduction band and results in a single Hamiltonian eigenvalue problem (1) with the effective mass mðkÞ given by the secular equation
detðHe kÞ ¼ 0: ð2Þ
With properly chosen basis functions, the effective Hamiltonian He can be block-diagonalized and the solution to the secular equation (2) gives the effective mass as a rational function of the energy:
mðr; kÞ ¼ m1ðkÞ in material 1ðthe dotÞ;
m2ðkÞ in material 2ðthe matrixÞ;
VðrÞ ¼ V1 in the dot;
V2 in the matrix;
ð3Þ
H
W W
Pyramid Dot Cuboid Matrix
Fig. 1. Structure schema of a pyramid quantum dot.
where 1 m‘ðkÞ¼P‘2
h2
2 kþ g‘ V‘
þ 1
kþ g‘ V‘þ d‘
; ‘¼ 1; 2: ð4Þ
Here P‘, g‘and d‘are the momentum matrix element, the conduction and spin–orbit split-off band gaps for the dot (‘¼ 1) and the matrix (‘ ¼ 2), respectively.
In the small wave number limit, mðkÞ can be further approximated by the effective mass at the Brillouin zone center. Therefore mðr; kÞ and V ðrÞ are both piecewise constant functions of r:
mðr; kÞ ¼ m1 in the dot;
m2 in the matrix;
VðrÞ ¼ V1 in the dot;
V2 in the matrix:
ð5Þ
The approximation (5) gives a quadratic dispersion relation between the energy and the wave number. Thus (5) is sometimes referred to as the parabolic approximation and (4) the non-parabolic approximation.
A detail introduction of the k p theory and the Kane model mentioned above can be found, for example, in [3,8].
An improvement over the non-parabolic approximation (4) is to take into account the effect of strain.
The strain tensor is obtained using linear elastic theory via minimization of the free energy for the structure [33]. The computed strain is then added to the effective Hamiltonian Hefollowing the Bir–Pikus theory [4].
The presence of the strain Hamiltonian changes the dispersion relation and the band gaps of the materials.
Since the strain is anisotropic and position dependent, the resulting effective mass is a tensor and the corresponding Schr€odinger equation becomes
h2 2
X
a;b¼x;y;z
@a
1
mabðr; kÞ@bþ V ðrÞ
!
u¼ ku: ð6Þ
A detailed derivation starting from the six-band Kane model with strain (neglecting the spin–orbit split- off coupling) can be found in [14]. Following the notations in [12,13], the perpendicular and in-plane components of the (diagonal) effective mass tensor are given by
1
mxxðr; kÞ¼ 1
myyðr; kÞ¼Eg m
0:25 k VlhðrÞ
þ 0:75 k VhhðrÞ
1
mzzðr; kÞ¼Eg m
1 k VlhðrÞ
;
ð7Þ
where mand Eg denote the bulk effective mass and band gap of InAs, and VhhðrÞ and VlhðrÞ the position dependent heavy hole and light-hole bands of the structure.
The potential term VðrÞ in (6) include the contributions from externally applied voltage, the conduction- band offset, the conduction-band strain potential, and the piezoelectric potential. The contribution from the strain tensor is also included in the band edge energies VcðrÞ, VhhðrÞ and VlhðrÞ [12,14].
A constant effective mass model derived from the six-band model with strain (excluding the conduction band) can be found in [8].
If in addition, the diffusion of Indium is taken into account, the material constants for InxGa1xAs can be taken as, for example, the interpolation of those of InAs and GaAs [8,12] and the effective mass tensor remains a rational function of the form (7). An explicit form of the Indium mole fraction distribution xðrÞ was proposed in [12].
In summary, the effective mass (3) and (7) can be expressed as a rational function in k of the form 1
mðr; kÞ¼Xn
i
ciðrÞ
k WiðrÞ; ð8Þ
where ciðrÞ, WiðrÞ and the potential term V ðrÞ are position dependent functions which are discontinuous across the heterojunction. This motivates us to propose an efficient numerical method for solving the generalized eigenvalue problem (1). Our method applies to the general case (8). For simplicity of presen- tation, we will only consider (4) and (5) throughout this paper.
Due to the complicated nature of the Schr€odinger equation (1), the analytic solutions can only be obtained for 1D (quantum wells) and 2D (quantum wires) models with symmetry. Relatively few re- search reports are focused on real 3D QDs. In such cases, neither analytical techniques nor asymptotic analysis gives useful information and numerical simulation becomes a very important access to this problem.
While numerical simulations play an important role [32,36,44], most existing methods are designed for 1D and 2D problems with symmetry where the matrix size of the eigenvalue problem is much smaller.
Besides, very few results can be applied by current computational methods for real 3D QDs [17, Section 11.6].
There are a few recent progress on various models of genuine 3D QD simulations. Pryor [33] studied the electronic structure of pyramidal shaped QDs. An eight-band strain-dependent k p Hamiltonian was discretized by finite differences and then solved by the Lanczos algorithm. El-Moghraby et al. [10] used a finite difference technique to solve the Schr€odinger equation of a pyramid QD with constant effective mass model. Instead of solving the eigenvalue problems, the energies (eigenvalues) were scanned over a range and each of the chosen energies are substituted into the discretized equation to form a linear system and de- termine if the corresponding matrix is singular. All of these systems are examined to obtain the candidates of the wave function (eigenvector).
In this paper, we propose to solve the linear and non-linear eigenvalue problems directly using a finite volume discretization and a Jacobi–Davidson based eigensolver.
Since the 3D problem results in a huge matrix eigenvalue problem, the discretization of the Schr€odinger equation much affects the performance of the overall scheme. From our experience in the simulation of cylindrical QDs, first order discretization gives relatively low accuracy, especially for non-linear eigenvalue problems [26,43]. As a result, local mesh refinement near the heterojunction is necessary in order to achieve reasonable accuracy. However, local mesh refinement also destroys the structure of the matrices. The di- agonal entries corresponding to the fine meshes are much larger than others. The variations among the diagonal entries made efficient preconditioning quite difficult. This is even more challenging for large scale non-linear problems.
Fortunately, thanks to the geometric structure of the pyramid, the QD can be embedded in a Cartesian coordinate with uniform mesh. We will show that, with a carefully chosen grids, the Schr€odinger equation can be discretized using a finite volume method that automatically builds in the interface condition and maintain global second order accuracy.
The discretized Schr€odinger equation then leads to the matrix polynomial eigenvalue problem AðkÞF Xs
i¼0
kiAi
!
F¼ 0; ð9Þ
where k2 C, F 2 CN, s is the degree of the matrix polynomial, Ai2 RNN, and N is the total number of unknowns. The matrix polynomial eigenvalue problem (9) is a linear (s¼ 1) eigenvalue problem for the constant effective mass model (5). It becomes a quintic or higher order polynomial eigenvalue problem (s P 5) provided the models (4) or (8) is used. The main computational challenge is that only several smallest positive eigenvalues of a very large-sized polynomial eigenvalue problem are relevant.
Traditional linear eigenvalue problem solvers, such as rational Krylov subspace methods [24,34,35], a linear system is solved approximately in each of the iterative steps. This process usually leads to intractable computational cost and sometimes inaccurate eigenpairs for large eigenvalue problems that typically arises in 3D problems. In contrast, recently developed Jacobi–Davidson methods [30,37–39] suggest a promising approach for linear and quadratic eigenvalue problems of very large size. Instead of inverting a large sparse matrix, the main source of computational cost for this approach is the matrix-vector multiplication. Only a much smaller (usually tens of unknowns) eigenvalue problem is actually solved. Inspired by the results reported in [30,37–39], we propose a modified version of the Jacobi–Davidson methods with locking and restarting deflation technique to solve the desired eigenpairs of (9).
Overall, we find the resulting spatial discretization to be second order accurate and the eigensolver converges extremely fast. The uniform placement of the grids not only reduces the programming labor, it is also crucial in accelerating the convergence of the Jacobi–Davidson step. From our simulation using 64 64 48 grids, the computed eigenvalues have only 0.2% to 0.2% relative error from its limiting value (h! 0). On the other hand, since the matrices are sparse, we are able to increase the matrix size to the machine memory limit (about 8 gigabytes) that corresponds to more than 32 million unknowns, the re- sulting eigenvalues and eigenvectors converges within 74–220 iterations, or 4800–8200 s on a decent workstation.
The rest of the paper is organized as follows. We first derive the finite volume scheme over a 2D tri- angular and a 3D pyramid domain in Sections 2.1 and 3. The Jacobi–Davidson type method for solving the corresponding eigenvalue problems is then described in Section 4. Numerical examples are reported and analyzed in Section 5 to explore the performance of the scheme. Finally we close the paper with a few concluding remarks in Section 6.
2. The discretization of the Schr€odinger equation
The Schr€odinger equations (1) and (10) with discontinuous effective mass is also known as the elliptic interface problem. Associated with the discontinuity in m, we have the following interface condition (usually refereed to as the Ben Daniel–Duke condition [3])
1 mðr; kÞ
ou on
oD
¼ 1
mðr; kÞ ou on
oD
þ
; ð10Þ
where D is the domain of the pyramid dot and n is the outward normal.
The interface problem, in general, may be solved by other second order methods such as the immersed interface method (IIM) [25], body fitting finite element method [5,7] or the immersed finite element method [11]. However, most of these methods are designed for piecewise constant coefficient problems and become quite complicated, sometimes inaccessible, for polynomial eigenvalue problems. There is also a class of finite difference schemes that replace the discontinuous coefficient by a smooth coefficient using an aver- aging process and then discretize the modified problem with standard centered difference. These scheme are known to have only first order accuracy and are thus not recommended for this problem.
Due to the special structure of the pyramid, we can embed the QD in a Cartesian coordinate with uniform mesh. We will show that the corresponding finite volume discretization results in a monotone operator with OðhÞ local truncation error on the interface and Oðh2Þ elsewhere. Therefore global second order accuracy is expected. This is indeed verified in our numerical examples. See Section 5.
For simplicity of presentation, we will give detail derivation of the scheme over 2D triangular domain and generalize the result to 3D pyramid domain.
2.1. 2D finite volume scheme for a triangular domain
Consider a 2D triangular QD with height H and base width W placed in rectangular domain. In practice, we are only interested in the confined eigenfunctions of (1) which decay exponentially outside the QD. We therefore choose a 2W 3H rectangle and impose the homogeneous Dirichlet boundary conditions u ¼ 0 on the outer boundary of the rectangle for simplicity. We fix the ratio of the mesh sizes Dx=Dy¼ W =2H as mesh refines, see Fig. 2. For convenience, the nodes in the exterior and interior of the triangle are labeled by 0 and 1, respectively. Those on the left, right, and bottom of the triangle are labeled by 2, 3, and 4; the top, left, and right corners are labeled by 5, 6, and 7, respectively.
For simplicity, we use the notation aþ h2
2m2
and a h2 2m1
ð11Þ to rewrite the Schr€odinger equation (1) and the interface condition (10) as
r ðaruÞ þ Vu ¼ ku; ð12Þ
where
a¼ a inside;
aþ outside;
V ¼ V ¼ V1 inside;
Vþ ¼ V2 outside
and aou
on
oD
¼ aþou on
oD
þ
; ð13Þ
respectively.
–2 –1.5 –1 –0.5 0 0.5 1 1.5 2
–1 –0.5 0 0.5 1 1.5 2
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 2 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 2 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 1 2 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 1 1 2 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 1 1 1 5 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 1 1 3 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 1 3 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 1 3 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 1 3 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 1 3 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 4 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Fig. 2. Schema of the uniform discretization scheme of a 2D domain.
Consider the standard centered difference method for points in the interior and the exterior
a uiþ1;j 2ui;jþ ui1;j
ðDxÞ2 þui;jþ1 2ui;jþ ui;j1
ðDyÞ2
!
þ Vui;j¼ kui;j: ð14Þ
Since all five points are on the same side of the heterojunction, a and V are therefore a constants among all five points and the local truncation error for (14) is Oðh2Þ from standard analysis.
The standard finite difference scheme (14) can also be rearranged as
a uiþ1;j ui;j
Dx Dy
þui1;j ui;j
Dx Dyþui;jþ1 ui;j
Dy Dxþui;j1 ui;j
Dy Dx
¼ Dx Dyðk V Þui;j; ð15Þ
which is a 2nd order approximation of the integral form of the Schr€odinger equation (1) over the control volume X of size Dx Dy
Z
oX
aou on¼
Z Z
X
ðk V Þu ð16Þ
as shown in (e) of Fig. 3.
This finite volume interpretation for (14) is the basis for the discretization of nodes on the heterojunction (labeled 2 through 7).
For example, consider the mesh points located at the left hypotenuse and the integral equation (16) over the control volume X as shown in (d) of Fig. 3. We denote the four portions ofoX as oXE,oXN,oXW, and oXS, respectively.
p q r
Ω p
q r
Ω Ω
Ω Ω
(a) (b) (c)
(d) (e)
Fig. 3. Discretization schema of the (a) top, (b) left, (c) right corners and (d) left hypotenuse. The solid lines represents the hetero- junctions. The solid points are the mesh points.
The line integral over, say,oXE can be approximated by Z
oXE
aunds¼ Z
oXE
auxds¼ au0iþ1
2;jDyþ OðDy3Þ ¼ auiþ1;j ui;j
Dx Dyþ Oðh3Þ;
where h¼ maxðDx; DyÞ and we can therefore approximate the left hand side of Eq. (16) by
Z
oX
aou on ¼
Z
oXE
aou on
Z
oXN
aou on
Z
oXW
aou on
Z
oXS
aou on
¼ auiþ1;j ui;j
Dx Dy aþui;jþ1 ui;j
Dy Dx aþui1;j ui;j
Dx Dy aui;j1 ui;j
Dy Dxþ Oðh3Þ:
ð17Þ On the other hand, the right hand side of Eq. (16) can be approximated by
Z Z
X
ðk V Þu ¼ Dx Dy
1
2V
þ1 2Vþ
ui;jþ kui;jþ OðhÞ
: ð18Þ
It follows that the finite volume discretization on the left hypotenuse is given by
1
ðDxÞ2aþui1;j
ðaþþ aÞui;jþ auiþ1;j
1
ðDyÞ2aui;j1
ðaþþ aÞui;jþ aþui;jþ1
þ 1
2V
þ1 2Vþ
ui;j¼ kui;jþ OðhÞ: ð19Þ
Similarly, we can derive the finite volume discretization on the right hypotenuse
1
ðDxÞ2aui1;j
ðaþþ aÞui;jþ aþuiþ1;j
1
ðDyÞ2aui;j1
ðaþþ aÞui;jþ aþui;jþ1
þ 1
2V
þ1 2Vþ
ui;j¼ kui;jþ OðhÞ ð20Þ
and on the top
1
ðDxÞ2aþui1;j
2aþui;jþ aþuiþ1;j
1
ðDyÞ2aui;j1
ðaþþ aÞui;jþ aþui;jþ1
þ 1
4V
þ3 4Vþ
ui;j¼ kui;jþ OðhÞ; ð21Þ
respectively.
The derivation for the finite volume approximation on the left and the right corners is similar. We divide oXE into two segments rq and qp (Fig. 3(b)) and then use the local expansions
uxðpÞ ¼ uxðqÞ þDy
2 uyxðqÞ þ Oðh2Þ;
uxðrÞ ¼ uþxðqÞ Dy
2 uþyxðqÞ þ Oðh2Þ
ð22Þ
to compute the east line integral. Here the superscript denotes the limiting value from inside ()) and outside (+), respectively.
From (22), we have Z
oXE
aou ond‘y¼
Z
oXE
auxd‘y¼ Z p
q
auxd‘yþ Z q
r
aþuxd‘y
¼ a uxðpÞ þ uxðqÞ 2
Dy
2 þ aþ uxðqÞ þ uxðrÞ 2
Dy
2 þ Oðh3Þ ðtrapezoidal ruleÞ
¼ a uxðqÞ
þuxðpÞ uxðqÞ 2
Dy
2 þ aþ uþxðqÞ
þuxðrÞ uþxðqÞ 2
Dy
2 þ Oðh3Þ
¼ að þ aþÞDy
2 uxðqÞ þ Dy 2
2
1 2
auyxðqÞ
aþuþyxðqÞ
þ Oðh3Þ
¼ðaþ aþÞDy
2 uxðqÞ þ Oðh3Þ
¼ðaþ aþÞ 2
uiþ1;j ui;j
Dx Dyþ Oðh3Þ;
where we have used the formula aþuþyxðqÞ ¼ auyxðqÞ;
which is a direct consequence of the interface condition aþuþyðqÞ ¼ auyðqÞ:
We can now approximate the left hand side of Eq. (16) at the left corner as
Z
oX
aou on ¼
Z
oXN
aou on
Z
oXW
aou on
Z
oXS
aou on
Z
oXE
aou on
¼ aþui;jþ1 ui;j
Dy Dx aþui1;j ui;j
Dx Dy aþui;j1 ui;j
Dy Dx
ðaþ aþÞ 2
uiþ1;j ui;j
Dx Dyþ Oðh3Þ: ð23Þ
Similarly, the right hand side of the equation is approximated by Z Z
X
ðk V Þu ¼ Dx Dy kui;j
1
8V
þ7 8Vþ
ui;jþ OðhÞ
ð24Þ
and we conclude with the finite volume discretization at the left corner as
1
ðDxÞ2 aþui1;j
aþ
þaþþ a 2
ui;jþaþþ a 2 uiþ1;j
1
ðDyÞ2aþui;j1
2aþui;jþ aþui;jþ1
þ 1
8V
þ7 8Vþ
ui;j¼ kui;jþ OðhÞ: ð25Þ
The derivation for the discretization at the right corner and the bottom points are quite similar. The results are given by (the right corner)
1 ðDxÞ2
aþ aþ 2 ui1;j
aþþ a 2
þ aþ
ui;jþ aþuiþ1;j
1
ðDyÞ2aþui;j1
2aþui;jþ aþui;jþ1
þ 1
8V
þ7 8Vþ
ui;j¼ kui;jþ OðhÞ ð26Þ
and (the bottom)
1
ðDxÞ2
aþ aþ 2 ui1;j
2ui;jþ uiþ1;j
1
ðDyÞ2aþui;j1
ðaþþ aÞui;jþ aui;jþ1
þ 1
2V
þ1 2Vþ
ui;j¼ kui;jþ OðhÞ; ð27Þ
respectively.
Remarks.
1. We can summarize the discretization (14), (19)–(21) and (25)–(27) as
rh ðarhuÞ þ V u¼ ku; ð28Þ
where aand V denote line average of a and area average of V over the control cell, respectively.
2. The finite volume scheme presented here can be derived as a finite difference scheme which incorporates the jump condition into the discretization without explicitly enforcing them. The finite difference inter- pretation of the scheme allows us the to extend the current scheme into higher order finite difference scheme. The detail derivation is given in Appendix A.
3. 3D scheme for the pyramid quantum dot
The finite volume scheme over the 3D domain can be easily extended from the 2D scheme presented in Section 2, we omit the detail derivation and only summarize them as
rh ðarhuÞ þV u ¼ ku; ð29Þ
where aandV denote the surface averages of a and the volume average of V over the controlled volume element, respectively.
We list the detailed formulas for representative points on the interior/exterior, surfaces, edges and corners as below:
• Points in the exterior of the pyramid:
1
ðDxÞ2 aþui1;j;k
2aþuijkþ aþuiþ1;j;k
þ 1
ðDyÞ2 aþui;j1;k
2aþuijkþ aþui;jþ1;k
þ 1
ðDzÞ2 aþui;j;k1
2aþuijkþ aþui;j;kþ1
¼ Vþuijk kuijk: ð30Þ
• Points in the interior of the pyramid:
1
ðDxÞ2 aui1;j;k
2auijkþ auiþ1;j;k
þ 1
ðDyÞ2 aui;j1;k
2auijkþ aui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
2auijkþ aui;j;kþ1
¼ Vuijk kuijk: ð31Þ
• Points on the western surface of the pyramid:
1
ðDxÞ2 aþui1;j;k
að þþ aÞuijkþ auiþ1;j;k
þ 1
ðDyÞ2
aþþ a
2 ui;j1;k 2uijkþ ui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
að þþ aÞuijkþ aþui;j;kþ1
¼Vþþ V
2 uijk kuijk: ð32Þ
• Points on the southern surface of the pyramid:
1 ðDxÞ2
aþþ a
2 ui1;j;k 2uijkþ uiþ1;j;k
þ 1
ðDyÞ2 aþui;j1;k
að þþ aÞuijkþ aui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
að þþ aÞuijkþ aþui;j;kþ1
¼Vþþ V
2 uijk kuijk: ð33Þ
• Points on the bottom surface of the pyramid:
1 ðDxÞ2
aþþ a
2 ui1;j;k 2uijkþ uiþ1;j;k
þ 1
ðDyÞ2
aþþ a
2 ui;j1;k 2uijkþ ui;jþ1;k
þ 1
ðDzÞ2 aþui;j;k1
að þþ aÞuijkþ aui;j;kþ1
¼Vþþ V
2 uijk kuijk: ð34Þ
• Points on the southwestern edge of the pyramid:
1
ðDxÞ2 aþui1;j;k aþþaþþ a 2
uijkþaþþ a 2 uiþ1;j;k
þ 1
ðDyÞ2
aþui;j1;k aþþaþþ a 2
uijkþaþþ a 2 ui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
að þþ aÞuijkþ aþui;j;kþ1
¼2Vþþ V
3 uijk kuijk: ð35Þ
• Points on the western edge at the bottom of the pyramid:
1
ðDxÞ2 aþui1;j;k aþþaþþ a 2
uijkþaþþ a 2 uiþ1;j;k
þ 1
ðDyÞ2
7aþþ a
8 ui;j1;k 2uijkþ ui;jþ1;k
þ 1
ðDzÞ2 aþui;j;k1
2aþuijkþ aþui;j;kþ1
¼7Vþþ V
8 uijk kuijk: ð36Þ
• Points on the southern edge at the bottom of the pyramid:
1 ðDxÞ2
7aþþ a
8 ui1;j;k 2uijkþ uiþ1;j;k
þ 1
ðDyÞ2 aþui;j1;k aþþaþþ a 2
uijk
þaþþ a 2 ui;jþ1;k
þ 1
ðDzÞ2 aþui;j;k1
2aþuijkþ aþui;j;kþ1
¼7Vþþ V
8 uijk kuijk: ð37Þ
• The southwestern corner at the bottom of the pyramid:
1
ðDxÞ2 aþui1;j;k aþþ7aþþ a 8
uijkþ7aþþ a 8 uiþ1;j;k
þ 1
ðDyÞ2
aþui;j1;k aþþ7aþþ a 8
uijkþ7aþþ a 8 ui;jþ1;k
þ 1
ðDzÞ2 aþui;j;k1
2aþuijkþ aþui;j;kþ1
¼23Vþþ V
24 uijk kuijk: ð38Þ
• The tip of the pyramid:
1
ðDxÞ2 aþui1;j;k
2aþuijkþ aþuiþ1;j;k
þ 1
ðDyÞ2 aþui;j1;k
2aþuijkþ aþui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
að þþ aÞuijkþ aþui;j;kþ1
¼5Vþþ V
6 uijk kuijk: ð39Þ
Remark.
1. The finite volume discretization (28) or (29) can be applied, in a straightforward manner, to arrays of QDs using Cartesian coordinates. It also extends easily to 2D or 3D truncated pyramid dots [12]. For ex- ample, the discretization on the left top corner of a 2D truncated dot (point P in Fig. 4(a)) is given by
P P
(a) (b)
Fig. 4. Structure schemes of truncated pyramid quantum dots in 2D and 3D.
1
ðDxÞ2 aþui1;j aþþaþþ a 2
ui;j
þaþþ a 2 uiþ1;j
1
ðDyÞ2aui;j1
ðaþ aþÞui;jþ aþui;jþ1 þ 3
8V
þ5 8Vþ
ui;j¼ kui;jþ OðhÞ:
The discretization on the southwestern corner at the top of a truncated 3D pyramid (point P in Fig. 4(b)) is given by
1
ðDxÞ2 aþui1;j;k aþþ5aþþ 3a 8
uijk
þ5aþþ 3a 8 uiþ1;j;k
þ 1
ðDyÞ2 aþui;j1;k aþþ5aþþ 3a 8
uijk
þ5aþþ 3a 8 ui;jþ1;k
þ 1
ðDzÞ2 aui;j;k1
ðaþ aþÞuijkþ aþui;j;kþ1
¼17Vþþ 7V
24 uijk kuijk:
2. For QDs of general shape, e.g. lens shaped ones, a more sophisticated discretization method is needed.
We have adopted the finite difference method on a body-fitting curvilinear coordinate system developed in [20]. The test results are quite encouraging and the details will be reported elsewhere.
4. Jacobi–Davidson based algorithm for eigenvalue problems
The finite volume discretization (29) results in various eigenvalue problems. For the constant coefficient case, a and V are piecewise constants, we have the linear eigenvalue problem
A0F¼ kA1F; ð40Þ
where A1is the identity matrix and A0is a sparse matrix with non-zero entries located in the main diagonal and six off-diagonals. Furthermore, since the matrix A0is symmetric and positive definite, the eigenvalues of (40) are all positive. In other words, the desired eigenvalues are the ones located in the lower end of the matrix spectrum.
For the non-parabolic model (4), VðrÞ is piecewise constant and a is a rational function of the k with piecewise constant coefficients. At each grid point, the discretization (29) results in a rational function of k.
By multiplying the common denominators of the rational functions ðk þ g1 V1Þðk þ g1 V1þ d1Þ
and
ðk þ g2 V2Þðk þ g2 V2þ d2Þ
for the grids in the interior and exterior, respectively, Eqs. (30) and (31) become cubic functions of k. For all interface grid points (32)–(39), we multiply (29) by the common denominator
ðk þ g1 V1Þðk þ g1 V1þ d1Þðk þ g2 V2Þðk þ g2 V2þ d2Þ:
The result is a quintic polynomial of k. With elementary algebraic manipulation, the discretization (29) is reduced to a quintic polynomial eigenvalue problem:
ðk5A5þ k4A4þ k3A3þ k2A2þ kA1þ A0ÞF ¼ 0: ð41Þ
For this polynomial eigenvalue problem, the target smallest positive eigenvalues are embedded in the interior of the spectrum of the matrix polynomial (41). This is merely an algebraic artifact as we have approximated the effective mass using rational functions of the unknown eigenvalues, hence raised the degree of the polynomial and produced extra roots located on the negative real axis. These negative ei- genvalues are non-physical and should not be taken as solutions.
To compute the smallest positive eigenvalue for the linear eigenvalue problem (40), variants of the rational Krylov subspace method [24,34,35] can be used to compute the eigenpairs. Theoretically, these methods can also be applied to the quintic polynomial (41) written as an enlarged linear eigenvalue problem
0 I 0 0 0
0 0 I 0 0
0 0 0 I 0
0 0 0 0 I
A0 A1 A2 A3 A4
2 66 66 4
3 77 77 5
F kF k2F k3F k4F 2 66 66 4
3 77 77 5¼ k
I 0 0 0 0
0 I 0 0 0
0 0 I 0 0
0 0 0 I 0
0 0 0 0 A5
2 66 66 4
3 77 77 5
F kF k2F k3F k4F 2 66 66 4
3 77 77
5: ð42Þ
However, the relevant eigenvalues are the smallest positive ones and they are located in the interior of the matrix eigenvalue spectrum. To solve them successively by these methods, the classical shift-and-invert technique is much too expensive for large scale problems. Moreover, the enlarged matrix is five times larger than the AiÕs. and does not have a good structure as the AiÕs do. The condition number for the enlarged matrix can be significantly larger than those of the AiÕs since it has a large set of admissible perturbations [42]. As a consequence, all these approach are not favorable in terms of accuracy and efficiency. To bypass these difficulties, it is desirable to have numerical methods that deal with the eigenvalue problems (40) and (41) directly by such as the Jacobi–Davidson methods. It is thus straightforwardly to use the linear Jacobi–
Davidson method [38] for the linear eigenvalue problem (40). However, rare numerical results were reported for the quintic eigenvalue problem (41), while variants of the Jacobi–Davidson methods are considered for the generalized linear [37], quadratic [19,39], and cubic [43] eigenvalue problems.
Besides, to compute successively all other desired eigenvalues, deflation or locking schemes are neces- sary. The implicit deflation techniques based on the Schur forms is known to have good performance for linear eigenvalue problems [1, Section 4.7 and 8.4]. However, it is not clear how to incorporate the implicit deflation technique in a quintic polynomial eigenvalue problem since the Schur form is not defined for a quintic polynomial matrix in general. Alternatively, an explicit deflation scheme is proposed in [43] for cubic eigenvalue problems. Once the smallest eigenvalue is obtained, the scheme transforms the computed eigenvalue to infinity and keeps all other eigenvalues unchanged. The second eigenvalue thus becomes the smallest eigenvalue of the newly transformed eigenvalue problem, which can be solved by the cubic Jacobi–
Davidson method. The technique is then applied repeatedly until all desired eigenvalues are determined. In addition to the deflation schemes, Meerbergen [30] proposed a quadratic eigensolver using the locking and restarting scheme based on the Schur form of the linearized problem, which is similar to (41). The method gives a linkage between the methods for solving the quadratic and the linearized eigenvalue problems.
Here we propose a general polynomial Jacobi–Davidson algorithm, shown in Fig. 6, to compute all the desired eigenvalues for the problem (9). This algorithm is based on a polynomial Jacobi–Davidson method modified directly from the quadratic Jacobi–Davidson method illustrated in [2, Section 9.2] and a locking technique similar to the one mentioned in [30]. This algorithm is capable of calculating the smallest positive eigenvalue first and then computing successively all other desired eigenvalues by suitably choosing the orthonormal searching space V ¼ ½VF; Vini:
• The algorithm computes the smallest positive eigenvalue k1and the associated eigenvector F1by the gen- eral Jacobi–Davidson method first. In this stage, VF is empty, Viniis any matrix satisfying ViniTVini¼ I (in Step (1) of Fig. 6).
• To compute the second smallest positive eigenvalue k2 (j¼ 2 in Step (2)), the algorithm normalizes the approximate eigenvalue F1 and adds it to the initial search space VF. That is, the algorithm computes VF¼kFF1
1kand chooses an orthonormal matrix Vini such that½VF; ViniT½VF; Vini ¼ I. By doing so, the first eigenvalue k1 will be included (or ‘‘locked’’) in the eigenpairs computed by
ðVTAðkÞV Þs ¼ 0 ð43Þ
in Step (2.2.i). Therefore, the algorithm chooses the second smallest positive eigenvalue from the ones in (43).
• The algorithm computes k^jin a similar way. The algorithm sets VF to be an orthonormal basis of the eigenspace spanned by F1; . . . ; F^j1 and then choose an orthonormal matrix Vini satisfying
½VF; ViniT½VF; Vini ¼ I. The eigenvalues k1 through k^j1 will be locked in the eigenpairs computed in Eq.
(43). The algorithm thus chooses the ^jth smallest positive eigenvalue from the ones obtained in (43) as an approximation of k^j.
To conclude this algorithm, we note that the correction equation Ipu
up
AðhÞðI uuÞt ¼ r
ð44Þ
in Step (2.2.v) is solved approximately by computing
t¼ MA1rþ eMA1p with e¼uMA1r
uMA1p: ð45Þ
The matrix MAis a preconditioner of AðhÞ. For example, in the SSOR preconditioning scheme [40] used in our numerical experiments,
AðhÞ MA¼ ðD xLÞD1ðD xU Þ;
where x is a scalar, AðhÞ ¼ D L U with D ¼ diagðAðhÞÞ, L and U are strictly lower and upper triangular of AðhÞ.
Remark. The polynomial Jacobi–Davidson algorithm outlined in Fig. 6 can be applied to general matrix polynomial eigenvalue problems. We have implemented and tested it up to s¼ 5 at very large scale (32 million unknowns) and are quite satisfied with its performance. As analyzed in Section 4.1, the performance of our scheme is related to the structures of the matrices Ai. The algorithm did not make use of any special properties related to the degree of the polynomial and we expect it to work as well for higher degree polynomial eigenvalue problems that arises from discretizing (7) over regular or truncated pyramid dots.
4.1. Analysis of the quintic polynomial eigenvalue problem
We have found the performance of the algorithm quite satisfactory for both the linear and quintic polynomial eigenvalue problems. See Section 5 for detail. Here we provide a heuristic explanations of the rapid convergence of the algorithms based on the structures of the coefficient matrices Ai. Fig. 5 shows the sparsity patterns of the matrices A0 through A5 withðL; M; N Þ ¼ ð8; 8; 6Þ.
For the constant effective mass model (40), it is easy to verify that A0 is symmetric and diagonal dominant.
For the non-parabolic effective mass approximation (41), the matrices Aihave the following properties:
4.1.1. Symmetry
The matrices Ai, i¼ 0; 1; 2; 3, are symmetric except the entries involving the grids on the interface. A4and A5 are diagonal. See Fig. 5(c).
4.1.2. Diagonal dominance
The matrix A0 is diagonal dominant. The matrices A1, A2, and A3are almost diagonal dominant except the rows involving the grids on the interface. These rows of A2and A3are indicated in part (d) of Fig. 5. The matrices A4 and A5 are diagonal dominant since they are diagonal matrices.
4.1.3. Asymptotic structure
• The structures of A0and A1 mentioned above are independent of h.
• The entries of A2 and A3 are of Oðh2Þ except those associated with the interface grids. For the matrices A2 and A3, we observe the main diagonal of the matrices are formed by Oðh2Þ entries, except the ones associated with the interface conditions. The rows associated with these interface grids has non-zero Oð1Þ entries on the main diagonal and six off-diagonals. The sparsity of A2 and A3 are shown in Fig. 5(b).
• The entries of the diagonal matrices A4 and A5 are all Oðh2Þ.
0 50 100 150 200
0
50
100
150
200
nz = 1477 Sparsity of A
0 and A
1
0 50 100 150 200
0
50
100
150
200
nz = 449 Sparsity of A
2 and A
3
0 50 100 150 200
0
50
100
150
200
nz = 34 Sparsity of A
4 and A
5
0 50 100 150 200
0
50
100
150
200
nz = 238
Sparsity of the entries involving the interfaces (a)
(c)
(b)
(d)
Fig. 5. Sparsity patterns of the matrices A0through A5are presented in part (a)–(c). In part (d), only the entries involving the interfaces are plotted. The notation ‘‘nz’’ in the figure stands for the number of non-zero entries in the matrices.
In summary, as h goes to zero, the quintic matrix polynomial can be viewed as a small perturbation of the generalized linear eigenvalue problem
A0F¼ kA1Fþ Oðh2Þ þ B;
with A0 being diagonal dominant and nearly symmetric. Here B is a low rank matrix. The only non-zero entries are those associated with interface grids.
Fig. 6. The polynomial Jacobi–Davidson algorithm designed to compute the r smallest positive eigenvalues and the associated eigenvectors of polynomial eigenvalue problems (9).
