typeset using JPSJ.sty <ver.0.8>

Time Dependent Schrödinger Equations with Vector Potentials - Numerical Calculations and Visualizations -

HiroshiNatori^{∗} and TomoMunehisa^{∗∗}

Faculty of Engineering, Yamanashi University Takeda Kofu, Japan 400

(Received October 14, 1996 )

We study algorithm to calculate the time dependent Schrödinger equation numerically and to visualize the results. We extend the product formula by the Trotter-Suzuki decomposition to cases of Hamiltonians with vector potentials. This method is unconditionally stable and it can be used for a time dependent Hamiltonian so that one can apply it to various problems without any difficulty. After examining numerical errors for time evolutions of one dimensional wave packets, we present the errors under a constant magnetic field. Also we discuss a cyclotron motion in quantum mechanics, where an advantage for the unconditional stability is demonstrated. In order to visualize the results we make animations of the time evolutions, where the phase of the wave function is represented by colors and its absolute value is shown by brightness. Here we suggest a new way to represent a complex number by colors.

KEYWORDS: time dependent Schrödinger equation, vector potential, Trotter-Suzuki decomposition, product for- mula, unconditionally stable, visualization, mesoscopic, Aharonov-Bohm effect

§1. Introduction

Recent developments of experiments in fields of mesoscopic scale physics^{1)} as well as those of high performance computers stimulate us to study efficient methods for the time dependent Schrödinger equation,

ih̄\frac{∂}{∂t}ψ(x, t) = Hψ(x, t).  (1.1)

Also the computer graphics attracts our attention to making it easy to understand physical contents of quantum mechanics problems in researches as well as in educations. Here we study a method of solving the time dependent Schrödinger equation numerically and visualizations of these results using the computer graphics.

Since the intermediate wave functions are required to make animations on time evolutions of the wave function, the global propagation schemes such as the Chebyshev one^{2)} are not suitable for our purpose, though this is well known to be a very accurate method. On many applications to various problems in quantum mechanics including cases of time-dependent Hamiltonians we would like to obtain qualitatively correct results even for rough approximations. For this reason we adopt the product formula^{3)} based on the Trotter-Suzuki decomposition^{4)}, which has a property of the unconditional stability.

For the time dependent Schrödinger equation, the formal solution is

ψ(x, t) = exp(\frac{−it}{h̄}H)ψ(x, 0).  (1.2)

The exponential of the operators H is an object to be approximated. In the Trotter-Suzuki decomposition, H

∗ E-mail: natori@hep.esb.yamanashi.ac.jp ∗∗ E-mail: munehisa@hep.esb.yamanashi.ac.jp

are divided into several operators.

H = H_{1} + H_{2} + ⋯ + H_{K},  (1.3)

where it is assumed that one can diagonalize H_{i} analytically. Then

exp(\frac{−iHt}{h̄}) = exp(\frac{−i(H_{1} + H_{2} + ⋯ + H_{K})t}{h̄})

= lim_{N→∞}[exp(\frac{−iH_{1}t}{Nh̄})exp(\frac{−iH_{2}t}{Nh̄})⋯ exp(\frac{−iH_{K}t}{Nh̄})]^{N}.  (1.4)

When H is hermitian, i.e. exp(\frac{−iHt}{h̄}) is unitary, the above approximated expression is unitary if every H_{i} is hermitian. This guarantees the unconditional stability.

Here we comment on the stability and a time step Δt = t/N that control calculations. If the stability is kept, theoretical estimations on errors by Δt are reliable. While if the stability is not kept, it is difficult to obtain even qualitative results. An example that the stability is difficult to be kept is a time evolution with confining potentials, where diverge at large |x|. In this case the schemes with the conditional stability require extremely small Δt in order to exclude unphysical effects at large distance, because the stability condition is that E_{effective}Δt/h̄ ≪ 1 and E_{effective} becomes large at large distance for the confining potentials. On the other hand in the schemes with unconditional stability we can obtain qualitative results even for not small Δt.

Also actual visualizations do not always need extremely accurate results of numerical calculations, rather sometimes qualitative results are enough. When one makes animations for more complex systems by limited resources of computers, this point becomes important so that the unconditional stability is indispensable.

In this paper we extend the product formula to cases of Hamiltonians with vector potentials in order to study phenomena such as the Aharonov-Bohm effect. Our method is effective for a vector potential with any gauge condition so that one can confirm results of calculations by changing the gauge condition and checking them.

Next we visualize time evolutions of wave functions by using techniques of the computer graphics. In visualizations the phase is represented by colors and the absolute value is given by brightness. Here we suggest a new way to represent a complex number by colors. Some results on the animations are reported.

