- Email: [email protected]
A new solution of Schrödinger equation based on symplectic algorithm
Computers and Mathematics with Applications 69 (2015) 1303–1312
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A new solution of Schrödinger equation based on symplectic algorithm Zhixiang Huang a,∗ , Jie Xu a , Bingbing Sun a , Bo Wu a , Xianliang Wu a,b,∗∗ a
Key Laboratory of Intelligent Computing and Signal Processing, Anhui University, Hefei 230039, China
b
Department of Physics and Electronic Engineering Hefei Normal University, Hefei 230061, China
article
info
Article history: Received 3 December 2013 Received in revised form 5 January 2015 Accepted 22 February 2015 Available online 9 April 2015 Keywords: Symplectic integrator Schrödinger equation Quantum mechanics High-order collocated difference
abstract An explicit three dimensional symplectic finite-difference time-domain method is introduced to solve the Schrödinger equation. The method is obtained by using a symplectic scheme with a propagator of exponential differential operators for the time direction approximation and fourth-order collocated finite differences for the space discretization. Firstly, a high-order symplectic framework for discretizing the Schrödinger equation is presented. Secondly, the comparisons on numerical stability, dispersion are also provided between the new symplectic scheme and other commonly used methods. Finally, numerical examples are given to evaluate the performance of the proposed method and the advantages on the accuracy and stability are also further demonstrated. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Developing an accurate and efficient numerical method to solve the Schrödinger equation is the main target in the regime of modeling and optimization of modern nanodevices [1]. Two typically available numerical methods are the Diffusion Monte Carlo (DMC) method [2–6] and Numerov technique [7]. In the DMC algorithm the resulting ‘‘dynamical’’ equations are transformed into an integral Green’s function form and then the resulting integral equation is computed using stochastic sampling. The method is highly inefficient unless importance sampling [8,9] is used. As for the matrix method, in cylindrical and spherical geometry it is possible to reduce the dimensionality by considering the symmetry of the problem, and rewriting it using a suitable coordinate system. The computational cost increases rapidly with matrix techniques when it is not possible to reduce the problem to a one-dimensional case. The finite difference time domain (FDTD) method has been another useful tool for solving the Schrödinger equation [10,11]. However, larger dispersion errors in a long time simulation is the bottleneck that has been unable to overcome yet. The eigenvalues and eigenfunctions of the device can be considered as the bases for analyzing the characteristics of current nanodevices. For example the ballistic transport of electrons depends mainly on the transverse eigenstates of the semiconducting channel, which is closely related to non-equilibrium Green’s function. Other interesting phenomena in quantum transport such as resonant tunneling effect [12], Fano resonance [13] are dominated by the excitement of the characteristic mode and coupling between different characteristic modes. So an accurate method for calculating the eigenvalues and eigenfrequencies of nanostructures plays an important role in the process of analysis and design of a device.
∗ ∗∗
Corresponding author. Corresponding author at: Key Laboratory of Intelligent Computing and Signal Processing, Anhui University, Hefei 230039, China. E-mail addresses: [email protected] (Z. Huang), [email protected] (X. Wu).
http://dx.doi.org/10.1016/j.camwa.2015.02.025 0898-1221/© 2015 Elsevier Ltd. All rights reserved.
1304
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
In addition, the extraction of the eigenvalues can meet the industrial demand of quantum CAD modeling which has profound effects on microscopic techniques. Many physical phenomena can be described by Hamilton differential equations which have time evolution and symplectic transformation [14–16]. Meanwhile time discretization methods were utilized in symplectic methods to preserve the global symplectic structure of the phase space for a Hamilton system. It is proven from many experiments that substantial advantages exit for a Hamilton system, especially in long-term simulations. A symplectic FDTD-quantum (SFDTD-Q) or namely SFDTD(4,4) scheme which is explicit third-order accurate in time was successfully introduced to solve the Schrödinger equation [17]. The currently available numerical methods fall into two categories. Methods in the first category rely on splitting the wave function into two parts: the real part and the imaginary part [18–20]. In the second category are other procedures that require dividing the Hamilton system into the kinetic energy operator and potential operator [21]. As for stationary Schrödinger equation, the SFDTD-Q can also be imposed under the condition of generalized coordinates and generalized velocity [22,23]. Moreover, it can be found that nonlinear Schrödinger equation can also be solved by the symplectic scheme [24,25]. We provided a two dimensional three-order symplectic FDTD method for solving Schrödinger equation in [17]. In this paper, we present a detailed formulation of the numerical SFDTD-Q(4, 4) technique for the general threedimensional Schrödinger equation with four order symplectic propagator technique. It is predicted that the method is quite general and has, we believe, some special advantages for certain problems. Then the stability and dispersion relation were discussed compared to other approaches. Finally, calculations with the SFDTD-Q(4, 4) have been carried out in two typical problems in which the advantages can be demonstrated. The paper is organized as follows. The formulation of the SFDTD-Q(4, 4) scheme is given in Section 2, followed by the comparisons on stability and dispersion error to other methods in Section 3. Two numerical experiments are presented in Section 4. Lastly, the summary is included in Section 5. 2. Symplectic framework for Schrödinger equation The time-dependent Schrödinger equation with a static potential V (r , t ) = V (r ) is given by: h¯ 2 2 ∂ψ(r , t ) =− ∇ ψ(r , t ) + V (r , t )ψ(r , t ) (1) ∂t 2m where ψ is the wave function that is a probability amplitude in quantum mechanism describing the quantum state of a ih¯
particle at position r and t, m is the mass of the particle being represented by the Schrödinger equation, usually an electron, 2
h¯ but which might be the effective mass of that particle in a particular semiconductor, − 2m ∇ 2 is the kinetic energy operator, 2
V (r ) is the potential of the structure being simulated, and − 2h¯∗m ∇ 2 + V is the Hamiltonian operator. To avoid using complex numbers, we will split ψ(r , t ) into its real and imaginary components:
ψ(r , t ) = ψR (r , t ) + j0 ψI (r , t ).
(2)
Inserting Eq. (2) into Eq. (1) and separating the real and imaginary parts result in the following coupled set of equations:
∂ψR (r , t ) h¯ 2 ∂ 2 ψI (r , t ) ∂ 2 ψI (r , t ) ∂ 2 ψI (r , t ) =− + + + V (r )ψI (r , t ) ∂t 2m ∂ 2x ∂ 2y ∂ 2z h¯ 2 ∂ 2 ψR (r , t ) ∂ 2 ψR (r , t ) ∂ 2 ψR (r , t ) ∂ψI (r , t ) h¯ = + + − V (r )ψR (r , t ). ∂t 2m ∂ 2x ∂ 2y ∂ 2z
h¯
(3)
(4)
Eqs. (3) and (4) can be rewritten in the form of an infinite-dimensional Hamiltonian system as
∂ ∂t
ψR ψI
0 0
A=
B=
ψR =L ψI
ψR = (A + B) ψI
(5a)
K 0
0 −K
K =−
h¯ 2m
0 0
∂2 ∂2 ∂2 + 2 + 2 2 ∂ x ∂ y ∂ z
+
V (r ) h¯
.
(5b)
Using the product of elementary symplectic mapping, the exact solution of Eq. (5) from t = 0 to t = 1t can be approximately constructed [26] exp(1t (A + B)) =
m l =1
exp(dl 1tB) exp(cl 1tA) + O(1t p+1 )
(6)
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1305
Table 1 Coefficients of the symplectic integrator propagators cl = cm+1−l (0 < l < m + 1); dl = dm−l (0 < l < m), dm = 0. Coefficients
2nd-order 2-stage
4th-order 5-stage
c1 d1 c2 d2 c3 d3
0.5 1.0 0.5 0.0 – –
0.17399689146541 0.62337932451322 −0.12038504121430 −0.12337932451322 0.89277629949778 d2
where cl and dl are the constant coefficients of the symplectic integrator and should be chosen in such a way to provide the highest possible value for p at a given integer number m. Here m is stage number needed in every integer time step, and p is the order of the approximation. Generally m > p. Coefficients cl and dl have been studied by many authors with various methods. In order to calculated the coefficients cl and dl , firstly, we expand both side of Eq. (7) into a matrix form. And then we compare the coefficients of the 1t p terms of expand form [15]. Therefore the pth-order conditions can be derived as follows:
cl1 dl2 cl3 dl4 . . . clp =
1≤l1 ≤l2
dl1 cl2 dl3 cl4 . . . dlp =
1≤l1
cl1 dl2 cl3 dl4 . . . dlp =
1≤l1 ≤l2
dl1 cl2 dl3 cl4 . . . clp =
1≤l1
1
p! 1
when p is odd integer .
(7a)
when p is even integer .
(7b)
p! 1
p! 1
p!
Because Eqs. (7a) and (7b) are not independent of each other, the time-reversible constrain are necessary as following cl = cm+1−l (0 < l < m + 1) . dl = dm−l (0 < l < m) , dm = 0
(7c)
Coefficients cl and dl can be obtained by solving Eqs. (7a)–(7c). The values of the coefficients for second- and fourth-order approximations are listed in Table 1. The exponential operator exp(dl 1tB) and exp(cl 1tA) can be analytically computed as follows: exp(dl 1tB) = I + dl 1tB
(8)
exp(cl 1tA) = I + cl 1tA.
(9)
Therefore, the pth-order decomposition method is applicable to the numerical integration of Schrödinger equation in the temporal domain as exp(1t (A + B)) =
m
(I + dl 1tB)(I + cl 1tA) + O(1t p+1 ).
(10)
l =1
A wave function of space and time evaluated at a discrete point in the Cartesian lattice and at a discrete stage in the time step can be notated as
l 1t ψ(x, y, z , t ) = ψ i1x, j1y, k1z , n + m
(11)
where 1x, 1y, 1z are the lattice space increments along the x, y, z direction. For the spatial direction, the explicit qth-order accurate centered difference expressions in conjugation with the collocated lattice are used to discretize the second-order spatial derivatives, shown in Eq. (12).
l
∂ 2 ψ n+ m ∂δ 2
q
=
2 q r =− 2
l
q +1
Wr ψ n+ m (h + r ) + O(∆δ
)
where δ = x, y, z , h = i, j, k, and Wr are the coefficients of spatial differences. The values of Wr are listed in Table 2.
(12)
1306
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312 Table 2 Coefficients of qth-order accurate centered differences. Order(q)
W−2
2 1 − 12
4
W−1
W0
W1
1
−2
1
4 3
− 52
4 3
W2
1 − 12
Based on the above description, the formulation of the real component of the wave function for the SFDTD-Q(4,q) scheme can be written as n+ ml
ψR
αx1 =
V (i, j, k)1t n+ l (i, j, k) + × ψI m (i, j, k) h¯ n+ ml n+ ml n+ ml − αx1 × ψI (i + 1, j, k) − 2ψI (i, j, k) + ψI (i − 1, j, k) n+ ml n+ ml n+ ml (i + 2, j, k) − 2ψI (i, j, k) + ψI (i − 2, j, k) − αx2 × ψI n+ ml n+ ml n+ ml − αy1 × ψI (i, j + 1, k) − 2ψI (i, j, k) + ψI (i, j − 1, k) n+ ml n+ ml n+ ml (i, j + 2, k) − 2ψI (i, j, k) + ψI (i, j − 2, k) − αy2 × ψI n+ ml n+ ml n+ ml − αz1 × ψI (i, j, k + 1) − 2ψI (i, j, k) + ψI (i, j, k − 1) n+ ml n+ ml n+ ml − αz2 × ψI (i, j, k + 2) − 2ψI (i, j, k) + ψI (i, j, k − 2) (l−1) n+ m
(i, j, k) = ψR
4 3
cl × Sx ,
αy1 =
4 3
cl × Sy ,
αz1 =
4 3
cl × Sz
−1
−1 −1 cl × Sx , αy2 = cl × Sy , αz2 = cl × Sz 12 12 12 h¯ 1t h¯ 1t h¯ 1t , Sy = , Sz = . Sx = 2m∆2x 2m∆2y 2m∆2z αx2 =
(13) (14) (15) (16)
For the cubic lattice case, ∆x = ∆y = ∆z = ∆δ , then Sx = Sy = Sz = Sδ . The uniform constant Sδ is defined similarly as the Courant–Friedrichs–Levy (CFL) number in electromagnetic field [27]. 3. Stability and dispersion error analysis We start the stability from the Neumann stability method. A superposition of plane-waves is the most general solution for a wave equation.
ψ(x, y, z , t ) = A0 exp(−j0 (i∆x kx + j∆y ky + k∆z kz ))
(17)
kx = k0 sin θ cos ϕ,
(18)
pm h¯
Here k0 =
ky = k0 sin θ sin ϕ,
kz = k0 cos θ .
is the wave number in quantum mechanics, pm is the momentum, and θ , ϕ are the spherical angles.
At the same time the operator ∇ 2 ψ is discretized by the qth-order collocated differences:
∂ 2ψ ∂ 2ψ ∂ 2ψ + + 2 2 ∂x ∂y ∂ z2 q 2 exp(−j0 rky ∆y ) exp(−j0 rkz ∆z ) exp(−j0 rkx ∆x ) + + = Wr ψ(i, j, k) ∆2x ∆2y ∆2z q
∇ 2ψ =
r =− 2
= (ηx + ηy + ηz )ψ = ηψ
(19)
where q
ηδ =
2 q
r =− 2
Wr
exp(−j0 rkδ ∆δ )
∆2δ
.
(20)
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1307
Table 3 Numerical stability for various algorithm. Algorithm
CFL-Q number
FDTD(2, 2)
1 d 0.8 d 1.5 d
FDTD(2, 4) SFDTD-Q(4, 4)
So the Schrödinger equation in a potential-free circumstance can be written in a matrix form as
0 ∂ ψR = h¯ ∂ t ψI
− η
h¯
η 2m ψR . ψ
(21)
I
0
2m Meanwhile the discretized evolution matrix S d can be obtained using the high-order symplectic integration scheme Sd =
S11 S21
S12 S22
=
m
1
0
h¯
ηdl 1t
1
l =1
1 0
−
h¯
ηcl 1t
2m 1
.
(22)
2m The eigenvalues λ of the above evolution matrix satisfy the following eigenequation
λ2 − (S11 + S22 )λ + (S11 S22 − S12 S21 ) = 0. (23) d d It is noticed that the trace is tr (S ) = S11 + S22 and the determinant of the evolution symplectic matrix is det(S ) = S11 S22 − S12 S21 = 1. Thus the eigenequation can be rewritten as
λ2 − tr (S d )λ + 1 = 0. The solutions of Eq. (24) are λ1,2 = From Eq. (22), we can get tr (S ) = 2 + d
m
(−1) gl l
l =1
gl =
h¯
tr (s)±j0
2
2m
(24)
√
4−tr (sd )2 , and the condition for stability is 2
|λ1,2 | = 1, that requires |tr (S d )| ≤ 2.
l η
∆2t 2
ci1 dj1 ci2 dj2 . . . cil djl +
1≤i1 ≤j1 ≤i2 ≤j2 <···
(25)
di1 cj1 di2 cj2 . . . dil cjl .
(26)
1≤i1 ≤j1 ≤i2 ≤j2 <···
The stability of the scheme is determined by the maximum of CFL in electromagnetics. Here we define CFL-Q = mh¯
1t . ∆2δ
The numerical results for CFL-Q are listed in Table 3 where d is the dimensional number. It is evident that the spatial high-order collocated differences decrease the CFL-Q number that can be improved by the high-order symplectic integrators. If we optimize the symplectic integrator, the stability of the SFDTD-Q(4, 4) scheme can become better. The dispersion relation of free electron is a paraboloid,
ω=
h¯ 2m
|k0 |2 .
(27)
h¯ We define a numerical velocity of the Schrödinger equation as v0 = ( 2m ). In symplectic scheme, the dispersion relation can be expressed as
cos(ω1t ) =
tr (S d )
.
2 The relative phase error is defined by
vp − v0 v
Err = 20 log10
(28)
(29)
0
where vp = ω2 and ω is calculated from Eq. (28). k0
The relative phase velocity error versus the propagating angle θ at φ = 90 is depicted in Fig. 1 while Fig. 2 shows the relative phase velocity error as a function of the propagating angle φ at θ = 0. It is observed from the figures that the fourth-order symplectic integrators have less dispersion than the standard FDTD(2, 2) but the same dispersion as second-order symplectic integrators FDTD(2, 4) which is quite different from the results studied in electromagnetics [28]. Fortunately, the high-order collocated differences allow coarse grids within a given error bound, which, in turn, results in shorter simulation time and less storage.
1308
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
Fig. 1. Dispersion curves versus the propagating angle θ at PPW = 7, CFL-Q = 0.25, φ = π/2.
Fig. 2. Dispersion curves versus the propagating angle φ at PPW = 7, CFL-Q = 0.25, θ = 0.
4. Numerical results 4.1. One-dimensional problem The following numerical experiments are based on a simulation of a particle which is initialized at a certain point interacting with a barrier of a slanting potential. The particle dynamically interact with the barrier results in parts of its kinetic energy for potential energy but the total energy remains the same. The reason for choosing to start with onedimensional (1D) problems is to demonstrate the advantage of the stability of the SFDTD-Q(4, 4) method. The validity and accuracy of the SFDTD-Q(4, 4) method will be evaluated by comparing the numerical results with the exact analytic solution on the scheme of error analysis in the following 3D problem. Here we will reconsider the parameter CFL-Q. The particle is given by a Gaussian wavepacket with a width of σ = 3 nm initialized at 10 nm point:
2 2π (k − k0 ) 1 k−k0 , ψR (k, 0) = e 2 σ cos λ 2 2π (k − k0 ) 1 k−k0 ψI (k, 0) = e 2 σ sin . λ Simulations are done with CFL-Q = 1.0, ∆x = 0.1 nm, m = 9.1 × 10−31 , h¯ = 1.05 × 10−34 . The potential is .2 Ve (x) = 0.2 − 040 x eV. We can see from Fig. 3 that under the condition of CFL-Q = 1.0, after a period of 120 time steps (Fig. 3(b)), the result from FDTD(2, 2) is unstable and begins to diverge while the SFDTD-Q(4, 4) algorithm still runs well and keeps stable and accurate even at a longer time. To verify the energy-conserving characteristic hold for the symplectic scheme, sum of the probability of finding the wavepacket along the x-axis is plotted in Fig. 4. Comparing with the fluctuation from traditional FDTD(2, 2), the SFDTD-Q(4, 4) scheme better preserves energy of the quantum system. 4.2. Three dimensional eigenvalue problem In this section, we first describe the theory for solving the Schrödinger eigenvalue problem, which utilizes a method that was developed earlier for determining the eigenvalues and eigenfunctions for the modes of waveguide from numerical
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
(a) 0 time step.
(c) 400 time step.
1309
(b) 120 time step.
(d) 1500 time step.
Fig. 3. CFL-Q = 1.0, the propagation of a particle interact with a slanting barrier using FDTD(2, 2) (a, b), and SFDTD-Q(4, 4) method (c, d) at different time step. The solid blue and dashed red curves are corresponding to the real and imaginary part of wave function, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Sum of probability at different time step.
solutions of the paraxial wave equation. Due to the similarity between paraxial wave equation and the Schrödinger equation, it is possible to apply the previously developed method to quantum mechanical problems with little necessary change. The eigenfrequencies are important quantities which other eigenvalues such as the eigenenergies and the eigenfunctions can be derived from. The method for calculating the eigenfrequencies and reconstructing the eigenfunctions is described in [29]. Here a similar procedure that involved in the application of the SFDTD-Q(4, 4) for the solution of the eigenfrequencies and eigenfunctions are given as follows:
1310
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
(a) Time domain results.
(b) Fourier transform results.
Fig. 5. Time domain data recorded at the source point and the corresponding Fourier transform results. Table 4 Eigenenergies (MeV) for 3D quantum infinite well. Quantum numbers
Analytic solution
FDTD(2, 2)
SFDTD-Q(4, 4)
(1, 1, 1) (1, 2, 1) (1, 3, 1) (2, 1, 1) (2, 2, 1) (2, 3, 1) (3, 1, 1) (3, 2, 1) (3, 3, 1) (4, 1, 1)
780 1080 1579 949 1248 1748 1230 1529 2029 1622
775 1072 1565 942 1242 1735 1220 1520 2010 1610
777 1073 1574 944 1246 1741 1227 1523 2024 1619
(a) A narrow test pulse is used as an excitation which is distributed in the problem space. As the SFDTD-Q(4, 4) simulation proceeds, time domain data is recorded at every time step. (b) When the SFDTD-Q(4, 4) simulation halted, Fourier transform is imposed on the recorded time domain data. The spectrum displays a set of sharp local maxima that are eigenfrequencies. (c) Eigenfunctions are reconstructed by a discrete Fourier transform (DFT) at each point in the problem space at each eigenfrequency in the process of a second SFDTD-Q(4, 4) simulation. Here we construct a problem domain which is 40 × 30 × 12 representing a 3D quantum well. The cell and the time step is chosen as 2.5 nm and 0.03 fs respectively. We directly set a test pulse in the middle of the problem domain. As the simulation proceeds, the outgoing and reflected wave function is mixed. The real part of the wave function is recorded every time step. Fig. 5 shows the recorded time data and its corresponding Fourier transform results. Peaks reveal the eigenenergies, where eigenfrequencies have been converted to energies by multiplication by Planck’s constant. It is different from the procedure done in [29] which simulated two times due to the source point we have chosen. The well symmetric characteristic ensures that the eigenfrequencies can be obtained as many as impossible. Meanwhile the analytic solution of the eigenenergies of a 3D infinite well can be obtained from the following equation: E (nx , ny , nz ) =
h¯ 2 2m
π nx Lx
2
+
π ny Ly
2
+
π nz Lz
2
.
(30)
Energy eigenvalues calculated from SFDTD-Q(4, 4) method are listed in Table 4 with values from FDTD(2, 2) and analytic solutions in [29] where available. The mode numbers are shown in the first column and the analytic solutions are in the second column. The simulated results of FDTD(2, 2) and SFDTD-Q(4, 4) are given in the third and fourth column respectively. From the table which contains the first 10 eigenvalues it is clearly seen that under the same condition the SFDTD-Q(4, 4) method can get a better agreement with the analytic values. For a 3D problem, there are so many modes. The above 10 eigenvalues are the lowest order levels but the most important levels. Fig. 6 illustrates several certain eigenfunctions calculated from the reconstructed procedure and the results are quite close to their analytic counterparts. Finally, the mean error between numerical eigenvalues and analytic solutions versus cell size are plotted in Fig. 7. The mean error is calculated averaging the error at 10 eigenenergies for a given cell size. It is clear that the rate of convergence has increased significantly with the use of SFDTD-Q(4, 4) scheme.
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1311
Fig. 6. Eigenfunctions distribution for the corresponding mode.
Fig. 7. Mean error at 10 eigenenergies for different cell size.
4.3. Dodecahedron potential problem As a nontrivial example of the broad applicability of the our method, we apply it to a potential that is a constant negative value of V = −100 inside a surface defined by a regular dodecahedron with the following 20 vertices
1 1 1 ± ,± ,± , 0, ± 2 , ±1 , ± 2 , ±1, 0 , ± 2 , 0, ±1 γ γ γ γ γ γ √ with γ = (1 + 5)/2. 1
1
1
(31)
We can numerically obtain the ground and first excited states extracted from a run on a 1283 lattice with a uniform lattice spacing of ∆δ = 0.1, a time step of 1t = 0.001 and particle mass of m = 1, h¯ = 1. We find that the energies of these two levels are E0 = −99.78049 and E1 = −99.55279. 5. Conclusions We have proposed a new method SFDTD-Q(p,q) to solve and determine the eigenvalues of the Schrödinger equation. The method relies on selecting a group of proper p and q and the principle here we follow is that high p should be coupled with high q (p = q = 4). Satisfying numerical stability and dispersion have been achieved by the SFDTD-Q(4, 4) scheme. From the perspective of numerical algorithm, it is qualified of SFDTD-Q(4, 4) method to serve as an effective and accurate tool in quantum simulations. Although we are concerned here only with some typical structures, it is hoped that our new scheme will be particularly useful in simulations with long time and large domain.
1312
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
Acknowledgments The work was supported by the NSFC (Nos. 61101064, 51277001, 61471001), Distinguished Natural Science Foundation (No. 1108085J01), Universities Natural Science Foundation of Anhui Province (No. KJ2012A013), DFMEC (No. 20123401110009) and NCET (NCET-12-0596) of China. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
S. Datta, Quantum Transport: Atom to Transistor, Cambridge University Press, New York, 2005. R. Grimm, R. Storer, Monte-Carlo solution of Schrödinger’s equation, J. Comput. Phys. 7 (1) (1971) 134–156. W. Foulkes, L. Mitas, R. Needs, G. Rajagopal, Quantum Monte Carlo simulations of solids, Rev. Modern Phys. 73 (1) (2001) 33–83. J. Anderson, A random-walk simulation of the Schrödinger equation: H+3, J. Chem. Phys. 63–67 (1975) 1499. L. Costa, F. Prudente, P. Acioli, J. Neto, J. Vianna, A study of confined quantum systems using the Woods–Saxon potential, J. Phys. B: At. Mol. Opt. Phys. 32 (1999) 2461–2470. I. Kosztin, B. Faber, K. Schulten, Introduction to the diffusion Monte Carlo method, Arxiv Preprint. F. Vesely, Computational Physics: An Introduction, Plenum Pub. Corp., 1994. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21 (1953) 1087–1092. C. Umrigar, Accelerated Metropolis method, Phys. Rev. Lett. 71 (3) (1993) 408–411. D. Sullivan, D. Citrin, Determination of the eigenfunctions of arbitrary nanostructures using time domain simulation, J. Appl. Phys. 91 (2002) 3219–3226. A. Soriano, E. Navarro, J. Portı, V. Such, Analysis of the finite difference time domain technique to solve the Schrödinger equation for quantum devices, J. Appl. Phys. 95 (12) (2004) 8011–8018. D. Griffiths, Introduction to Quantum Mechanics, Benjamin-Cummings Pub. Co., 2005. Y. Joe, A. Satanin, C. Kim, Classical analogy of Fano resonances, Phys. Scr. 74 (2006) 259–266. J. Sanz-Serna, M. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. W.E. Sha, Z. Huang, M. Chen, X. Wu, Survey on symplectic finite-difference time-domain schemes for Maxwell’s equations, IEEE Trans. Antennas and Propagation 56 (2) (2008) 493–500. K. Feng, M. Qin, Symplectic Geometry Algorithms for Hamiltonian Systems, Zhejiang Publ. United Group, Zhejiang Science and Technology Publ. House, 2010. J. Shen, W. Sha, Z. Huang, M. Chen, X. Wu, High-order symplectic FDTD scheme for solving time-dependent Schrödinger equation, Comput. Phys. Comm. 184 (2012) 480–492. S. Gray, D. Manolopoulos, Symplectic integrators tailored to the time-dependent Schrödinger equation, J. Chem. Phys. 104 (1996) 7099–7102. X. Liu, P. Ding, J. Hong, L. Wang, Optimization of explicit symplectic schemes for time-dependent Schrödinger equations, Comput. Math. Appl. 50 (3–4) (2005) 637–644. S. Blanes, F. Casas, A. Murua, Symplectic splitting operator methods for the time-dependent Schrödinger equation, J. Chem. Phys. 124 (2006) 234105–234114. S. Chin, C. Chen, Gradient symplectic algorithms for solving the Schrödinger equation with time-dependent potentials, J. Chem. Phys. 117 (2002) 1409–1415. X. Liu, X. Liu, Z. Zhou, P. Ding, S. Pan, Numerical solution of one-dimensional time-independent Schrödinger equation by using symplectic schemes, Int. J. Quantum Chem. 79 (6) (2000) 343–349. T. Monovasilis, Z. Kalogiratou, T. Simos, Computation of the eigenvalues of the Schrödinger equation by symplectic and trigonometrically fitted symplectic partitioned Runge–Kutta methods, Phys. Lett. A 372 (5) (2008) 569–573. A. Islas, D. Karpeev, C. Schober, Geometric integrators for the nonlinear Schrödinger equation, J. Comput. Phys. 173 (1) (2001) 116–148. T. Wang, T. Nie, L. Zhang, Analysis of a symplectic difference scheme for a coupled nonlinear Schrödinger system, J. Comput. Appl. Math. 231 (2) (2009) 745–759. H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (5) (1990) 262–268. A. Taflove, S. Hagness, Computational Electrodynamics, Artech House, Boston, 1995. W. Sha, Z. Huang, X. Wu, M. Chen, Application of the symplectic finite-difference time-domain scheme to electromagnetic simulation, J. Comput. Phys. 225 (1) (2007) 33–50. D. Sullivan, D. Citrin, Determining quantum eigenfunctions in three-dimensional nanoscale structures, J. Appl. Phys. 97 (2005) 104305–104310.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
A new solution of Schrödinger equation based on symplectic algorithm Zhixiang Huang a,∗ , Jie Xu a , Bingbing Sun a , Bo Wu a , Xianliang Wu a,b,∗∗ a
Key Laboratory of Intelligent Computing and Signal Processing, Anhui University, Hefei 230039, China
b
Department of Physics and Electronic Engineering Hefei Normal University, Hefei 230061, China
article
info
Article history: Received 3 December 2013 Received in revised form 5 January 2015 Accepted 22 February 2015 Available online 9 April 2015 Keywords: Symplectic integrator Schrödinger equation Quantum mechanics High-order collocated difference
abstract An explicit three dimensional symplectic finite-difference time-domain method is introduced to solve the Schrödinger equation. The method is obtained by using a symplectic scheme with a propagator of exponential differential operators for the time direction approximation and fourth-order collocated finite differences for the space discretization. Firstly, a high-order symplectic framework for discretizing the Schrödinger equation is presented. Secondly, the comparisons on numerical stability, dispersion are also provided between the new symplectic scheme and other commonly used methods. Finally, numerical examples are given to evaluate the performance of the proposed method and the advantages on the accuracy and stability are also further demonstrated. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction Developing an accurate and efficient numerical method to solve the Schrödinger equation is the main target in the regime of modeling and optimization of modern nanodevices [1]. Two typically available numerical methods are the Diffusion Monte Carlo (DMC) method [2–6] and Numerov technique [7]. In the DMC algorithm the resulting ‘‘dynamical’’ equations are transformed into an integral Green’s function form and then the resulting integral equation is computed using stochastic sampling. The method is highly inefficient unless importance sampling [8,9] is used. As for the matrix method, in cylindrical and spherical geometry it is possible to reduce the dimensionality by considering the symmetry of the problem, and rewriting it using a suitable coordinate system. The computational cost increases rapidly with matrix techniques when it is not possible to reduce the problem to a one-dimensional case. The finite difference time domain (FDTD) method has been another useful tool for solving the Schrödinger equation [10,11]. However, larger dispersion errors in a long time simulation is the bottleneck that has been unable to overcome yet. The eigenvalues and eigenfunctions of the device can be considered as the bases for analyzing the characteristics of current nanodevices. For example the ballistic transport of electrons depends mainly on the transverse eigenstates of the semiconducting channel, which is closely related to non-equilibrium Green’s function. Other interesting phenomena in quantum transport such as resonant tunneling effect [12], Fano resonance [13] are dominated by the excitement of the characteristic mode and coupling between different characteristic modes. So an accurate method for calculating the eigenvalues and eigenfrequencies of nanostructures plays an important role in the process of analysis and design of a device.
∗ ∗∗
Corresponding author. Corresponding author at: Key Laboratory of Intelligent Computing and Signal Processing, Anhui University, Hefei 230039, China. E-mail addresses: [email protected] (Z. Huang), [email protected] (X. Wu).
http://dx.doi.org/10.1016/j.camwa.2015.02.025 0898-1221/© 2015 Elsevier Ltd. All rights reserved.
1304
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
In addition, the extraction of the eigenvalues can meet the industrial demand of quantum CAD modeling which has profound effects on microscopic techniques. Many physical phenomena can be described by Hamilton differential equations which have time evolution and symplectic transformation [14–16]. Meanwhile time discretization methods were utilized in symplectic methods to preserve the global symplectic structure of the phase space for a Hamilton system. It is proven from many experiments that substantial advantages exit for a Hamilton system, especially in long-term simulations. A symplectic FDTD-quantum (SFDTD-Q) or namely SFDTD(4,4) scheme which is explicit third-order accurate in time was successfully introduced to solve the Schrödinger equation [17]. The currently available numerical methods fall into two categories. Methods in the first category rely on splitting the wave function into two parts: the real part and the imaginary part [18–20]. In the second category are other procedures that require dividing the Hamilton system into the kinetic energy operator and potential operator [21]. As for stationary Schrödinger equation, the SFDTD-Q can also be imposed under the condition of generalized coordinates and generalized velocity [22,23]. Moreover, it can be found that nonlinear Schrödinger equation can also be solved by the symplectic scheme [24,25]. We provided a two dimensional three-order symplectic FDTD method for solving Schrödinger equation in [17]. In this paper, we present a detailed formulation of the numerical SFDTD-Q(4, 4) technique for the general threedimensional Schrödinger equation with four order symplectic propagator technique. It is predicted that the method is quite general and has, we believe, some special advantages for certain problems. Then the stability and dispersion relation were discussed compared to other approaches. Finally, calculations with the SFDTD-Q(4, 4) have been carried out in two typical problems in which the advantages can be demonstrated. The paper is organized as follows. The formulation of the SFDTD-Q(4, 4) scheme is given in Section 2, followed by the comparisons on stability and dispersion error to other methods in Section 3. Two numerical experiments are presented in Section 4. Lastly, the summary is included in Section 5. 2. Symplectic framework for Schrödinger equation The time-dependent Schrödinger equation with a static potential V (r , t ) = V (r ) is given by: h¯ 2 2 ∂ψ(r , t ) =− ∇ ψ(r , t ) + V (r , t )ψ(r , t ) (1) ∂t 2m where ψ is the wave function that is a probability amplitude in quantum mechanism describing the quantum state of a ih¯
particle at position r and t, m is the mass of the particle being represented by the Schrödinger equation, usually an electron, 2
h¯ but which might be the effective mass of that particle in a particular semiconductor, − 2m ∇ 2 is the kinetic energy operator, 2
V (r ) is the potential of the structure being simulated, and − 2h¯∗m ∇ 2 + V is the Hamiltonian operator. To avoid using complex numbers, we will split ψ(r , t ) into its real and imaginary components:
ψ(r , t ) = ψR (r , t ) + j0 ψI (r , t ).
(2)
Inserting Eq. (2) into Eq. (1) and separating the real and imaginary parts result in the following coupled set of equations:
∂ψR (r , t ) h¯ 2 ∂ 2 ψI (r , t ) ∂ 2 ψI (r , t ) ∂ 2 ψI (r , t ) =− + + + V (r )ψI (r , t ) ∂t 2m ∂ 2x ∂ 2y ∂ 2z h¯ 2 ∂ 2 ψR (r , t ) ∂ 2 ψR (r , t ) ∂ 2 ψR (r , t ) ∂ψI (r , t ) h¯ = + + − V (r )ψR (r , t ). ∂t 2m ∂ 2x ∂ 2y ∂ 2z
h¯
(3)
(4)
Eqs. (3) and (4) can be rewritten in the form of an infinite-dimensional Hamiltonian system as
∂ ∂t
ψR ψI
0 0
A=
B=
ψR =L ψI
ψR = (A + B) ψI
(5a)
K 0
0 −K
K =−
h¯ 2m
0 0
∂2 ∂2 ∂2 + 2 + 2 2 ∂ x ∂ y ∂ z
+
V (r ) h¯
.
(5b)
Using the product of elementary symplectic mapping, the exact solution of Eq. (5) from t = 0 to t = 1t can be approximately constructed [26] exp(1t (A + B)) =
m l =1
exp(dl 1tB) exp(cl 1tA) + O(1t p+1 )
(6)
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1305
Table 1 Coefficients of the symplectic integrator propagators cl = cm+1−l (0 < l < m + 1); dl = dm−l (0 < l < m), dm = 0. Coefficients
2nd-order 2-stage
4th-order 5-stage
c1 d1 c2 d2 c3 d3
0.5 1.0 0.5 0.0 – –
0.17399689146541 0.62337932451322 −0.12038504121430 −0.12337932451322 0.89277629949778 d2
where cl and dl are the constant coefficients of the symplectic integrator and should be chosen in such a way to provide the highest possible value for p at a given integer number m. Here m is stage number needed in every integer time step, and p is the order of the approximation. Generally m > p. Coefficients cl and dl have been studied by many authors with various methods. In order to calculated the coefficients cl and dl , firstly, we expand both side of Eq. (7) into a matrix form. And then we compare the coefficients of the 1t p terms of expand form [15]. Therefore the pth-order conditions can be derived as follows:
cl1 dl2 cl3 dl4 . . . clp =
1≤l1 ≤l2
dl1 cl2 dl3 cl4 . . . dlp =
1≤l1
cl1 dl2 cl3 dl4 . . . dlp =
1≤l1 ≤l2
dl1 cl2 dl3 cl4 . . . clp =
1≤l1
1
p! 1
when p is odd integer .
(7a)
when p is even integer .
(7b)
p! 1
p! 1
p!
Because Eqs. (7a) and (7b) are not independent of each other, the time-reversible constrain are necessary as following cl = cm+1−l (0 < l < m + 1) . dl = dm−l (0 < l < m) , dm = 0
(7c)
Coefficients cl and dl can be obtained by solving Eqs. (7a)–(7c). The values of the coefficients for second- and fourth-order approximations are listed in Table 1. The exponential operator exp(dl 1tB) and exp(cl 1tA) can be analytically computed as follows: exp(dl 1tB) = I + dl 1tB
(8)
exp(cl 1tA) = I + cl 1tA.
(9)
Therefore, the pth-order decomposition method is applicable to the numerical integration of Schrödinger equation in the temporal domain as exp(1t (A + B)) =
m
(I + dl 1tB)(I + cl 1tA) + O(1t p+1 ).
(10)
l =1
A wave function of space and time evaluated at a discrete point in the Cartesian lattice and at a discrete stage in the time step can be notated as
l 1t ψ(x, y, z , t ) = ψ i1x, j1y, k1z , n + m
(11)
where 1x, 1y, 1z are the lattice space increments along the x, y, z direction. For the spatial direction, the explicit qth-order accurate centered difference expressions in conjugation with the collocated lattice are used to discretize the second-order spatial derivatives, shown in Eq. (12).
l
∂ 2 ψ n+ m ∂δ 2
q
=
2 q r =− 2
l
q +1
Wr ψ n+ m (h + r ) + O(∆δ
)
where δ = x, y, z , h = i, j, k, and Wr are the coefficients of spatial differences. The values of Wr are listed in Table 2.
(12)
1306
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312 Table 2 Coefficients of qth-order accurate centered differences. Order(q)
W−2
2 1 − 12
4
W−1
W0
W1
1
−2
1
4 3
− 52
4 3
W2
1 − 12
Based on the above description, the formulation of the real component of the wave function for the SFDTD-Q(4,q) scheme can be written as n+ ml
ψR
αx1 =
V (i, j, k)1t n+ l (i, j, k) + × ψI m (i, j, k) h¯ n+ ml n+ ml n+ ml − αx1 × ψI (i + 1, j, k) − 2ψI (i, j, k) + ψI (i − 1, j, k) n+ ml n+ ml n+ ml (i + 2, j, k) − 2ψI (i, j, k) + ψI (i − 2, j, k) − αx2 × ψI n+ ml n+ ml n+ ml − αy1 × ψI (i, j + 1, k) − 2ψI (i, j, k) + ψI (i, j − 1, k) n+ ml n+ ml n+ ml (i, j + 2, k) − 2ψI (i, j, k) + ψI (i, j − 2, k) − αy2 × ψI n+ ml n+ ml n+ ml − αz1 × ψI (i, j, k + 1) − 2ψI (i, j, k) + ψI (i, j, k − 1) n+ ml n+ ml n+ ml − αz2 × ψI (i, j, k + 2) − 2ψI (i, j, k) + ψI (i, j, k − 2) (l−1) n+ m
(i, j, k) = ψR
4 3
cl × Sx ,
αy1 =
4 3
cl × Sy ,
αz1 =
4 3
cl × Sz
−1
−1 −1 cl × Sx , αy2 = cl × Sy , αz2 = cl × Sz 12 12 12 h¯ 1t h¯ 1t h¯ 1t , Sy = , Sz = . Sx = 2m∆2x 2m∆2y 2m∆2z αx2 =
(13) (14) (15) (16)
For the cubic lattice case, ∆x = ∆y = ∆z = ∆δ , then Sx = Sy = Sz = Sδ . The uniform constant Sδ is defined similarly as the Courant–Friedrichs–Levy (CFL) number in electromagnetic field [27]. 3. Stability and dispersion error analysis We start the stability from the Neumann stability method. A superposition of plane-waves is the most general solution for a wave equation.
ψ(x, y, z , t ) = A0 exp(−j0 (i∆x kx + j∆y ky + k∆z kz ))
(17)
kx = k0 sin θ cos ϕ,
(18)
pm h¯
Here k0 =
ky = k0 sin θ sin ϕ,
kz = k0 cos θ .
is the wave number in quantum mechanics, pm is the momentum, and θ , ϕ are the spherical angles.
At the same time the operator ∇ 2 ψ is discretized by the qth-order collocated differences:
∂ 2ψ ∂ 2ψ ∂ 2ψ + + 2 2 ∂x ∂y ∂ z2 q 2 exp(−j0 rky ∆y ) exp(−j0 rkz ∆z ) exp(−j0 rkx ∆x ) + + = Wr ψ(i, j, k) ∆2x ∆2y ∆2z q
∇ 2ψ =
r =− 2
= (ηx + ηy + ηz )ψ = ηψ
(19)
where q
ηδ =
2 q
r =− 2
Wr
exp(−j0 rkδ ∆δ )
∆2δ
.
(20)
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1307
Table 3 Numerical stability for various algorithm. Algorithm
CFL-Q number
FDTD(2, 2)
1 d 0.8 d 1.5 d
FDTD(2, 4) SFDTD-Q(4, 4)
So the Schrödinger equation in a potential-free circumstance can be written in a matrix form as
0 ∂ ψR = h¯ ∂ t ψI
− η
h¯
η 2m ψR . ψ
(21)
I
0
2m Meanwhile the discretized evolution matrix S d can be obtained using the high-order symplectic integration scheme Sd =
S11 S21
S12 S22
=
m
1
0
h¯
ηdl 1t
1
l =1
1 0
−
h¯
ηcl 1t
2m 1
.
(22)
2m The eigenvalues λ of the above evolution matrix satisfy the following eigenequation
λ2 − (S11 + S22 )λ + (S11 S22 − S12 S21 ) = 0. (23) d d It is noticed that the trace is tr (S ) = S11 + S22 and the determinant of the evolution symplectic matrix is det(S ) = S11 S22 − S12 S21 = 1. Thus the eigenequation can be rewritten as
λ2 − tr (S d )λ + 1 = 0. The solutions of Eq. (24) are λ1,2 = From Eq. (22), we can get tr (S ) = 2 + d
m
(−1) gl l
l =1
gl =
h¯
tr (s)±j0
2
2m
(24)
√
4−tr (sd )2 , and the condition for stability is 2
|λ1,2 | = 1, that requires |tr (S d )| ≤ 2.
l η
∆2t 2
ci1 dj1 ci2 dj2 . . . cil djl +
1≤i1 ≤j1 ≤i2 ≤j2 <···
(25)
di1 cj1 di2 cj2 . . . dil cjl .
(26)
1≤i1 ≤j1 ≤i2 ≤j2 <···
The stability of the scheme is determined by the maximum of CFL in electromagnetics. Here we define CFL-Q = mh¯
1t . ∆2δ
The numerical results for CFL-Q are listed in Table 3 where d is the dimensional number. It is evident that the spatial high-order collocated differences decrease the CFL-Q number that can be improved by the high-order symplectic integrators. If we optimize the symplectic integrator, the stability of the SFDTD-Q(4, 4) scheme can become better. The dispersion relation of free electron is a paraboloid,
ω=
h¯ 2m
|k0 |2 .
(27)
h¯ We define a numerical velocity of the Schrödinger equation as v0 = ( 2m ). In symplectic scheme, the dispersion relation can be expressed as
cos(ω1t ) =
tr (S d )
.
2 The relative phase error is defined by
vp − v0 v
Err = 20 log10
(28)
(29)
0
where vp = ω2 and ω is calculated from Eq. (28). k0
The relative phase velocity error versus the propagating angle θ at φ = 90 is depicted in Fig. 1 while Fig. 2 shows the relative phase velocity error as a function of the propagating angle φ at θ = 0. It is observed from the figures that the fourth-order symplectic integrators have less dispersion than the standard FDTD(2, 2) but the same dispersion as second-order symplectic integrators FDTD(2, 4) which is quite different from the results studied in electromagnetics [28]. Fortunately, the high-order collocated differences allow coarse grids within a given error bound, which, in turn, results in shorter simulation time and less storage.
1308
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
Fig. 1. Dispersion curves versus the propagating angle θ at PPW = 7, CFL-Q = 0.25, φ = π/2.
Fig. 2. Dispersion curves versus the propagating angle φ at PPW = 7, CFL-Q = 0.25, θ = 0.
4. Numerical results 4.1. One-dimensional problem The following numerical experiments are based on a simulation of a particle which is initialized at a certain point interacting with a barrier of a slanting potential. The particle dynamically interact with the barrier results in parts of its kinetic energy for potential energy but the total energy remains the same. The reason for choosing to start with onedimensional (1D) problems is to demonstrate the advantage of the stability of the SFDTD-Q(4, 4) method. The validity and accuracy of the SFDTD-Q(4, 4) method will be evaluated by comparing the numerical results with the exact analytic solution on the scheme of error analysis in the following 3D problem. Here we will reconsider the parameter CFL-Q. The particle is given by a Gaussian wavepacket with a width of σ = 3 nm initialized at 10 nm point:
2 2π (k − k0 ) 1 k−k0 , ψR (k, 0) = e 2 σ cos λ 2 2π (k − k0 ) 1 k−k0 ψI (k, 0) = e 2 σ sin . λ Simulations are done with CFL-Q = 1.0, ∆x = 0.1 nm, m = 9.1 × 10−31 , h¯ = 1.05 × 10−34 . The potential is .2 Ve (x) = 0.2 − 040 x eV. We can see from Fig. 3 that under the condition of CFL-Q = 1.0, after a period of 120 time steps (Fig. 3(b)), the result from FDTD(2, 2) is unstable and begins to diverge while the SFDTD-Q(4, 4) algorithm still runs well and keeps stable and accurate even at a longer time. To verify the energy-conserving characteristic hold for the symplectic scheme, sum of the probability of finding the wavepacket along the x-axis is plotted in Fig. 4. Comparing with the fluctuation from traditional FDTD(2, 2), the SFDTD-Q(4, 4) scheme better preserves energy of the quantum system. 4.2. Three dimensional eigenvalue problem In this section, we first describe the theory for solving the Schrödinger eigenvalue problem, which utilizes a method that was developed earlier for determining the eigenvalues and eigenfunctions for the modes of waveguide from numerical
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
(a) 0 time step.
(c) 400 time step.
1309
(b) 120 time step.
(d) 1500 time step.
Fig. 3. CFL-Q = 1.0, the propagation of a particle interact with a slanting barrier using FDTD(2, 2) (a, b), and SFDTD-Q(4, 4) method (c, d) at different time step. The solid blue and dashed red curves are corresponding to the real and imaginary part of wave function, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 4. Sum of probability at different time step.
solutions of the paraxial wave equation. Due to the similarity between paraxial wave equation and the Schrödinger equation, it is possible to apply the previously developed method to quantum mechanical problems with little necessary change. The eigenfrequencies are important quantities which other eigenvalues such as the eigenenergies and the eigenfunctions can be derived from. The method for calculating the eigenfrequencies and reconstructing the eigenfunctions is described in [29]. Here a similar procedure that involved in the application of the SFDTD-Q(4, 4) for the solution of the eigenfrequencies and eigenfunctions are given as follows:
1310
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
(a) Time domain results.
(b) Fourier transform results.
Fig. 5. Time domain data recorded at the source point and the corresponding Fourier transform results. Table 4 Eigenenergies (MeV) for 3D quantum infinite well. Quantum numbers
Analytic solution
FDTD(2, 2)
SFDTD-Q(4, 4)
(1, 1, 1) (1, 2, 1) (1, 3, 1) (2, 1, 1) (2, 2, 1) (2, 3, 1) (3, 1, 1) (3, 2, 1) (3, 3, 1) (4, 1, 1)
780 1080 1579 949 1248 1748 1230 1529 2029 1622
775 1072 1565 942 1242 1735 1220 1520 2010 1610
777 1073 1574 944 1246 1741 1227 1523 2024 1619
(a) A narrow test pulse is used as an excitation which is distributed in the problem space. As the SFDTD-Q(4, 4) simulation proceeds, time domain data is recorded at every time step. (b) When the SFDTD-Q(4, 4) simulation halted, Fourier transform is imposed on the recorded time domain data. The spectrum displays a set of sharp local maxima that are eigenfrequencies. (c) Eigenfunctions are reconstructed by a discrete Fourier transform (DFT) at each point in the problem space at each eigenfrequency in the process of a second SFDTD-Q(4, 4) simulation. Here we construct a problem domain which is 40 × 30 × 12 representing a 3D quantum well. The cell and the time step is chosen as 2.5 nm and 0.03 fs respectively. We directly set a test pulse in the middle of the problem domain. As the simulation proceeds, the outgoing and reflected wave function is mixed. The real part of the wave function is recorded every time step. Fig. 5 shows the recorded time data and its corresponding Fourier transform results. Peaks reveal the eigenenergies, where eigenfrequencies have been converted to energies by multiplication by Planck’s constant. It is different from the procedure done in [29] which simulated two times due to the source point we have chosen. The well symmetric characteristic ensures that the eigenfrequencies can be obtained as many as impossible. Meanwhile the analytic solution of the eigenenergies of a 3D infinite well can be obtained from the following equation: E (nx , ny , nz ) =
h¯ 2 2m
π nx Lx
2
+
π ny Ly
2
+
π nz Lz
2
.
(30)
Energy eigenvalues calculated from SFDTD-Q(4, 4) method are listed in Table 4 with values from FDTD(2, 2) and analytic solutions in [29] where available. The mode numbers are shown in the first column and the analytic solutions are in the second column. The simulated results of FDTD(2, 2) and SFDTD-Q(4, 4) are given in the third and fourth column respectively. From the table which contains the first 10 eigenvalues it is clearly seen that under the same condition the SFDTD-Q(4, 4) method can get a better agreement with the analytic values. For a 3D problem, there are so many modes. The above 10 eigenvalues are the lowest order levels but the most important levels. Fig. 6 illustrates several certain eigenfunctions calculated from the reconstructed procedure and the results are quite close to their analytic counterparts. Finally, the mean error between numerical eigenvalues and analytic solutions versus cell size are plotted in Fig. 7. The mean error is calculated averaging the error at 10 eigenenergies for a given cell size. It is clear that the rate of convergence has increased significantly with the use of SFDTD-Q(4, 4) scheme.
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
1311
Fig. 6. Eigenfunctions distribution for the corresponding mode.
Fig. 7. Mean error at 10 eigenenergies for different cell size.
4.3. Dodecahedron potential problem As a nontrivial example of the broad applicability of the our method, we apply it to a potential that is a constant negative value of V = −100 inside a surface defined by a regular dodecahedron with the following 20 vertices
1 1 1 ± ,± ,± , 0, ± 2 , ±1 , ± 2 , ±1, 0 , ± 2 , 0, ±1 γ γ γ γ γ γ √ with γ = (1 + 5)/2. 1
1
1
(31)
We can numerically obtain the ground and first excited states extracted from a run on a 1283 lattice with a uniform lattice spacing of ∆δ = 0.1, a time step of 1t = 0.001 and particle mass of m = 1, h¯ = 1. We find that the energies of these two levels are E0 = −99.78049 and E1 = −99.55279. 5. Conclusions We have proposed a new method SFDTD-Q(p,q) to solve and determine the eigenvalues of the Schrödinger equation. The method relies on selecting a group of proper p and q and the principle here we follow is that high p should be coupled with high q (p = q = 4). Satisfying numerical stability and dispersion have been achieved by the SFDTD-Q(4, 4) scheme. From the perspective of numerical algorithm, it is qualified of SFDTD-Q(4, 4) method to serve as an effective and accurate tool in quantum simulations. Although we are concerned here only with some typical structures, it is hoped that our new scheme will be particularly useful in simulations with long time and large domain.
1312
Z. Huang et al. / Computers and Mathematics with Applications 69 (2015) 1303–1312
Acknowledgments The work was supported by the NSFC (Nos. 61101064, 51277001, 61471001), Distinguished Natural Science Foundation (No. 1108085J01), Universities Natural Science Foundation of Anhui Province (No. KJ2012A013), DFMEC (No. 20123401110009) and NCET (NCET-12-0596) of China. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
S. Datta, Quantum Transport: Atom to Transistor, Cambridge University Press, New York, 2005. R. Grimm, R. Storer, Monte-Carlo solution of Schrödinger’s equation, J. Comput. Phys. 7 (1) (1971) 134–156. W. Foulkes, L. Mitas, R. Needs, G. Rajagopal, Quantum Monte Carlo simulations of solids, Rev. Modern Phys. 73 (1) (2001) 33–83. J. Anderson, A random-walk simulation of the Schrödinger equation: H+3, J. Chem. Phys. 63–67 (1975) 1499. L. Costa, F. Prudente, P. Acioli, J. Neto, J. Vianna, A study of confined quantum systems using the Woods–Saxon potential, J. Phys. B: At. Mol. Opt. Phys. 32 (1999) 2461–2470. I. Kosztin, B. Faber, K. Schulten, Introduction to the diffusion Monte Carlo method, Arxiv Preprint. F. Vesely, Computational Physics: An Introduction, Plenum Pub. Corp., 1994. N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21 (1953) 1087–1092. C. Umrigar, Accelerated Metropolis method, Phys. Rev. Lett. 71 (3) (1993) 408–411. D. Sullivan, D. Citrin, Determination of the eigenfunctions of arbitrary nanostructures using time domain simulation, J. Appl. Phys. 91 (2002) 3219–3226. A. Soriano, E. Navarro, J. Portı, V. Such, Analysis of the finite difference time domain technique to solve the Schrödinger equation for quantum devices, J. Appl. Phys. 95 (12) (2004) 8011–8018. D. Griffiths, Introduction to Quantum Mechanics, Benjamin-Cummings Pub. Co., 2005. Y. Joe, A. Satanin, C. Kim, Classical analogy of Fano resonances, Phys. Scr. 74 (2006) 259–266. J. Sanz-Serna, M. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, London, 1994. W.E. Sha, Z. Huang, M. Chen, X. Wu, Survey on symplectic finite-difference time-domain schemes for Maxwell’s equations, IEEE Trans. Antennas and Propagation 56 (2) (2008) 493–500. K. Feng, M. Qin, Symplectic Geometry Algorithms for Hamiltonian Systems, Zhejiang Publ. United Group, Zhejiang Science and Technology Publ. House, 2010. J. Shen, W. Sha, Z. Huang, M. Chen, X. Wu, High-order symplectic FDTD scheme for solving time-dependent Schrödinger equation, Comput. Phys. Comm. 184 (2012) 480–492. S. Gray, D. Manolopoulos, Symplectic integrators tailored to the time-dependent Schrödinger equation, J. Chem. Phys. 104 (1996) 7099–7102. X. Liu, P. Ding, J. Hong, L. Wang, Optimization of explicit symplectic schemes for time-dependent Schrödinger equations, Comput. Math. Appl. 50 (3–4) (2005) 637–644. S. Blanes, F. Casas, A. Murua, Symplectic splitting operator methods for the time-dependent Schrödinger equation, J. Chem. Phys. 124 (2006) 234105–234114. S. Chin, C. Chen, Gradient symplectic algorithms for solving the Schrödinger equation with time-dependent potentials, J. Chem. Phys. 117 (2002) 1409–1415. X. Liu, X. Liu, Z. Zhou, P. Ding, S. Pan, Numerical solution of one-dimensional time-independent Schrödinger equation by using symplectic schemes, Int. J. Quantum Chem. 79 (6) (2000) 343–349. T. Monovasilis, Z. Kalogiratou, T. Simos, Computation of the eigenvalues of the Schrödinger equation by symplectic and trigonometrically fitted symplectic partitioned Runge–Kutta methods, Phys. Lett. A 372 (5) (2008) 569–573. A. Islas, D. Karpeev, C. Schober, Geometric integrators for the nonlinear Schrödinger equation, J. Comput. Phys. 173 (1) (2001) 116–148. T. Wang, T. Nie, L. Zhang, Analysis of a symplectic difference scheme for a coupled nonlinear Schrödinger system, J. Comput. Appl. Math. 231 (2) (2009) 745–759. H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (5) (1990) 262–268. A. Taflove, S. Hagness, Computational Electrodynamics, Artech House, Boston, 1995. W. Sha, Z. Huang, X. Wu, M. Chen, Application of the symplectic finite-difference time-domain scheme to electromagnetic simulation, J. Comput. Phys. 225 (1) (2007) 33–50. D. Sullivan, D. Citrin, Determining quantum eigenfunctions in three-dimensional nanoscale structures, J. Appl. Phys. 97 (2005) 104305–104310.