Jacobi Davidson method and its applications

(1)

師大

Jacobi Davidson method and its applications

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University, Taiwan

October 31, 2008

T.M. Huang (Taiwan Normal Univ.) JD method and applications October 31, 2008 1 / 71

(2)

師大

Outline

1 Some basic theorems

2 Jacobi-Davidson method for linear eigenvalue problems

3 Polynomial eigenvalue problems

4 Locking method for computing λ2, . . . , λ`

5 Non-equivalence deflation of quadratic eigenproblems

6 Non-equivalence deflation for cubic polynomial eigenproblems

7 Applications

8 Numerical experiments for linear eigenproblems

9 Numerical experiments for polynomial eigenproblems

(3)

師大

Some basic theorems

Eigenproblems

Standard eigenproblems: Ax = λx Generalized eigenproblems: Ax = λBx

Higher order poly. eigenproblems: (A0+ λA1+ .... + λⁿAn)x = 0 Eigenproblems of λ-matrices: F (λ)x = 0

What do we care ?

(i) In theory: eigenstructure, spectrum decomposition, canonical form, . . . , etc.

(ii) In computation: eigenvalues, eigenvectors, invariant subspaces, . . . , etc.

(4)

師大

Theorem (Fischer)

Let the Hermitian matrix A have eigenvalues λ₁ ≥ λ₂ ≥ · · · ≥ λ_n. Then

λi = min dim(W)=n−i +1

w ∈W,kw kmax2=1w^HAw

and

λ_i = max dim(W)=i

min

w ∈W,kw k2=1w^HAw

.

(5)

師大

Theorem

Let X be an eigenspace of A and let X be a basis for X . Then there is a unique matrix L such that

AX = XL.

The matrix L is given by

L = XÎAX , where XÎ is a matrix satisfying XÎX = I .

If (λ, x ) is an eigenpair of A with x ∈ X , then (λ, X^Ix ) is an eigenpair of L.

Conversely, if (λ, u) is an eigenpair of L, then (λ, Xu) is an eigenpair of A.

(6)

師大

Proof:

Let

X = [x₁· · · x_k] and Y = AX = [y₁· · · y_k] .

Since y_i ∈ X and X is a basis for X , there is a unique vector `_i such that y_i = X `_i.

If we set L = [`₁· · · `_k], then AX = XL and L = X^IXL = X^IAX .

Now let (λ, x ) be an eigenpair of A with x ∈ X . Then there is a unique vector u such that x = Xu. However, u = X^Ix . Hence

λx = Ax = AXu = XLu ⇒ λu = λX^Ix = Lu.

Conversely, if Lu = λu, then

A(Xu) = (AX )u = (XL)u = X (Lu) = λ(Xu), so that (λ, Xu) is an eigenpair of A.

(7)

師大

Theorem (Optimal residuals) Let [X X⊥] be unitary. Let

R = AX − XL and S^H = X^HA − LX^H. Then kRk and kS k are minimized when

L = X^HAX , in which case

(a) kRk = kX_⊥^HAX k, (b) kSk = kX^HAX⊥k, (c) X^HR = 0.

(8)

師大

Proof:

Set

X^H X_⊥^H

A

X X⊥ =

ˆL H

G M

. Then

X^H X_⊥^H

R =

ˆL H

G M

X^H X_⊥^H

X −

X^H X_⊥^H

XL =

L − Lˆ G

. It implies that

kRk =

X^H X_⊥^H

R

=

ˆL − L G

, which is minimized when L = ˆL = X^HAX and

min kRk = kG k = kX_⊥^HAX k.

The proof for S is similar. If L = X^HAX , then

X^HR = X^HAX − X^HXL = X^HAX − L = 0.

(9)

師大

Definition

Let X be of full column rank and let X^I be a left inverse of X . Then X^IAX is a Rayleigh quotient of A.

Theorem

Let X be orthonormal, A be Hermitian and R = AX − XL.

If `1, . . . , `_k are the eigenvalues of L, then there are eigenvalues λ_j₁, . . . , λ_j_k of A such that

|`_i − λ_j_i| ≤ kRk₂ and v u u t