In the next section we make detailed descriptions of the product formula with vector potentials. Numerical results and error estimations of one dimensional wave packets are given in §3. Also in this section we present results of time evolutions of wave functions under a constant magnetic field. In §4 we discuss representations of a complex number by colors. The final section is devoted to conclusions and discussions.

§2. Product Formula with Vector Potentials

The Hamiltonian H which we consider hereafter is

H = \frac{1}{2m}(−ih̄∇ − eA(x))^{2}+V (x)
= \frac{1}{2m}(−h̄^{2}∇^{2} + ih̄e∇ · A(x) + ih̄eA(x) · ∇ + e^{2}A(x)^{2})+V (x).  (2.1)

Here A(x), V (x), and e denote vector potentials, scaler potentials and electric charges, respectively. In order of differentials, we divide H into three parts, T, Q and V̄.

T = −\frac{1}{2m}h̄^{2}∇^{2},  (2.2)
Q = \frac{ih̄e}{2m}(∇ · A(x) + A(x) · ∇),  (2.3)
V̄ = V (x) + \frac{e^{2}}{2m}|A(x)|^{2}.  (2.4)

For a given initial wave function ψ(x, 0) , the formal solution is given by eq. (1.2).

We will solve numerically the eq. (1.2) by the product formula. First the wave function is approximated by a discrete spatial representation, where the coordinate is given by a discrete number i and the space differences are Δx, Δy and Δz. If the space in one direction is represented by uniform L points, a wave vector of dimension L^{D}, which corresponds to the wave function, is

ψ(x_{i}, y_{j}, z_{k}, t) = ( ψ(x_{1}, y_{1}, z_{1}, t) ; ψ(x_{2}, y_{1}, z_{1}, t) ; ⋮ ; ψ(x_{L}, y_{1}, z_{1}, t) ; ψ(x_{1}, y_{2}, z_{1}, t) ; ψ(x_{2}, y_{2}, z_{1}, t) ; ⋮ ; ⋮ ; ψ(x_{L}, y_{L}, z_{L}, t) ) .  (2.5)

Here D is a number of space dimension and usually D = 1, 2, 3. Then the Hamiltonian becomes a L^{D} × L^{D} Hermitian matrix. In further expressions, variables y_{j}, z_{k} will be fixed.

First kinetic energy terms of the Hamiltonian are

T = T_{x} + T_{y} + T_{z},  (2.6)

T_{x} = −\frac{h̄^{2}}{2m} \frac{∂^{2}}{∂x^{2}} , T_{y} = −\frac{h̄^{2}}{2m} \frac{∂^{2}}{∂y^{2}} , T_{z} = −\frac{h̄^{2}}{2m} \frac{∂^{2}}{∂z^{2}} .

(2.7)

The second order expression of T_{x}ψ(x, t) in the discrete space representation is

− \frac{h̄^{2}}{2m} \frac{∂^{2}ψ(x, t)}{∂x^{2}} |_{x=x_{i},y=y_{j},z=z_{k}}= −\frac{h̄^{2}}{2m} \frac{ψ(x_{i+1}, y_{j}, z_{k}, t) − 2ψ(x_{i}, y_{j}, z_{k}, t) + ψ(x_{i−1}, y_{j}, z_{k}, t)}{(Δx)^{2}} + O((Δx)^{2}). (2.8)

Then a matrix for T_{x}, which hereafter [ ] denotes is

[T_{x}] = \frac{h̄^{2}}{2mΔx^{2}} ( \begin{matrix}2 & −1 & 0 & ⋯ & 0 & −1 \\ −1 & 2 & −1 & 0 & ... & 0 \\ 0 & −1 & 2 & −1 & ... & 0 \\ ⋮ & ... & ... & ... & ... & ⋮ \\ 0 & ... & ... & −1 & 2 & −1 \\ −1 & 0 & ⋯ & 0 & −1 & 2 \end{matrix} ) .

(2.9)

Next the terms of the first order differentials are

Q = Q_{x} + Q_{y} + Q_{z},  (2.10)

Q_{x} = \frac{ih̄e}{2m}(\frac{∂A_{x}}{∂x} + A_{x} \frac{∂}{∂x}),  (2.11)

Q_{y} = \frac{ih̄e}{2m}(\frac{∂A_{y}}{∂y} + A_{y} \frac{∂}{∂y} ),  (2.12)

Q_{z} = \frac{ih̄e}{2m}(\frac{∂A_{z}}{∂z} + A_{z} \frac{∂}{∂z} ).  (2.13)

The discrete representations of these differential operations are

\frac{∂A_{x}ψ(x, t)}{∂x}|_{x=x_{i}} = \frac{A_{x}(x_{i+1}, y_{j}, z_{j}, t)ψ(x_{i+1}, y_{j}, z_{k}, t) − A_{x}(x_{i−1}, y_{j}, z_{k}, t)ψ(x_{i−1}, y_{j}, z_{k}, t)}{2Δx} + O(Δx^{2}), (2.14)

