Dynamics Quantum Basics
In this chapter, we will simulate the dynamics of a particle, such as an electron, which follows the law of quantum mechanics [1]. Basics of the quantum-dynamics (QD) method [2-5] are described, along with corresponding data structures in program, qd.c.
Schrdinger Wave Equation
WAVE FUNCTION The state of an electron at time t is specified by a complex-valued wave function, ( r ,t) = Re( r ,t) + iIm( r ,t) C (where i = 1), which is spread in the 3-dimensional space, r = (x, y,z) R3. Given the wave function, we can calculate various physical properties such as: 2 2 2 P( r ,t) = * ( r ,t)( r ,t) = ( r ,t) = Re( r ,t) + Im( r ,t) : The probability to find the electron at position r at time t. 2 r (t) = dx dy dz ( r ,t) r : The expected position of the electron at time t. Here, * ( r ,t) = Re( r ,t) iIm( r ,t) is the complex conjugate of the wave function. Normalization: The electron wave function must be normalized such that the electron must be found somewhere in the entire space with probability 1, i.e., 2 (1) dx dy dz ( r ,t) = 1. WAVE EQUATION The time evolution of the electron state, subjected to a time-independent, real-valued potential, V ( r ), is described by the following partial differential equation:
i
2 ( r ,t) = 2 + V ( r )( r ,t), t 2m
(2)
where = 1.0545710-27 gcm2/s is the Planck constant, m = 9.1093810-28 g is the electron mass, and 2 2 2 2 = 2 + 2 + 2 is the Laplacian operator. x y z
Dimensionless equation: In the following, we measure length (x, y, z) in unit of 2 /me 2 = 0.52917710-8 cm, time (t) in unit of 3 /me 4 = 2.4188910-17 s, and energy (V) in unit of me 4 / 2 = 4.3597410-11 gcm2/s2, where e = 4.8032110-10 esu is the electron charge in the CGS unit. Substituting 2 2 2 x= x , y = y ,z = z me 2 me 2 me 2 3 t (3) t = 4 me me 4 V = 2 V in Eq. (2), we obtain
1
2 2 2 4 me 4 me 2 me i 3 ( r , t ) = + 2 V ( r ) ( r , t ), 2m 2 t
or the dimensionless equation,
2 ( r , t ) = + V ( r )( r , t ). t 2
i
(4)
In the following, we will use the dimensionless variables discussed above but omit the symbol for brevity. TWO-DIMENSIONAL ELECTRON As a specific example, the program qd.c simulates the time evolution of an electron confined in the 2-dimensional plane (z = 0). Such electrons are common at the interface between two heterogeneous materials in semiconductor devices. The electron state is now specified by the 2-dimensional wave function, (x, y,t) , where 0 x Lx and 0 y Ly (Lx and Ly are the system sizes in the x and y directions, respectively), and its time evolution is governed by the 2-dimensional Schrdinger equation,
In Eq. (5), the Hamiltonian operator, H, is defined as
i
(x, y,t) = H(x, y,t). t
1 2 1 2 + V (x, y) 2 x 2 2 y 2 .
(5)
H =
(6)
= Tx + Ty + V
Boundary condition: We impose the periodic boundary condition on the wave function such that
(x + L x , y) = (x, y) . (x, y + Ly ) = (x, y)
(7)
DISCRETIZATION The wave function is discretized on a regular mesh of size x and y in the x and y directions, respectively. Here x = Lx/Nx y = Ly/Ny, where Nx and Ny are the numbers of mesh points in the x and and y directions, respectively. We denote the discretized wave function as j,k = (jx, ky), see the figure below.
Figure: 2-dimensional mesh, on which the wave function is discretized.
2
In the discretized form, the Hamiltonian operator in Eq. (6) acts as
(H ) j,k = (Tx ) j,k + (Ty ) j,k + (V ) j,k ,
where
(8)
1 j1,k 2 j,k + j +1,k (Tx ) j,k = 2 (x) 2 1 j,k1 2 j,k + j,k +1 = Ty j,k 2 (y) 2 (V ) j,k = V j,k j,k
()
(9)
and the potential-energy function is discretized as Vj,k = V(jx, ky). Note that the discretized H is a mapping from an NxNy array (j,k is its element at the j-th column and k-th row) to another NxNy array H. The (j, k)-th element, (H)j,k, of the output array, H, is a linear combination of the input array with different indices, as specified in Eqs. (8) and (9).
Numerical Integration of Schrdinger Equation
The time evolution of the wave function is formally written down as
(t + t) = exp(iHt )(t) ,
(10)
where we omit the indices for simplicity. Here the exponential function of an operator is defined as a series expansion, 1 n exp(iHt ) = (iHt ) . (11) n=0 n!
We can show that Eq. (10) is the solution of the Schrdinger equation, Eq. (5), as follows:
(t + t) (t) i = i lim t0 t t exp(iHt )(t) (t) = i lim t0 t 1 2 1 iHt + (iHt) + (t) (t) 2 = i lim t0 t
= i lim iHt(t) + O (t) 2 t
(
)
t0
= H(t)
We can also show that Eq. (10) satisfies the correct initial condition at t = 0:
(t + 0) = exp(iH 0)(t) = 1(t) = (t) . Since Eq. (10) satisfies both the first-order differential equation and its initial condition, it is the solution of Eq. (5). // In the split-operator method [2-5], the wave function is propagated for a small time interval, t, as
3
(t + t) = exp(iVt /2) exp(iTx t ) exp(iTy t ) exp(iVt /2)(t)
+O [t ]
(
3
)
()
.
(12)
Note that, for operators (or for matrices) A and B, exp(A+B) exp(A)exp(B) in general. However, if both operators are multiplied by a small factor , the following symmetric splitting is a very accurate approximation: (13) exp((A + B)) = exp(A /2) exp(B) exp(A /2) + O 3 .
This can be easily proven by expanding both sides in and comparing the corresponding terms of the same power of : 2 exp((A + B)) = 1+ (A + B) + (A 2 + AB + BA + B 2 ) + O 3 2
()
A 2 A 2 A 2 A 2 2 B 2 exp(A /2) exp(B) exp(A /2) = 1+ + + 1+ B + + 1+ + + 2 8 2 2 8 A 2 B 2 A 2 AB BA A 2 A A = 1+ + B + + 2 + + + + + + O 3 2 2 8 2 8 2 2 4
()
Equation (12) is obtained by repeatedly applying this symmetric splitting operation. //
In Eq. (12), the application of the potential propagator, exp(-iVt/2), is straightforward.
(exp(iVt /2) )
j,k
= j,k = j,k
it (it /2) 2 (V ) j,k + 2! V 2 j,k + 2
()
2 it (it /2) 2 V j,k j,k + V j,k j,k + . 2 2! = exp iV j,k t /2 j,k
()
(14)
(
)
Note that, for real number a,
exp(ia) = 1+ ia +
1 1 1 1 (a 2 ) + (ia 3 ) + (a 4 ) + (ia 5 ) + 2! 3! 4! 4! 1 1 1 1 = 1 a 2 + a 4 + + i a a 3 + a 5 + . 2! 4! 3! 5! = cos(a) + isin(a)
(15)
Using Eq. (15) in Eq. (14),
(exp(iVt /2) ) j,k = [cos(V j,kt /2) isin(V j,kt /2)][Re j,k + iIm j,k ]
= cos V j,k t /2 Re j,k + sin V j,k t /2 Im j,k
) ( ) [( ] +i[cos(V j,k t /2) Im j,k sin(V j,k t /2) Re j,k ]
.
(16)
Many algorithms have been proposed to apply the kinetic propagators such as exp(-iTxt). Among these algorithms, the space-splitting method (SSM) [4,5] is highly scalable on massively parallel computers. To understand the SSM, we first note that the operation of Tx on j,k is expressed as Tx j,k = b j1,k + 2a j,k + b j +1,k (17) 4
where
a = 1/2 x 2 () b = 1/2(x ) 2
(18)
For each index k, the operation of Tx on j,k thus amounts to the multiplication of a tridiagonal matrix,
2a b b 2a b b 2a b Tx = . b 2a b b 2a b b 2a
In the SSM, this tridiagonal matrix is expressed as a direct sum of 2 2 submatrices,
2a b Tx = a b 1 = 2 b 2a b b 2a b b 2a b b 2a b b 2a a a b b a ab b a ba ab ab 1 ba ba + + 2 ab b a ba b a b a a b a
(19)
(20)
where we have omitted the index k. The exponential of a 22 matrix on the right-hand side of Eq. (20) is calculated analytically as follows:
(half) (full) (half) exp(itTx ) = U x U x U x + O (t) 3 =
(
)
+ + + 2 2 2 1 2 + + + 1 1 2 2 2 2 + + + 2 2 1 1 2 2 (21) + + 2 2 2 2 + 1 1 + + + 2 1 1 2 2 2 + + + 2 1 2 2 2
5
where
+ 1 it it (a + b) + exp (a b) n = exp n 2 n it it = 1 exp (a + b) exp (a b) n 2 n n
(22)
The operation of exp(-iTyt) is executed in a similar manner.
Equation (22) can be derived by utilizing the exponential of an 22 matrix, a b , b a
which can be calculated in terms of its eigenvalues, . The eigenvalues are the nontrivial solutions of the eigenvalue equation, a bu u (23) = , b av v or equivalently a b u 0 (24) = , b av 0 where u v is the corresponding eigenvector. In fact, Eq. (23) has a trivial solution, u = v = 0, for general values of , for which the determinant, a b D= = ( a) 2 b 2 , (25) b a takes a non-zero value. In such a case, Eq. (24) can be inverted as
u a b 10 a b 0 0 1 = = = . a0 0 v b a 0 ( a) 2 b 2 b
The nontrivial solutions, i.e., the eigenvalues, on the other hand, make the determinant D = ( - a)2 - b2 = 0. Let us denote the two solutions of this equation, D = 0, as = a b. The corresponding eigenvectors obtained by examining Eq. (24): are
+ a b u+ b bu+ 0 u+ 1/ 2 = = = b + av + b b v + 0 v + 1/ 2 . b u b bu 0 u 1/ 2 a = = = b av b bv 0 v 1/ 2
(26)
In Eq. (26), we have normalized the eigenvectors such that their norm is 1. Or the eigenvalue equation, Eq. (23), for these specific values reads
6
a b a b
1/ 2 b1/ 2 = (a + b) a1/ 2 1/ 2 . 1/ 2 b 1/ 2 = (a b) a1/ 2 1/ 2
(27)
The two matrix equations in Eq. (27) can be combined as
a b1/ 2 1/ 2 1/ 2 1/ 2 a + b 0 = , a b b a1/ 2 1/ 2 1/ 2 1/ 2 0
(28)
or equivalently
a b 1/ = b a 1/ 1/ = 1/
0 1/ 2 1/ 2 a + b a b1/ 2 1/ 2 0 0 1/ 2 1/ 2 a + b a b1/ 2 1/ 2 0
1 2 1/ 2 2 1/ 2 2 1/ 2 . 2 1/ 2
(29)
= UEU 1
Using Eq. (23), we can show that (a + b) n a bn n = UEU 1 = UEU 1UEU 1 UEU 1 = UE nU 1 = U b a 0
(
)
0 1 U . (a b) n
(30)
By using Eq. (30), we can now calculate the exponential of the 22 matrix as follows:
n a b n a bn n (a + b) n 0 1 exp U U = = (a b) n b a n=0 n! b a n=0 n! 0 n n (a + b) n 0 1 = U U n=0 n! 0 (a b) n exp( (a + b)) 1 0 = U . U 0 exp( (a b)) 1/ 2 1/ 2 exp( (a + b)) 1/ 2 1/ 2 0 = 0 exp( (a b))1/ 2 1/ 2 1/ 2 1/ 2 exp( (a + b)) + exp( (a b)) exp( (a + b)) exp( (a b)) 2 2 = exp( (a + b)) exp( (a b)) exp( (a + b)) + exp( (a b)) 2 2
(31)
Equation (21) is a straightforward application of the formula, Eq. (31). //
7
Data Structures of qd.c
NX, NY: Number of mesh points in the x and y directions. psi[NX+2][NY+2][2]: psi[i][j][0|1] is the real|imaginary
part of the wave function on mesh point
(i, j) in the xy plane. The wave function to be simulated is in the range, 1 i NX and 1 j NY. To simplify the operation of the finite-difference operators considering the periodic boundary condition, the wave function values at the edge are duplicated as follows:
for (sy=1; sy<=NY; sy++) for (s=0; s<=1; s++) { psi[0][sy][s] = psi[NX][sy][s]; psi[NX+1][sy][s] = psi[1][sy][s]; } for (sx=1; sx<=NX; sx++) for (s=0; s<=1; s++) { psi[sx][0][s] = psi[sx][NY][s]; psi[sx][NY+1][s] = psi[sx][1][s]; } v[NX+2][NY+2]: v[i][j] is the potential energy at mesh point (i, j). u[NX+2][NY+2][2]: u[i][j][0|1] is the real|imaginary part of the
potential propagator at mesh point
(i, j). The potential propagator, exp(-iVt/2), is operated in qd.c as follows, see Eq. (16).
for (sx=1; sx<=NX; sx++) for (sy=1; sy<=NY; sy++) { wr=u[sx][sy][0]*psi[sx][sy][0]-u[sx][sy][1]*psi[sx][sy][1]; wi=u[sx][sy][0]*psi[sx][sy][1]+u[sx][sy][1]*psi[sx][sy][0]; psi[sx][sy][0]=wr; psi[sx][sy][1]=wi; }
The program qd.c simulates an electron incident on a potential barrier of height BH and width BW, see the figure below. In addition, an edge potential of height EH is applied at i = 1 or NX or j = 1 or NY in v[i][j]. The potential v[i][j] = 0 at all the other mesh points.
Figure: potential The energy function.
8
In classical mechanics, a particle coming from one side of the potential barrier with a higher kinetic energy than BH will pass through the barrier to the other side; otherwise, the particle will bounce back at the barrier. In quantum mechanics, a part of the electron wave function is transmitted through the barrier and the other part is reflected at the barrier.
ND: The number of spatial dimensions = 2. al[ND][2][2]: al[0|1][0|1][0|1] is the
x|y-direction, half (t/2)|full (t)-step, real|imaginary-part of the diagonal element of the kinetic propagator, see Eqs. (21) and (22). bux|y[2][NX+2|Y+2][2]: bux|y[0|1][i][0|1] is the x|y-direction, half (t/2)|full (t)-step, real|imaginary-part of the upper off-diagonal kinetic propagator on mesh i, see Eqs. (21) and (22). blx|y[2][NX+2|Y+2][2]: blx|y[0|1][i][0|1] is the x|y-direction, half (t/2)|full (t)-step, real|imaginary-part of the lower off-diagonal kinetic propagator on mesh i, see Eqs. (21) and (22). The 22 block-diagonal form of the kinetic propagator, Eq. (21), can be handled
(U x(half) )i, j = mod(i,2),0i1, j + 2+i, j + mod(i,2),1i +1, j 2 2 (U x(full))i, j = 1mod(i,2),1i1, j + 1+i, j + 1mod(i,2),0i +1, j
(32) (33)
where mod(i,2),0 = 1 (if mod(i, 2) = 0) and 0 (else), etc. The above kinetic propagator is used in qd.c to update the wave function as follows, where d (= 0 for x; 1 for y) is the direction and t (= 0 for t/2 half; 1 for tfull) is the time step:
/* WRK|PSI holds the new|old wave function */ for (sx=1; sx<=NX; sx++) for (sy=1; sy<=NY; sy++) { wr=al[d][t][0]*psi[sx][sy][0]-al[d][t][1]*psi[sx][sy][1]; wi=al[d][t][0]*psi[sx][sy][1]+al[d][t][1]*psi[sx][sy][0]; if (d==0) { wr+=(blx[t][sx][0]*psi[sx-1][sy][0]-blx[t][sx][1]*psi[sx-1][sy][1]); wi+=(blx[t][sx][0]*psi[sx-1][sy][1]+blx[t][sx][1]*psi[sx-1][sy][0]); wr+=(bux[t][sx][0]*psi[sx+1][sy][0]-bux[t][sx][1]*psi[sx+1][sy][1]); wi+=(bux[t][sx][0]*psi[sx+1][sy][1]+bux[t][sx][1]*psi[sx+1][sy][0]); } else if (d==1) { wr+=(bly[t][sy][0]*psi[sx][sy-1][0]-bly[t][sy][1]*psi[sx][sy-1][1]); wi+=(bly[t][sy][0]*psi[sx][sy-1][1]+bly[t][sy][1]*psi[sx][sy-1][0]); wr+=(buy[t][sy][0]*psi[sx][sy+1][0]-buy[t][sy][1]*psi[sx][sy+1][1]); wi+=(buy[t][sy][0]*psi[sx][sy+1][1]+buy[t][sy][1]*psi[sx][sy+1][0]); } wrk[sx][sy][0]=wr; wrk[sx][sy][1]=wi; } /* Copy the new wave function back to PSI */ for (sx=1; sx<=NX; sx++) for (sy=1; sy<=NY; sy++) for (s=0; s<=1; s++) psi[sx][sy][s]=wrk[sx][sy][s];
INITIAL WAVE FUNCTION We choose the initial wave function at time t = 0 to be
y (x x ) 2 0 (x, y,t = 0) = C exp exp(ik0 x ) sin . L 4 2 y
To understand the meaning of this wave function, consider
(34)
9
(x, y) = C exp(ik0 x ) ,
where C = L x Ly is the normalization constant such that
|(x, y) |2 = 0 x dx 0 y dyC 2 (cos 2 (k0 x) + sin 2 (k0 x)) = 1. 0 x dx 0 y dy
Then
L L L L
(35)
(36)
1 2 2 H(x, y) = 2 + 2 C exp(ik 0 x ) 2 x y C d2 exp(ik 0 x ) 2 dx 2 ik C d = 0 exp(ik0 x ) 2 dx =
,
(37)
2 k0 C k2 exp(ik 0 x ) = 0 (x, y) 2 2 and thus the expectation value of its energy is
=
E 0 = 0 x dx 0 y dy * (x, y)H(x, y) = 0 x dx 0 y dy * (x, y) =
2 k0 2 L L 2 k0 (x, y) . 2
L
L
(38)
or k0 = 2E 0 . We can show that the following traveling wave function is a solution to the time-dependent Schrdinger equation, if the potential energy function is 0:
(x, y) = C exp(ik0 ( x v 0 t )) ,
where v 0 = E 0 /k0 is the wave speed.
(39)
The last factor in Eq. (34) also satisfies the Schrdinger equation and follows the boundary = condition, (x,0) = (x,Ly ) 0 , which is required if there is a very high potential barrier at y = 0 and Ly (the electron then cannot exist). Finally, the first Gaussian factor in Eq. (34) acts to localize the wave function around x = x0 and spread .
10
Quantum Dynamics Basics IIOne Dimensional System
This lecture note explains the quantum dynamics (QD) simulation program, qd1.c, which simulates the time evolution of a wave function in one spatial dimension. Here, we only highlight differences from the two-dimensional QD simulation program, qd.c, explained in the lecture note on Quantum Dynamics Basics.
Schrdinger Wave Equation
WAVE EQUATION The program qd1.c simulates the time evolution of a wave function in one-dimensional space, 0 x Lx (Lx is the system size). The one-dimensional wave function, (x,t) = Re(x,t) + iIm(x,t) C (where i = 1), satisfies the normalization condition: (1) The time evolution of the wave function, subjected to a time-independent, real-valued potential, V (x) , is described (in the dimensionless form as explained in the lecture note) by the time-dependent Schrdinger equation: (2) i (x,t) = H(x,t) . t In Eq. (2), the Hamiltonian operator, H, is defined as
1 2 + V (x) . 2 x 2 = Tx + V
2 dx (x,t) = 1.
H =
(3)
Boundary condition: We impose the periodic boundary condition on the wave function such that DISCRETIZATION
(x + L x ) = (x).
(4)
The wave function is discretized on a regular mesh of size x. Here x = Lx/Nx, and Nx is the number of mesh points. We denote the discretized wave function as j = (jx), see the figure below.
Figure: One-dimensional mesh, on which the wave function is discretized.
In the discretized form, the Hamiltonian operator in Eq. (3) acts as
(H ) j = (Tx ) j + (V ) j ,
where
(5)
11
1 j1 2 j + j +1 (Tx ) j = 2 (x) 2 (V ) = V j j j
(6)
and the potential-energy function is discretized as Vj = V(jx).
Numerical Integration of Schrdinger Equation
The time evolution of the wave function is formally written down as
(t + t) = exp(iHt )(t) ,
(7)
where we omit the indices for simplicity. Here the exponential function of an operator is defined as a series expansion, 1 n exp(iHt ) = (iHt ) . (8) n=0 n! In the split-operator method, the wave function is propagated for a small time interval, t, as
(t + t) = exp(iVt /2) exp(iTx t ) exp(iVt /2)(t) + O [t ] .
In Eq. (9), the application of the potential propagator, exp(-iVt/2), is straightforward.
(
3
)
(9)
(exp(iVt /2) ) j
it (it /2) 2 2 = j V (V ) j + 2 2!
( ) j +
(10)
2 it (it /2) 2 = j V j j + Vj j +. 2 2! = exp iV j t /2 j
()
(
)
Note that, for real number a,
exp(ia) = 1+ ia +
1 1 1 1 (a 2 ) + (ia 3 ) + (a 4 ) + (ia 5 ) + 2! 3! 4! 4! 12 14 1 1 = 1 a + a + + i a a 3 + a 5 + . 2! 4! 3! 5! = cos(a) + isin(a)
(11)
Using Eq. (11) in Eq. (10),
(exp(iVt /2) ) j = [cos(V jt /2) isin(V jt /2)][Re j + iIm j ]
= cos V j t /2 Re j + sin V j t /2 Im j
) ( )] [( +i[cos(V j t /2) Im j sin(V j t /2) Re j ]
.
(12)
To apply the kinetic propagator, exp(-iTxt), we use the space-splitting method (SSM). To understand the SSM, we first note that the operation of Tx on j is expressed as Tx j = b j1 + 2a j + b j +1 (13) where
12
a = 1/2 x 2 () b = 1/2(x ) 2
(14)
The operation of Tx on j thus amounts to the multiplication of a tridiagonal matrix,
2a b b 2a b b 2a b Tx = . b 2...