k

X

i =1

(`_i − λ_j_i)² ≤√

2kRk_F.

(10)

師大

Jacobi’s orthogonal component correction (JOCC), 1846

Consider the eigenvalue problem A

1 z

≡

α c^T b F

1 z

= λ

1 z

, (1)

where A is diagonal dominant and α is the largest diagonal element.

(1) is equivalent to

λ = α + c^Tz, (F − λI )z = −b.

Jacobi iteration : (with z1 = 0)

θ_k = α + c^Tz_k,

(D − θ_kI )z_k+1 = (D − F )z_k − b (2) where D = diag(F ).

(11)

師大

Davidson’s method (1975)

Algorithm (Davidson’s method) Given unit vector v , set V = [v ] Iterate until convergence

Compute desired eigenpair (θ, s) ofV^TAV. Computeu = Vs andr = Au − θu.

If (k r k₂< ε), stop.

Solve (D_A− θI )t = r.

Orthog. t ⊥ V → v , V = [V , v ] end

(12)

師大

Let u_k = (1, z_k^T)^T. Then

rk = (A − θkI )uk =

α − θ_k+ c^Tz_k (F − θ_kI )z_k + b

Substituting the residual vector r_k into linear systems (DA− θ_kI )tk = −rk, where DA =

α 0

0 D

, we get

(D − θ_kI )y_k = −(F − θ_kI )z_k − b

= (D − F )z_k − (D − θ_kI )z_k− b From (2) and above equality, we see that

(D − θ_kI )(z_k+ y_k) = (D − F )z_k − b = (D − θ_kI )z_k+1

This implies that z_k+1 = z_k+ y_k as one step of JOCC starting with z_k.

(13)

師大

Polynomial eigenvalue problems

Polynomial eigenvalue problems:

A(λ)x ≡ λ^τA_τ + · · · + +λ²A₂+ λA₁+ A₀ x = 0. (3) Enlarged linear eigenvalue problem: (e.g. cubic polynomial)





0 I 0

0 0 I

A₀ A₁ A₂







 x λx λ²x



= λ





I 0 0

0 I 0

0 0 −A₃







 x λx λ²x



. Disadvantages:

I The order of the larger matrices are tripled

I The condition number of eigenvalues and eigenvectors may increase

(14)

師大

Consider the quadratic eigenvalue problem¹ Q(λ)x ≡ λ²M + λC + K x = 0 with

M = 1

2In, K = diag_{1≤j ≤n} j²π²(j²π²+ τ − κv²)/2 , C = −C^T = (c_ij) with c_ij =

( _4ij

j²−i²v if i + j is odd 0 otherwise.

Take v = 10, κ = 0.8 and τ = 77.9.

Enlarged linear eigenvalue problem:

K 0 0 I

x λx

= λ

−C −M

I 0

x λx

1F. Tisseur and K. Meerbergen, SIAM Rev. Vol. 43, pp. 235-286, 2001

(15)

師大

−1 0 1

x 10⁻⁹

−4

−3

−2

−1 0 1 2 3

4x 10⁴ n = 60

polyeig exact sol

Figure: The spectrum of Q(λ)

(16)

師大

Polynomial Jacobi-Davidson method

(θ_k, u_k): approx. eigenpair of A(λ), θ_k ≈ λ, with

u_k = V_ks_k, V_k^TA(λ)V_ks_k = 0 and k s_k k₂= 1.

Let

r_k = A(θ_k)u_k. Then

u_k^Tr_k = u^T_kA(θ_k)u_k = s_k^TV_k^TA(θ_k)V_ks_k = 0 ⇒ r_k ⊥ u_k Find the correction t such that

A(λ)(uk+ t) = 0.

That is

A(λ)t = −A(λ)uk = −rk + (A(θk) − A(λ))uk.

(17)

師大

Let

p_k = A⁰(θ_k)u_k ≡

τ

X

i =1

i θ^{i −1}_k Ai

! u_k.

A(λ) = A − λI :

pk = −uk,

(A(θ_k) − A(λ))u_k = (λ − θ_k)u_k = (θ_k − λ_k)p_k A(λ) = A − λB:

pk = −Buk,

(A(θ_k) − A(λ))u_k = (λ − θ_k)Bu_k = (θ_k− λ)p_k A(λ) =Pτ

i =0λⁱA_i with τ ≥ 2:

(A(θ_k) − A(λ))u_k =

(θ_k− λ)A⁰(θ_k) −1

2(θ_k − λ)²A⁰⁰(ξ_k)

u_k

= (θk − λ)p_k −1

2(θk − λ)²A⁰⁰(ξk)uk

(18)

師大

Hence

A(λ)t = −rk+ (θk − λ)p_k− 1

2(θk− λ)²A⁰⁰(ξk)uk

Since r_k ⊥ u_k, we have

I −p_ku_k^T u_k^Tpk

A(λ)t = −r_k −1

2(θ_k− λ)²

I − p_ku_k^T u_k^Tpk

A⁰⁰(ξ_k)u_k. Correction equation:

I − pku^T_k u_k^Tp_k

A(θ_k)(I − u_ku_k^T)t = −r_k and t ⊥ u_k,

or

I − pku_k^T u_k^Tp_k

(A − θ_kB)

I − ukp^T_k p_k^Tu_k

t = −r_k and t ⊥_B u_k,

with symmetric positive definite matrix B.

(19)

師大

Solving correction vector t

Correction equation:

I − pku^T_k u_k^Tp_k

A(θ_k)(I − u_ku_k^T)t = −r_k. (4) Method I:

Use preconditioning iterative approximations, e.g., GMRES, to solve (4).

Use a preconditioner M_p≡

I − p_ku^T_k u_k^Tpk

M

I − uku^T_k

≈

I − p_ku_k^T u_k^Tpk

A(θk)

I − uku_k^T

, where M is an approximation of A(θ_k).

In each of the iterative steps, it needs to solve the linear system

M_pt = y , t ⊥ u_k (5)

for a given y .

(20)

師大

Since t ⊥ u_k, Eq. (5) can be rewritten as

I − p_ku^T_k u_k^Tp_k

Mt = y ⇒ Mt = u_k^TMt

u_k^Tp_k p_k+ y ≡ η_kp_k + y . Hence

t = M⁻¹y + η_kM⁻¹p_k, where

ηk = − u_k^TM⁻¹y u_k^TM⁻¹p_k.

SSOR preconditioner: Let A(θ_k) = L + D + U. Then M = (D + ωL)D⁻¹(D + ωU).

(21)

師大

Method II: Since t ⊥ u_k, Eq. (4) can be rewritten as A(θk)t = u^T_kA(θ_k)t

u_k^Tp_k pk − r_k ≡ εp_k − r_k. (6) Let t1 and t2 be approximated solutions of the following linear

systems:

A(θ_k)t = −r_k and A(θ_k)t = p_k, respectively. Then the approximated solution ˜t for (6) is

˜t = t₁+ εt₂ for ε = −u_k^Tt₁ u_k^Tt2

. The approximated solution ˜t for (6) is

˜t = −M⁻¹r_k+ εM⁻¹p_k for ε = u^T_kM⁻¹r_k u_k^TM⁻¹p_k, where M is an approximation of A(θk).

(22)

師大

Method III:

Eq. (6) implies that

t = εA(θ_k)⁻¹p_k − A(θ_k)⁻¹r_k = εA(θ_k)⁻¹p_k − u_k. Let t1 be approximated solution of the following linear system:

A(θ_k)t = p_k. Then the approximated solution ˜t for (6) is

˜t = εt₁− u_k for ε =

u_k^Tt₁−1

.

(23)

師大

Algorithm (Jacobi-Davidson Algorithm for solving A(λ)x = 0) Choose an n-by-m orthonormal matrix V0

Do k = 0, 1, 2, · · ·

Compute all the eigenpairs of V_k^TA(λ)Vk = 0.

Select the desired (target) eigenpair (θ_k, s_k) with ks_kk₂= 1.

Computeu_k = V_ks_k,r_k = A(θ_k)u_k andp_k = A⁰(θ_k)u_k. If (krkk₂ < ε), λ = θk, x = uk, Stop

Solve (approximately) a t_k ⊥ u_k from (I − ^p^k^u

T k

u^T_kpk)A(θk)(I − uku^T_k)t = −rk. Orthogonalize t_k ⊥ V_k → v_k+1, V_k+1 = [V_k, v_k+1]

(24)

師大

Restarting

A double precision real variable needs to 8 bytes to save in Memory.

125 double precision real variables ≈ 1 KB 125, 000 double precision real variables ≈ 1 MB

Keep the locked Schur vectors as well as the Schur vectors of interest in the subspace and throw away those we are not interested.

(25)

師大

(v) Solve correct equation (approximately) to obtain a t ⊥ u_k by the method determined below.

If ( krkk₂> 0.1 and k ≤ 9 ) then

Use {BiCGSTAB, No precond., 7, 10⁻³} else

Use {GMRES, SSOR, 30, 10⁻³} End if

Figure: The heuristic strategy for computing the first target eigenvalue.

(26)

師大

Locking for Ax = λx

V_k with V_k^∗V = I_k are convergentSchur vectors, i.e., AV_k = V_kT_k

for someupper triangular Tk. SetV = [Vk, Vq]with V^∗V = Ik+q in k + 1-th iteration of Jacobi-Davidson Algorithm. Then

V^∗AV =

V_k^∗AV_k V_k^∗AV_q V_q^∗AVk V_q^∗AVq

=

T_k V_k^∗AV_q V_q^∗VkTk V_q^∗AVq

=

T_k V_k^∗AV_q 0 V_q^∗AVq

.

(27)

師大

Locking Polynomial Jacobi–Davidson Method (0) Given A(λ) =Pτ

i =0λⁱA_i and the number of desired eigenvalues σ.

(1) Initialize V = [V_ini] as an orthonormal matrix and V_x = [].

(2) For j = 1, 2, ..., σ

(2.1) Iterate until convergence

(i) Compute thej th desired eigenpairs (θ, s) of V^TA(θ)V . (ii) Compute u = Vs, p = A⁰(θ)u, and r = A(θ)u.

(iii) If (||r ||2 < ε), Set λj = θ, xj = u, Goto Step (2.2).

(iv) Solve (approximately) a t ⊥ u from (I − ^pu_uT^Tp)A(θ)(I − uu^T)t = −r . (v) Orthogonalize t ⊥ V → v , V = [V , v ] (2.2) Orthogonalizexj ⊥ V_x → x_j; Vx = [Vx, xj]

(2.3) Choose an orthonormal matrixVini ⊥ V_x; Set V = [Vx, Vini] Ref: G. L. G. Sleijpen, G. L. Booten, D. R. Fokkema and H. A. van der Vorst, BIT, 36:595-633, 1996

(28)

師大

(v) Solve correct equation (approximately) to obtain a t ⊥ u_k by the method determined below.

If ( k r_k k₂> 0.1 and k < 10 ) then

Use {BiCGSTAB, No precond., 7, 10⁻³} else if ( k r_k k₂≥ 0.1 and k > 14 ) then

Use {GMRES, SSOR, 30, 10⁻³}

else if ( k r_k k₂< 0.1 and k r_k−1 k₂/ k r_k k₂< 4) then Set j = min(30, j + 2) and use {GMRES, SSOR, j, 10⁻³} else

Use {GMRES, SSOR, j, 10⁻³} End if

Figure: The heuristic strategy for computing eigenvalues other than the first convergent eigenvalue.

(29)

師大

Non-equivalence deflation of quadratic eigenproblems

Let λ1 be a real eigenvalue of Q(λ) and x1, z1 be the associated right and left eigenvectors, respectively, with z₁^TKx₁ = 1. Let

θ1 = (z₁^TMx1)⁻¹. We introduce a deflated quadratic eigenproblem

Q(λ)x ≡e h

λ²M + λ ee C + eKi x = 0, where

Me = M − θ₁Mx₁z₁^TM, Ce = C + θ1

λ1

(Mx1z₁^TK + Kx1z₁^TM), Ke = K − θ1

λ²₁Kx1z₁^TK .

(30)

師大

Complex deflation

Let λ₁ = α₁+ i β₁ be a complex eigenvalue of Q(λ) and x₁ = x_1R+ ix_1I, z₁ = z_1R+ iz_1I be the associated right and left eigenvectors, respectively, such that

Z₁^TKX1 = I2, where X1 = [x1R, x1I] and Z1 = [z1R, z1I]. Let

Θ₁ = (Z₁^TMX₁)⁻¹.

Then we introduce a deflated quadratic eigenproblem with Me = M − MX₁Θ₁Z₁^TM,

Ce = C + MX₁Θ₁Λ^−T₁ Z₁^TK + KX₁Λ⁻¹₁ Θ^T₁Z₁^TM, Ke = K − KX1Λ⁻¹₁ Θ1Λ^−T₁ Z₁^TK

in which Λ₁ =

α₁ β₁

−β₁ α₁

.

(31)

師大

Theorem

(i) Let λ1 be a simple real eigenvalue of Q(λ). Then the spectrum of Q(λ) is given bye

(σ(Q(λ)){λ1}) ∪ {∞}

provided that λ²₁ 6= θ₁.

(ii) Let λ₁ be a simple complex eigenvalue of Q(λ). Then the spectrum of eQ(λ) is given by

σ(Q(λ)){λ1, ¯λ₁} ∪ {∞, ∞}

provided that Λ₁Λ^T₁ 6= Θ₁.

Furthermore, in both cases (i) and (ii), if λ₂ 6= λ₁ and (λ₂, x₂) is an eigenpair of Q(λ) then the pair (λ2, x2) is also an eigenpair of eQ(λ).

Ref: T.-M. Hwang, W.-W. Lin and V. Mehrmann, SIAM J. Sci. Comput.

(32)

師大

Suppose that M, C , K are symmetric. Given an eigenmatrix pair

(Λ1, X1) ∈ R^{r ×r} × R^n×r of Q(λ), where Λ1 is nonsingular and X1 satisfies X₁^TKX1 = Ir, Θ1:= (X₁^TMX1)⁻¹.

We define ˜Q(λ) := λ²M + λ ˜˜ C + ˜K , where M˜ := M − MX₁Θ₁X₁^TM,

C˜ := C + MX₁Θ₁Λ^−T₁ X₁^TK + KX₁Λ⁻¹₁ Θ₁X₁^TM, K˜ := K − KX1Λ⁻¹₁ Θ1Λ^−T₁ X₁^TK .

Theorem

Suppose that Θ₁− Λ₁Λ^T₁ is nonsingular. Then the eigenvalues of the real symmetric quadratic pencil ˜Q(λ) are the same as those of Q(λ) except that the eigenvalues of Λ1, which are closed under complex conjugation, are replaced by r infinities.

(33)

師大

Proof:

Since (Λ₁, X₁) is an eigenmatrix pair of Q(λ), i.e., MX₁Λ²₁+ CX₁Λ₁+ KX₁ = 0, we have

Q(λ)˜ = Q(λ) + [MX1(λIr + Λ1) + CX1] Θ1Λ^−T₁ (X₁^TK − λΛ^T₁X₁^TM)

= Q(λ) + Q(λ)X1(λIr − Λ₁)⁻¹Θ1Λ^−T₁ (X₁^TK − λΛ^T₁X₁^TM).

By using the identity

det(I_n+ RS ) = det(I_m+ SR), where R, S^T ∈ R^n×m, we have

det[ ˜Q(λ)]

= det[Q(λ)]det[I + X1(λIr − Λ₁)⁻¹Θ1Λ^−T₁ (X₁^TK − λΛ^T₁X₁^TM)]

= det[Q(λ)]det[Ir + (λIr − Λ₁)⁻¹Θ1Λ^−T₁ (Ir − λΛ^T₁Θ⁻¹₁ )]

= det[Q(λ)]

det(λIr − Λ₁)det(Θ1Λ^−T₁ − Λ₁).

(34)

師大

Since (Θ1− Λ₁Λ^T₁) ∈ R^{r ×r} is nonsingular, we have det(Θ₁Λ^−T₁ − Λ₁) 6= 0.

Therefore, ˜Q(λ) has the same eigenvalues as Q(λ) except that r eigenvalues of Λ₁ are replaced by r infinities.

(35)

師大

Non-equivalence deflation for cubic polynomial eigenproblems

Let (Λ, V_u) ∈ R^{r ×r} × R^n×r be an eigenmatrix pair of A(λ) with V_u^TVu= Ir and 0 /∈ σ(Λ), i.e.,

A₃V_uΛ³+ A₂V_uΛ²+ A₁V_uΛ + A₀V_u= 0. (7) Define a new deflated cubic eigenvalue problem by

A(λ)u = (λe ³A˜3+ λ²A˜2+ λ˜A1+ ˜A0)u = 0, (8) where











A˜0 = A0,

A˜₁ = A₁− (A₁V_uV_u^T+ A₂V_uΛV_u^T + A₃V_uΛ²V_u^T), A˜₂ = A₂− (A₂V_uV_u^T+ A₃V_uΛV_u^T),

A˜3 = A3− A₃VuV_u^T.

(9)

(36)

師大

Lemma

Let A(λ) and eA(λ) be cubic pencils given by (3) and (8), respectively.

Then it holds

A(λ) = A(λ)e

I_n− λV_u(λI_r − Λ)⁻¹V_u^T

. (10)

Theorem

Let (Λ, V_u) be an eigenmatrix pair of A(λ) with V_u^TV_u= I_r. Then (i) (σ(A(λ))σ(Λ)) ∪ {∞} = σ(eA(λ)).

(ii) Let (µ, z) be an eigenpair of A(λ) with k z k₂= 1 and µ /∈ σ(Λ).

Define

˜

z = (I_n− µV_uΛ⁻¹V_u^T)z ≡ T (µ)z. (11) Then (µ, ˜z) is an eigenpair of eA(λ).

(37)

師大

Proof of Lemma:

Using (9) and (7), and the fundamental matrix calculation, we have A(λ)e = A(λ) − λ

λ²A3VFV_F^T + λA2VFV_F^T + λA3VFΛV_F^T+ A1VFV_F^T +A2V_FΛV_F^T + A3V_FΛ²V_F^T

= A(λ) − λ

A₃V_F(λI_r − Λ)³(λI_r − Λ)⁻¹V_F^T +3A₃V_FΛ(λI_r − Λ)²(λI_r − Λ)⁻¹V_F^T +3A3V_FΛ²(λIr − Λ)(λI_r − Λ)⁻¹V_F^T +A2VF(λIr − Λ)²(λIr − Λ)⁻¹V_F^T +2A₂V_FΛ(λI_r − Λ)(λI_r − Λ)⁻¹V_F^T

+A₁V_F(λI_r − Λ)(λI_r − Λ)⁻¹V_F^T

(38)

師大

A(λ)e = A(λ) − λA₃V_F(λ³I_r − Λ³) + A₂V_F(λ²I_r − Λ²) +A1VF(λIr − Λ) + A₀VF − A₀VF] (λIr − Λ)⁻¹V_F^T

o

= A(λ) − λh

A(λ)V_F(λI_r − Λ)⁻¹V_F^Ti

= A(λ)h

I_n− λV_F(λI_r − Λ)⁻¹V_F^Ti .

(39)

師大

Proof of Theorem

: (i) Using the identity

det(I_n+ RS ) = det(I_m+ SR) and Lemma 8, we have

det( eA(λ)) = det(A(λ)) det

In− λV_F(λIr − Λ)⁻¹V_F^T

= det(A(λ)) det In− λ(λI_r − Λ)⁻¹

= det(A(λ)) det(λI_r − Λ)⁻¹det(−Λ).

Since 0 /∈ σ(Λ), det(−Λ) 6= 0. Thus, eA(λ) and A(λ) have the same finite spectrum except the eigenvalues in σ(Λ). Furthermore, dividing Eq. (8) by λ³ and using the fact that

Ae3V_F = (A3− A₃V_FV_F^T)V_F = 0,

we see that (diag_r{∞, · · · , ∞}, V_F) is an eigenmatrix pair of eA(λ) corresponding to infinite eigenvalues.

(40)

師大

(ii) Since µ /∈ σ(Λ), the matrix T (µ) = (I − µV_FΛ⁻¹V_F^T) in (11) is invertible with the inverse

T (µ)⁻¹= I_n− µV_F(µI_r − Λ)⁻¹V_F^T. (12) From Lemma 8, we have

A(µ)˜e z = A(µ) h

In− µV_F(µIr − Λ)⁻¹V_F^T i h

In− µV_FΛ⁻¹V_F^T i

z = 0.

This completes the proof.

(41)

師大

Applications

Quantum well (2 dim.)

Quantum wire (1 dim.)

Quantum dot (0 dim.)

(42)

師大

Nanometer

1nm = 10⁻⁹m

Nano-scale ≈ 1–100 nm A semiconductor QD ≈ 10 nm QD : hair ≈ 1 : 10000

Why consider quantum effects?

Small devices imply significant quantum effect

(43)

師大

Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

150 Å

300 Å 150 Å

300 Å

(44)

師大

Molecular beam epitaxy

(45)

師大

Quantum dot

Cross-sections of hetero-structure InAs/GaAs QDs by Transmission Electron Microscope [Schoenfeld, 00]

150 Å

300 Å 150 Å

300 Å

(46)

師大

Energy levels (eigenvalues)

Confined Discrete Energy Levels Confined Discrete

Energy Levels

(47)

師大

Numerical experiments for linear eigenproblems

TheSchr¨odinger equationfor Semiconductor:

−∇ · (α∇u) + Vu = λu,

where α =

( α⁻≡ _2m^~²

1 inside, α⁺≡ _2m^~²

2 outside, V =

V⁻= V1 inside, V⁺= V₂ outside

~: Plank constant m_`: parabolic effective mass V_`: confinement potential λ: total energy

Interface condition:

α⁻∂u

∂n ∂D−

= α⁺∂u

∂n ∂D+

Dirichlet boundary conditions

(48)

師大

Symmetric eigenvalue problem

Ax = λx where A is a symmetric matrix.

Reference:

Tsung-Min Hwang, Wen-Wei Lin, Wei-Cheng Wang and Weichung Wang, Numerical simulation of three dimensional pyramid quantum dot, Journal of Computational Physics, Vol 196, pp. 208-232, 2004.

Figure: PRB, 54, 8743, (1996)

(49)

師大

Three dimensional pyramid quantum dot

H

W W

Pyramid Dot Cuboid Matrix i

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 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

(50)

師大

Finite volume discretized scheme Symmetric eigenvalue problems:

Ax = λx

Second order convergent rate Correction vector with SSOR preconditioner:

t = −M_A⁻¹r + εM_A⁻¹p

with

ε = u^TM_A⁻¹r u^TM_A⁻¹p,

whereM_A= (D + ωL)D⁻¹(D + ωU).

H

W W

Pyramid Dot Cuboid Matrix i

Figure: Structure schema of a pyramid quantum dot.

(51)

師大

The width of the QD (matrix) base is 12.4 nm (24.8 nm); the height of the QD (matrix) is 6.2 nm (18.6 nm).

InAs QD: α1 = 0.024me and V1 = 0.0 GaAs matrix: α2= 0.067me and V2= 0.70.

Stopping criteria: residual < 10⁻¹⁰ Convergent rate

(L,M,N) Mtx. dim. λ₁ Rate λ₂ Rate

(16, 16, 12) 2,475 0.4226 - 0.6527 - (32, 32, 24) 22,103 0.4001 - 0.6423 - (64, 64, 48) 186,543 0.3934 1.744 0.6391 1.708 (128,128, 96) 1,532,255 0.3916 1.905 0.6383 1.866 (256,256,192) 12,419,775 0.3911 1.954 0.6380 1.912 Table: Convergent rate =log₂ (λ^(4h)− λ^(2h))/(λ^(2h)− λ^(h))

(52)

師大

dim(A) = 1, 532, 255on PC with P4 1.8 GHz CPUand1 GB of main memory.

Value Ite. no. CPU time (sec.)

λ1 0.3916 72 278.0

λ₂ 0.6383 72 284.3

λ3 0.6383 125 521.4

dim(A) = 32, 401, 863on HP workstation with a 1.0 GHz Intel Itanium II CPUand 12 GBsof main memory.

Value Ite. no. CPU time (sec.)

λ₁ 0.3910 138 5,852

λ2 0.6380 133 5,354

λ3 0.6380 220 8,511

(53)

師大

Unsymmetric eigenvalue problem

Ax = λx where A is a unsymmetric matrix.

Reference:

Tsung-Min Hwang, Wei-Hua Wang and Weichung Wang, Efficient numerical schemes for electronic states in coupled quantum dots, accepted for publication in Journal of Nanoscience and

Nanotechnology.

JAP, 90-12, (2001)-

Figure: JAP, 90-12, (2001)

(54)

師大

Vertically aligned quantum dot array

(a) Nonuniform mes h (b) Uniform mesh

(55)

師大

Reduce three-dimensional systems to two-dimensional systems by Fourier transformation.

Finite volume discretized scheme Unsymmetric eigenvalue problems:

Ax = λx

Second order convergent rate

Method II with SSOR preconditioner:

M_A= (D + ωL)D⁻¹(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

(56)

師大

Ax = λx

M_A= (D + ωL)D⁻¹(D + ωU). _Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

(57)

師大

Ax = λx

M_A= (D + ωL)D⁻¹(D + ωU).

Figure: Structure schema of a cylindrical vertically aligned quantum dot array.

(58)

師大

gap = 6nm and dim(A) = 12, 288, 000on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ1 0.1587 83 20,287 0.1587 84 23,840

λ₂ 0.35553 109 26,030 0.3526 91 31,612

λ3 0.3558 53 13,831 0.3555 61 25,706

gap = 3nm and dim(A) = 11, 059, 200on HP workstation with a 1.0 GHz Intel Itanium II CPU and12 GBsof main memory.

r2 = 3.198 r2 = 3.223

Value Ite. no. Time (sec.) Value Ite. no. Time (sec.)

λ₁ 0.1586 85 20,072 0.1586 83 18,455

λ₂ 0.3540 257 56,811 0.3515 130 29,506

λ3 0.3564 53 12,257 0.3558 54 13,082

(59)

師大

Generalized eigenvalue problem

Ax = λBx

where A is symmetric positive definiteand B is a positive diagonalmatrix.

Reference:

Tsung-Min Hwang, Wei-Cheng Wang and Weichung Wang, Numerical schemes for three dimensional irregular shape quantum dots over curvilinear coordinate systems, accepted for publication in Journal of Computational Physics.

(60)

師大

Three dimensional arbitrary shape quantum dots

Quantum dot Matrix

Z_A ZB

ZC

1 M

Rmtx

0 H

N1

N2

N

H HC

HRmtx

R adial coordinate

Axial coordinate

(61)

師大

Curvilinear coordinate system

Jump condition capturing scheme (nine points)

Generalized eigenvalue problems:

Ax = λBx , A issymmetric positive definite and B is apositive diagonal matrix.

Method I with SSOR preconditioner:

M_A= (D + ωL)D⁻¹(D + ωU).

Quantum dot Matrix

Figure: Structure schema of the quantum dot model.