A_{x} \frac{∂ψ(x, t)}{∂x}|_{x=x_{i}} = \frac{A_{x}(x_{i}, y_{j}, z_{k}, t)ψ(x_{i+1}, y_{j}, z_{k}, t) − A_{x}(x_{i}, y_{j}, z_{k}, t)ψ(x_{i−1}, y_{j}, z_{k}, t)}{2Δx} + O(Δx^{2}). (2.15)

Then Q_{x} is represented by

[Q_{x}] = \frac{ieh̄}{4mΔx}×

( \begin{matrix}0 & A_{x}^{(1)} + A_{x}^{(2)} & 0 & ⋯ & 0 & −A_{x}^{(L)} − A_{x}^{(1)} \\ −A_{x}^{(1)} − A_{x}^{(2)} & 0 & A_{x}^{(2)} + A_{x}^{(3)} & 0 & ... & 0 \\ 0 & −A_{x}^{(2)} − A_{x}^{(3)} & 0 & ... & ... & 0 \\ ⋮ & ... & ... & ... & A_{x}^{(L−2)} + A_{x}^{(L−1)} & ⋮ \\ 0 & ... & ... & −A_{x}^{(L−2)} − A_{x}^{(L−1)} & 0 & A_{x}^{(L−1)} + A_{x}^{(L)} \\ A_{x}^{(L)} + A_{x}^{(1)} & 0 & ⋯ & 0 & −A_{x}^{(L−1)} − A_{x}^{(L)} & 0 \end{matrix} ) , (2.16)

where indices of y, z are omitted. A_{x}^{(i)} mean A_{x}(x_{i}, y_{j}, z_{k}). For terms of the scaler potentials and squared vector potentials the matrix for V̄ is diagonal.

[V̄] = ( \begin{matrix}V^{(1)} + \frac{e^{2}|A^{(1)}|^{2}}{2m} & 0 & ⋯ & 0 \\ 0 & V^{(2)} + \frac{e^{2}|A^{(2)}|^{2}}{2m} & ... & ⋮ \\ ⋮ & ... & ... & 0 \\ 0 & ⋯ & 0 & V^{(L)} + \frac{e^{2}|A^{(L)}|^{2}}{2m} \end{matrix} ) .  (2.17)

By adding these matrices, the Hamiltonian matrices of x component are obtained. The matrices for y and z components are calculated likewise.

As a result, the matrix of the Hamiltonian H is given by

[H] = [H_{x}] + [H_{y}] + [H_{z}],  (2.18)

[H_{x}] = [T_{x}] + [Q_{x}] + [V]/3,
[H_{y}] = [T_{y}] + [Q_{y}] + [V]/3,
[H_{z}] = [T_{z}] + [Q_{z}] + [V]/3.  (2.19)

Where we divide [V̄] into three equal parts in order to treat symmetrically of the x−, y− and z−components. The explicit form of [H_{x}] is

[H_{x}] = ( \begin{matrix}2v_{1} & b + ic_{1} & 0 & ⋯ & 0 & b − ic_{L} \\ b − ic_{1} & 2v_{2} & b + ic_{2} & 0 & ⋯ & 0 \\ 0 & b − ic_{2} & 2v_{3} & b + ic_{3} & 0 & ... 0 \\ ⋮ & ... & ... & ... & ... & ⋮ \\ 0 & ⋯ & 0 & b − ic_{L−2} & 2v_{L−1} & b + ic_{L−1} \\ b + ic_{L1} & 0 & ⋯ & 0 & b − ic_{L−1} & 2v_{L} \end{matrix} ) ,  (2.20)

where

2v_{i} = \frac{h̄^{2}}{mΔx^{2}} + \frac{1}{3}(V^{(i)} + \frac{e^{2}|A^{(i)}|^{2}}{2m} ),
b = −\frac{h̄^{2}}{2mΔx^{2}} ,
c_{m} = \frac{eh̄}{4mΔx}(A_{x}^{(m)} + A_{x}^{(m+1)}).  (2.21)

Here m = 1, .., L and L + 1 is identified with 1 due to the periodicity of the space. But the exponential of [H_{x}] is impossible to be calculated except for small L cases. In the product formula, this exponentiation is approximated by products of exponentials of some matrices, each of which are easy to be calculated. Here we split the matrix [H_{x}] eq. (2.20) into two parts [H_{x1}] and [H_{x2}] so that both matrices are direct products of 2 × 2 matrices. They are

[H_{x1}] = ( \begin{matrix}v_{1} & b + ic_{1} & 0 & 0 & ⋯ & 0 & 0 \\ b − ic_{1} & v_{2} & 0 & 0 & ⋯ & ⋯ & 0 \\ 0 & 0 & v_{3} & b + ic_{3} & 0 & ... & 0 \\ 0 & 0 & b − ic_{3} & v_{4} & 0 & ... & 0 \\ ⋮ & ... & ... & ... & ... & ... & ⋮ \\ 0 & ⋯ & ⋯ & 0 & 0 & v_{L−1} & b + ic_{L−1} \\ & 0 & ⋯ & ⋯ & 0 & b − ic_{L−1} & v_{L} \end{matrix} ) ,  (2.22)

[H_{x2}] = ( \begin{matrix}v_{1} & 0 & 0 & ⋯ & ⋯ & 0 & b − ic_{L} \\ 0 & v_{2} & b + ic_{2} & 0 & ⋯ & ⋯ & 0 \\ 0 & b − ic_{2} & v_{3} & 0 & 0 & ⋯ & 0 \\ ⋮ & ... & ... & ... & ... & ⋮ \\ 0 & ⋯ & ⋯ & 0 & v_{L−2} & b + ic_{L−2} & 0 \\ 0 & ⋯ & ⋯ & 0 & b − ic_{L−2} & v_{L−1} & 0 \\ b + ic_{L} & 0 & ⋯ & ⋯ & ⋯ & 0 & v_{L} \end{matrix} ) .  (2.23)

For [H_{x1}] and [H_{x2}] basic form of 2 × 2 matrices [M^{(i)}] are the same, namely

[M]^{(i)} ≡ ( \begin{matrix}v_{i} & b + ic_{i} \\ b − ic_{i} & v_{i+1} \end{matrix} ) = ( \begin{matrix}a_{i} & 0 \\ 0 & a_{i} \end{matrix} ) + ( \begin{matrix}d_{i} & b + ic_{i} \\ b − ic_{i} & −d_{i} \end{matrix} ) ,  (2.24)

a_{i} = \frac{v_{i} + v_{i+1}}{2} ,  d_{i} = \frac{v_{i} − v_{i+1}}{2} .  (2.25)

Then exp(−i[H_{xi}]Δt/h̄) is described by a direct product of 2 × 2 matrices [U]^{(i)} defined by

[U]^{(i)} ≡ exp(−i[M]^{(i)}Δt/h̄)  (2.26)
= ( \begin{matrix}\frac{(λ−d)^{2} exp(i(λ−a)Δt)+(b^{2}+c^{2}) exp(−i(λ+a)Δt)}{2λ(λ−d)} & −\frac{(b+ic) exp(i(λ−a)Δt)−(b+ic) exp(−i(λ+a)Δt)}{2λ} \\ −\frac{(b−ic) exp(i(λ−a)Δt)−(b−ic) exp(−i(λ+a)Δt)}{2λ} & \frac{(b^{2}+c^{2}) exp(i(λ−a)Δt)+(λ−d)^{2} exp(−i(λ+a)Δt)}{2λ(λ−d)} \end{matrix} ) , (2.27)

λ = \sqrt{b^{2} + c^{2} + d^{2}},  (2.28)

where indices of a, c and d are omitted. exp(\frac{−iΔt}{h̄}[H_{x1}]) expressed by [U^{(i)}] is

exp(−i[H_{x1}]Δt/h̄)

= ( \begin{matrix}[U^{(1)}] & 0 & . & . & 0 \\ 0 & [U^{(3)}] & . & . & 0 \\ . & . & ... & . & 0 \\ . & . & . & ... & 0 \\ 0 & 0 & 0 & 0 & [U^{(L−1)}] \end{matrix} ) . (2.29)

Thus we obtain an explicit expression in the product formula.

exp(\frac{−iΔt}{h̄}[H_{x}]) = exp(\frac{−iΔt}{h̄}[H_{x1}])× exp(\frac{−iΔt}{h̄}[H_{x2}])  (2.30)
+O(Δt^{2}).  (2.31)

After including contributions of [H_{y}], [H_{z}], the time evolution of the wave function is given by

ψ(x, t + Δt) = exp(−i[H_{x1}]Δt/h̄) exp(−i[H_{x2}]Δt/h̄)
× exp(−i[H_{y1}]Δt/h̄) exp(−i[H_{y2}]Δt/h̄)
× exp(−i[H_{z1}]Δt/h̄) exp(−i[H_{z2}]Δt/h̄)
× ψ(x, t).  (2.32)

Finally we would like to comment on the symmetric scheme of the product formula, which is

exp(\frac{−iHΔt}{h̄}) = exp(\frac{−i(H_{1} + H_{2} + ⋯ + H_{K})Δt}{h̄})

= exp(\frac{−iH_{1}Δt}{2h̄})exp(\frac{−iH_{2}Δt}{2h̄})⋯ exp(\frac{−iH_{K}Δt}{h̄})⋯ exp(\frac{−iH_{2}Δt}{2h̄})exp(\frac{−iH_{1}Δt}{2h̄})
+O(Δt^{3}).  (2.33)

We do not use this symmetric scheme explicitly because even in the non-symmetric scheme, we obtain the same order of approximations. In the symmetric scheme terms [H_{i}, H_{j}] are canceled out, while in the non-symmetric scheme, they are not. However in actual applications of the product formula,

exp(\frac{−iHt}{h̄})≈ [exp(\frac{−iHΔt}{h̄})]^{N},  (2.34)

where Δt = t/N, N is large so that [H_{i}, H_{j}] at different t are canceled out. This is demonstrated by numerical calculations. We do not adopt the symmetric scheme, because the symmetric scheme needs more times for the matrix operations to obtain the same order of the approximation.

§3. Numerical results

3.1 One dimensional wave packet

First we study a free wave packet in one dimensional space to examine accuracy of the product formula. Since analytical solutions in discrete space are well known in this case, errors are only due to a finite time step. The Schrödinger equation in the discrete space without potential terms is

− \frac{h̄^{2}}{2m} \frac{ψ(x_{i+1}, t) − 2ψ(x_{i}, t) + ψ(x_{i−1}, t)}{(Δx)^{2}} = \frac{∂ψ(x_{i}, t)}{∂t}.

(3.1)

The solution of this equation under the periodic boundary condition is

ψ(x_{m}, t) = \frac{1}{\sqrt{L}} ∑_{n=0}^{L−1} C_{In}C_{Tn} exp(i2πmn/L), (3.2)

C_{In} = \frac{1}{\sqrt{L}} ∑_{n=0}^{L−1} exp(−i2πmn/L)ψ(x_{m}, 0), (3.3)
C_{Tn}= exp(−iE_{n}/h̄t),  (3.4)
E_{n} = \frac{sin^{2}(πnΔx/L)}{(Δx)^{2}}.  (3.5)

Two kinds of errors, the phase and norm error are often used for error estimations^{2)},^{5)}. They are defined as

ϵ_{phase} = |arg[< ψ_{num}|ψ_{exact} >]|,  (3.6)

ϵ_{norm} = 1 − | < ψ_{num}|ψ_{exact} > |.  (3.7)

The phase error is sensitive to derivations on the phase of the wave functions, while the norm error is sensitive to derivations on their amplitude. An initial state of the wave packet is

ψ(x, 0) = exp[−(\frac{x}{σ_{x}} )^{2}+ip_{x}x],  (3.8)

and numerical parameters are shown in Table I.

Here h̄ and mass m is set to be units. Numerical results of the phase and norm errors are shown in Fig. 1, as a function of the propagation time t. One finds that the phase errors are proportional to it, ϵ_{phase} ∝ t, in agreement with results by ref.5. The norm errors in Fig. 1 are smaller than the phase errors and show a weak dependence on t. Fig. 2 plots errors for various values of Δt, where the propagation time is fixed to be 100 so that a step number N increases as Δt becomes small. This figure shows that the errors are linear functions of (Δt)^{2}, which implies that the non-symmetric product formula we adopt is actually the approximation of the same order of the symmetric scheme, which is discussed in the end of a previous section. These results show that by the product formula we can evaluate enough accurate solution.

3.2 Wave functions under a constant magnetic field

Next we present numerical results of two-dimensional states in a constant magnetic field, where eigen states are called Landau states. A vector potential is fixed to be

A_{x} = −By, A_{y} = 0,  (3.9)

where B is a magnetic field. Because we can not obtain exact solutions in a discrete space, we compare our results with exact solutions in a continuous space. Then Landau states are represented by

ϕ_{n} = exp(ip_{x}x/h̄) exp[−\frac{|eB|}{2h̄}(y − \frac{p_{x}}{eB} )^{2}]H_{n}[\sqrt{\frac{|eB|}{h̄}}(y − \frac{p_{x}}{eB} )],  (3.10)

where H_{n} is the Hermite polynomial. These solutions are periodic in the x direction, while their values go to zero as |y| becomes large. Hence we use a periodic boundary condition in the x direction and a fixed boundary condition in the y direction. Also numerical values of parameters are shown in Table II.

We show the phase and norm errors in Fig. 3 as functions of the propagation time. Results are similar to those in the case of the free wave packet. The phase errors are linear on the propagation time, while the norm errors are smaller than the phase errors. The latter does not show a simple t−dependence. Next Fig.4 will give us the errors for various values of Δt. Here the propagation time is fixed to be 2π/E_{2}, where the energy E_{2} = 5/2. When Δt ≥ 0.3 × (Δx)^{2} , the errors show dependences on (Δt)^{2}, which implies that the product formula is same order of the symmetric scheme as pointed out in the previous subsection. As the time step becomes smaller than 0.1×(Δx)^{2} the errors do not decrease. A reason for this is that the errors due to Δx and Δy dominate those due to finite Δt. To confirm this consideration, we estimate the errors for L_{x} × L_{y} = 200 × 800, that is Δx = Δy = 0.02. The phase errors for Δx = Δy = 0.02 in Δt < 2 × 10^{−5} become the half of them for Δx = Δy = 0.04, while the norm errors become the fourth. This fact shows us that calculations with extremely small Δt are not useful, or if one wants to obtain accurate calculations, Δx and Δy should be small as well as Δt.

In the previous example, the absolute value |ϕ_{n}| of the wave function does not move because it is a eigen state. Here we study a wave packet in a constant magnetic field. A corresponding classical problem is a cyclotron motion, where any particle shows a circle motion and this period T_{c} is 2πm/eB. In the quantum problem the circle motion and the periodicity are realized. Any wave packet shows a circle motion and comes back to the initial state after 2T_{c}. The absolute value |ψ| comes back after T_{c}.

Here the initial state and vector potential are

ψ(x, 0) = exp[−(\frac{x − x_{0}}{σ_{x}} )^{2}−(\frac{y − y_{0}}{σ_{y}} )^{2}+i(p_{x}x + p_{y}y)],

(3.11)

A_{x} = \frac{−By}{2} , A_{y} = \frac{Bx}{2}.  (3.12)

Also numerical values of parameters are shown in Table III.

Fig. 5 shows expectation values of locations in the wave packet during a rotation. Points move clockwise and their trajectory shows an exact circle. The locations after T_{c} are denoted by symbols of a filled diamond and a filled square. One finds the circle motion even for rough approximations of Δt = 1.6 × 10^{−3}. This gives us a clear advantage of the product formula. In Fig. 6 we can see differences between the average locations after T_{c} and at the initial time. For Δt ≥ 1.6 × 10^{−4} the differences are larger than 1. These figures suggest us that the approximated wave function ψ_{A}(x, t) with finite Δt equals to the exact one ψ(x, t_{e}) with an effective time t_{e} = t − t_{D}, t_{D} ≥ 0, i.e. the error can be represented by a delay of time.

§4. Visualizing wave functions

In this section we would like to present the time evolutions of the wave functions in a visual manner. The basic idea to represent a complex wave function ψ =| ψ | exp(iθ) is to make | ψ | correspond to brightness and θ to a color. For concreteness, the red, green and blue color represent θ = 0, 2π/3 and 4π/3, respectively.

In general, colors of computer displays are given by normalized magnitudes [0, 1] of three colors, red, green and blue. A large magnitude gives the bright color while a small one does the dark one. So any color can be corresponded with a point of a cube [0, 1]^{3} (see Fig. 7), that is called a color cube. In order to clarify definitions on colors, let us call a point inside the cube a color. Also the maximally monochromatic color is called the maximum color, which is a point on sides that connect the red [1, 0, 0] and the yellow [1, 1, 0], the yellow and the green [0, 1, 0], the green and the cyan [0, 1, 1], the cyan and the blue [0, 0, 1], the blue and the magenta [1, 0, 1] and the magenta and the red.

A problem for the visualization is to find mapping of a point in the complex plane to a point in the cube. Clearly any simple function does not exist. One method for the mapping has been proposed in ref.6, where a distance from the original point [0, 0, 0] corresponds to the amplitude | ψ | and it uses most regions of the color cube that contain a point [1, 1, 1] that is the maximum white. There are two problems on this method. One is that it is difficult to distinguish the phases by colors for wave functions with large amplitudes because characters of colors near the white point are unclear even though they are bright. Also human recognitions of colors have some thresholds so that there might be a case that even for two colors, namely, ”a” and ”b” which have the same distance from the origin, one may recognize one color ”a” as a color while he may think the other color ”b” as the black. As a result, wave functions with small amplitudes can not give us the equal brightness for various phases.

In order to solve the above two problems, we suggest a hybrid use of two kinds of functions, one of which is suitable for the large amplitude and another is for the small one. For the large amplitude, we use only colors on a surface of a cone in the color cubes, where its vertices is [0, 0, 0], its center of the circle is [1/2, 1/2, 1/2] and the radius is \sqrt{6}/4. The circle passes through points [1, 1/4, 1/4], [1/4, 1, 1/4], [1/4, 1/4, 1], [3/4, 3/4, 0], [3/4, 0, 3/4] and [0, 3/4, 3/4](see Fig.7(b)) which will be ” red ”, ”green”, ”blue”, ”yellow”, ”magenta ” and ”cyan” on the cone surface, respectively. This method does not use colors near the white[1, 1, 1] and use the colors that are nearest to the maximum color. Then a distance from the origin corresponds to the amplitude while the angle in the circle does to the phase.

A point on the cone surface will be represented by three vectors.

n = [\frac{1}{2}, \frac{1}{2}, \frac{1}{2}], i = [\frac{1}{2}, −\frac{1}{4}, −\frac{1}{4}], j = [0, \frac{\sqrt{3}}{4}, −\frac{\sqrt{3}}{4}].

(4.1)

Here n denotes a position vector for a center of the cube and i and j are orthogonal vectors from the center to the cone circle. By using these vectors we make the amplitude and the phase to correspond with a point on the cone surface.

X_{c} = V_{c}{n + cos(H_{c})i + sin(H_{c})j}.

Here H_{c} and V_{c} are given by

H_{c} = arg(ψ), V_{c} = 1 − exp(−C|ψ|).

Here ”red” n + i = [1, 1/4, 1/4] on the cone surface is the nearest color to the red [1, 0, 0]. The constant C is fixed to correspond the maximum amplitude |ψ|_{max} with the maximum V_{c} within the color performance. If the display has N_{c} digital magnitude C = log(N_{c})/|ψ|_{max}.

For small amplitudes, the dark yellow will be recognized almost as black even if the dark red with the same amplitude is visible. In order to increase brightness of the dark yellow, dark cyan and dark magenta, we use only sides of the cube, which are lines between [V_{c}, 0, 0] and [0, V_{c}, 0], [0, V_{c}, 0] and [0, 0, V_{c}]and [0, 0, V_{c}] and [V_{c}, 0, 0].

\{ \begin{matrix}R_{2} = & V_{c} \\ G_{2} = & \frac{3}{π} H_{c}V_{c} \\ B_{2} = & 0 \end{matrix}  (0 ≤ H_{c} < π/3)

\{ \begin{matrix}R_{2} = & \frac{3}{π} (2π/3 − H_{c})V_{c} \\ G_{2} = & V_{c} \\ B_{2} = & 0 \end{matrix} (π/3 ≤ H_{c} < 2π/3)

\{ \begin{matrix}R_{2} = & 0 \\ G_{2} = & V_{c} \\ B_{2} = & \frac{3}{π} (H_{c} − 2π/3)V_{c} \end{matrix} (2π/3 ≤ H_{c} < π)

\{ \begin{matrix}R_{2} = & 0 \\ G_{2} = & \frac{3}{π} (4π/3 − H_{c})V_{c} \\ B_{2} = & V_{c} \end{matrix} (π ≤ H_{c} < 4π/3)

\{ \begin{matrix}R_{2} = & \frac{3}{π} (H_{c} − 4π/3)V_{c} \\ G_{2} = & 0 \\ B_{2} = & V_{c} \end{matrix} (4π/3 ≤ H_{c} < 5π/3)

\{ \begin{matrix}R_{2} = & V_{c} \\ G_{2} = & 0 \\ B_{2} = & \frac{3}{π} (2π − H_{c})V_{c} \end{matrix}  (5π/3 ≤ H_{c} < 2π)

(4.2)

This contribution is added to colors on the cone surface with a suppression factor (1 − V_{c}).

R = R_{1} + (1 − V_{c})R_{2},  (4.3)
G = G_{1} + (1 − V_{c})G_{2},  (4.4)
B = B_{1} + (1 − V_{c})B_{2}.  (4.5)

Based the above mapping, it is easy to represent our numerical results on the display if one uses available tools of computers even though some tedious and technical works are required. We have made animations for motions of two-dimensional wave packets for square-well potentials and those with spins under magnetic fields, the cyclotron motions and simulations of the Aharonov-Bohm effects, which are reported in ref.7, and have been recorded by video tapes^{8)}. Also one can see some of our results in the www server;

http://hep.esb.yamanashi.ac.jp/quantum/

§5. Conclusions and discussions

We studied the product formula by the Trotter-Suzuki decomposition of the operators, which is unconditionally stable, for numerical calculations of the time dependent Schrödinger equation and visualizations of these results. The product formula has been extended to cases for Hamiltonians with vector potentials. Our method can be applied to calculations with any gauge. By this extended formula one can study dynamical phenomena in the Aharonov-Bohm effect and others.

The extensive study of the numerical error for the time evolution of the one dimensional wave packets has confirmed us that the product formula is fully efficient. Through the study of the time evolutions under a constant magnetic field we have noticed that in available computer resource errors due to the finite time step Δt can be less than those due to the finite space difference Δx. In order to decrease the space difference much more memory of computers as well as the CPU time is required. This fact shows that it is quite important to improve the algorithms to decrease the errors due to Δx. Also we discussed the cyclotron motion in quantum mechanics. This study showed clearly the advantage of the unconditional stability in the product formula. Even for the rough approximation, the characteristic property, the circle motion of the wave packets, is obtained. This encourages us to make animations on many problems by this formula.

To visualize the numerical results we made the animations on the time evolution of the wave function where the phase is represented by colors and the absolute value is shown by brightness. We suggested the new way to represent a complex number by colors in order to solve the problems that the color with high brightness is too white and that the yellow with low brightness is too dark. Our method uses colors on surfaces of the cone and the cube in the color cube. Comparisons with the previous representation of the complex number can be seen in the www server.

We have made animations on motions of two-dimensional wave packets for square-well potentials and those with spins under magnetic fields, the cyclotron motions and simulations of the Aharonov-Bohm effects. We will discuss in other paper how much effect the visualizations have on our understanding contents of physics.

Acknowledgement We thank Dr. Y. Munehisa for her valuable comments on our work.

1) For example, see Mesoscopic Phenomena in Solids, ed. B. L. Altshuler, P. A. Lee and R. A. Webb: (North-Holland, Amsterdam, 1991).

2) C. Leforestier et al.: Comput. Phys. 94(1991)59-80.

3) H. De Raedt: Comput. Phys. Rep. 7(1987)1-72.

4) H. F. Trotter: Proc. Am. Math. Soc. 10(1959)545, M. Suzuki: Prog. Theor. Phys. 56(1976)1454.

5) T. Iitaka: Phys. Rev. E 5(1994)4684-4690.

6) J. L. Richardson: Comput. Phys. Commun. 63(1991)84-94.

7) H. Natori and T. Munehisa: Phys. Educ. in Univ. 3(1995)35 (Japanese).

8) Anyone can get the video tape on the animations by asking for authors.

Fig. 1. Norm and phase errors as a function of propagation time in the one dimensional free wave packets. Open and filled symbols are data for time step Δt = 0.01 and 0.001.

Fig. 2. Norm and phase errors for various values of the time step Δt. Every error is calculated at the same propagation time t = 100.0.

Fig. 3. Norm and phase errors as a function of propagation time in the Landau state. Open and filled symbols are data for time step Δt = 0.1 × 10^{−2}, 0.1 × 10^{−3} and L_{x} × L_{y} = 100 × 400.

Fig. 4. Norm and phase errors for various values of the time step Δt in the Landau state. Every error is calculated at the same propagation time t = 2π/E_{2} = 2.513. Data of Δx = Δy = 0.04 are given by L_{x} × L_{y} = 100 × 400 while data of Δx = Δy = 0.02 are calculated by L_{x} × L_{y} = 200 × 800.

Fig. 5. Orbits of average locations of the wave packet for the quantum cyclotron motion for Δt = 0.16×10^{−3} and Δt = 0.16× 10^{−4}. The initial location is (x, y) = (0, π). Data are plotted with the interval T_{c}/16 and data of t = T_{c} are shown by filled symbols.

Fig. 6. Differences between the average locations of the wave function after T_{c} and at the initial time for various values of the time step Δt. Every data is calculated at the same propagation time T_{c}.

Fig. 7. The color cube and the cone . (a) Any color can uniquely correspond to a point in the cube. The black, red, blue, green and white are points of [0,0,0],[0,0,1],[0,1,0],[0,0,1] and [1,1,1] respectively. (b) The center of the circle is [1/2,1/2/,1/2] which is pointed by n. The points R[1,1/4,1/4], Y[3/4,3/4,0], G[1/4,1,1/4], C[0,3/4,3/4], B[1/4,1/4,1] and M[3/4,0,3/4] are ”red”,”yellow”,”green”,”cyan”,”blue” and ”magenta” on the cone surface, respectively. A complex number is represented by a point on the surface of the cone side.

Table I. Values of parameters in calculations on one dimensional free wave packets. x_{min} and x_{max} are the minimum value and the maximum of x-space.

Parameter  Value
x_{min}  −50
x_{max}  49
Δx  1
p  π/5
σ_{x}  10
L  100

Table II. Values of parameters in calculations on the Landau state. We fixed values of x_{min} and x_{max} in order to satisfy that (x_{max} − x_{min})p_{x} = 2π. y_{min} and y_{max} are fixed under a

condition \frac{|eB|}{2h̄}(y_{min(max)} − \frac{p_{x}}{eB} )^{2} ≫ 1.

Parameter  Value
p_{x}  π/2
x_{min}, x_{max}  −2, 2
y_{min}, y_{max}  −8 − π/2, 8 − π/2
n  2
L_{x} × L_{y}  100 × 400, 200 × 800
B  1

Table III. Values of parameters in calculations on the quantum cyclotron motion.

Parameter  Value
p_{x}  π
p_{y}  0
σ_{x}  1
σ_{y}  1
x_{min}, x_{max}  −12, 12
y_{min}, y_{max}  −12, 12
x_{0}  0
y_{0}  π
L_{x} × L_{y}  600 × 600