Numerical solutions of the time-dependent Schrodinger equation in two dimensions

TL;DR

This paper presents a numerical approach using an adapted Crank-Nicolson method combined with high-order finite differences to solve the two-dimensional time-dependent Schrödinger equation efficiently and accurately.

Contribution

It introduces an adapted alternating-direction implicit method with high-order spatial discretization for improved numerical solutions of 2D Schrödinger equations.

Findings

The method achieves high accuracy with controlled step sizes.

It demonstrates efficiency in solving various 2D quantum systems.

The approach allows systematic analysis of precision and computational cost.

Abstract

The generalized Crank-Nicolson method is employed to obtain numerical solutions of the two-dimensional time-dependent Schrodinger equation. An adapted alternating-direction implicit method is used, along with a high-order finite difference scheme in space. Extra care has to be taken for the needed precision of the time development. The method permits a systematic study of the accuracy and efficiency in terms of powers of the spatial and temporal step sizes. To illustrate its utility the method is applied to several two-dimensional systems.

Figures

Tables

Table 1. Table 1: Summary of CPU time as r 𝑟 r is varied assuming a constant error for example 1. The fixed parameters are M = 5 𝑀 5 M=5 , N = 100 𝑁 100 N=100 , and Δ ​ t = 0.01 Δ 𝑡 0.01 \Delta t=0.01 .

$r$	$𝒥$	$h$	CPU(s)	$e_2\left({10}^{-7}\right)$
3	290	0.1379	897	8.27
4	200	0.1000	396	8.94
5	165	0.2424	317	9.97
6	150	0.2667	302	9.27
7	140	0.2587	306	9.63
8	135	0.2963	317	8.39
9	130	0.3077	331	9.96
10	130	0.3077	364	7.59
11	125	0.3200	353	8.27
12	123	0.3253	392	6.37
13	123	0.3253	424	4.98
14	123	0.3253	457	3.94
15	122	0.3279	466	9.75
16	122	0.3279	535	7.78
17	122	0.3279	577	7.14
18	120	0.3333	581	9.33
19	120	0.3333	629	8.62
20	120	0.3333	693	7.98

Table 2. Table 2: Summary of CPU time when M 𝑀 M is varied for example 1. The fixed parameters are r = 6 𝑟 6 r=6 , 𝒥 = 150 𝒥 150 {\cal J}=150 , and h = 0.2667 ℎ 0.2667 h=0.2667 .

$M$	$N$	$\Delta t$	CPU(s)	$e_2\left({10}^{-7}\right)$
1	2500	0.0004	398	10.00
2	130	0.0077	86.0	9.36
3	30	0.0333	76.8	9.74
4	25	0.0400	78.4	9.35
5	20	0.0500	92.6	9.27
6	18	0.0556	108	9.28
7	15	0.0667	128	9.28
8	14	0.0714	142	9.28
9	12	0.0833	153	9.27
10	11	0.0909	181	9.27
12	10	0.1000	215	9.27
15	8	0.1250	266	9.27
17	7	0.1429	303	9.27
20	6	0.1667	353	9.27

Table 3. Table 3: Parameters used for the error calculations of Figs.  4 and 5 .

$\hslash =m=1,n_x=n_y=0$	$k_x=k_y=0$
$x_0=y_0=-15,x_{𝒥}=y_{ℳ}=15$	${\alpha }_x={\beta }_x={\alpha }_y={\beta }_y=1$
$t_{\max }=2\pi ,N=100,dt=2\pi /N$	$A_x=A_y=2$
$𝒥=ℳ=100$

Table 4. Table 4: Parameters used for calculation of the oscillating, pulsating wave function shown in Fig.  6 .

$M=4,r=14,n_x=2,n_y=1$	$\hslash =m=1,k_x=0,k_y=5$
$x_0=y_0=-15,x_{𝒥}=y_{ℳ}=15$	${\alpha }_x=\sqrt{2},{\alpha }_y=1$
$t_{\max }=2\pi ,N=557,dt=2\pi /N$	${\beta }_x=2{\alpha }_x,{\beta }_y=2{\alpha }_y$
$𝒥=ℳ=200$	$A_x=-5,A_y=0$

Table 5. Table 5: Parameters used for the calculations of the free packet in Fig.  7 .

$\hslash =2m=1$,   $n_x=2,n_y=1$	$k_x=k_y=32$
$x_0=y_0=-2.5,x_{𝒥}=y_{ℳ}=2.5$	${\alpha }_x={\alpha }_y$
$N=80,dt={10}^{-4}$	$A_x=A_y=0.25$
$𝒥=ℳ=200,r=20,M=10$	${\alpha }_x=12.5\sqrt{2}$

Table 6. Table 6: Parameters used for the calculations of the wave packet passing through a single slit as shown in Figs.  8 and 9 .

$\hslash =2m=1$,   $n_x=n_y=0$	$k_x=32,k_y=0$
$x_0=y_0=-2,x_{𝒥}=y_{ℳ}=2$	${\alpha }_x={\alpha }_y$
$N=79,dt={10}^{-4}$	$A_x=-0.25,A_y=0$
$𝒥=ℳ=200,r=14,M=6$	${\alpha }_x=12.5\sqrt{2}$

Equations

$$\left(\hat{H}-i\hslash \frac{\partial }{\partial t}\right)\Psi (x,y,t)=0,$$

$$\begin{array}{rl}\hat{H}& =-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial y^2}+\hat{V}(x,y)\\ & ={\hat{K}}_x+{\hat{K}}_y+\hat{V}(x,y),\end{array}$$

$$\Psi (x,y,0)=\Phi (x,y).$$

$$\Psi (t+\Delta t)=e^{-i\hat{H}\Delta t/\hslash }\Psi (t),$$

$$e^{-i\hat{H}\Delta t/\hslash }=\prod _{s=1}^M{\hat{\Theta }}_s^{(M)}+O[(\Delta t{)}^{2M+1}],$$

$${\hat{\Theta }}_s^{(M)}\equiv \frac{1+(i\hat{H}\Delta t/\hslash )/z_s^{(M)}}{1-(i\hat{H}\Delta t/\hslash )/{\bar{z}}_s^{(M)}},$$

$${\Psi }^{(n+1)}=\prod _{s=1}^M{\hat{\Theta }}_s^{(M)}{\Psi }^{(n)}.$$

$${\Psi }^{+}=\left(\frac{1+(i\hat{H}\Delta t/\hslash )/z}{1-(i\hat{H}\Delta t/\hslash )/\bar{z}}\right){\Psi }^0,$$

$$\left(1-(i\hat{H}\Delta t/\hslash )/\bar{z}\right){\Psi }^{+}=\left(1+(i\hat{H}\Delta t/\hslash )/z\right){\Psi }^0,$$

$$\begin{array}{rl}& \left[1-i({\hat{K}}_x+{\hat{K}}_y+\hat{V})\frac{\Delta t}{\hslash \bar{z}}\right]{\Psi }^{+}\\ & =\left[1+i({\hat{K}}_x+{\hat{K}}_y+\hat{V})\frac{\Delta t}{\hslash z}\right]{\Psi }^0\end{array}$$

$$\begin{array}{rl}\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)& \left(1-i{\hat{K}}_y\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^{+}+{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^{+}\\ & =\left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)\left(1+i{\hat{K}}_y\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^0+{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^0+i\hat{V}\frac{\Delta t}{\hslash \bar{z}}{\Psi }^{+}+i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0.\end{array}$$

$$\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)\tilde{\Psi }=\left(1+i{\hat{K}}_y\frac{\Delta t}{\hslash z}\right){\Psi }^0+i\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0.$$

$$\begin{array}{rl}\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)\left(1-i{\hat{K}}_y\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^{+}=& \left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash z}\right)\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)\tilde{\Psi }-\left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash z}\right)\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0\\ & -{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^{+}+{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{\Psi }^0+iV\frac{\Delta t}{\hslash \bar{z}}{\Psi }^{+}+i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0.\end{array}$$

$$\begin{array}{rl}\left(1-i{\hat{K}}_y\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^{+}& =\left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash z}\right)\tilde{\Psi }-\left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash z}\right)i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0\\ & -{\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{-1}{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^{+}+{\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{-1}{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{\Psi }^0\\ & +{\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{-1}i\hat{V}\frac{\Delta t}{\hslash \bar{z}}{\Psi }^{+}+{\left(1-i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{-1}i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0.\end{array}$$

$$\begin{array}{rl}\left(1-i{K}_y\frac{\Delta t}{\hslash \bar{z}}-iV\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^{+}& =\left(1+i{K}_x\frac{\Delta t}{\hslash z}\right)\tilde{\Psi }+\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{K}_{x}V{\Psi }^0\\ & -\sum _{l=0}^{2(M-1)}{\left(i{K}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^l{K}_x{K}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^{+}+\sum _{l=0}^{2(M-1)}{\left(i{K}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^l{K}_x{K}_y\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{\Psi }^0\\ & +\sum _{l=1}^{2M-1}{\left(i{K}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{l}iV\frac{\Delta t}{\hslash \bar{z}}{\Psi }^{+}+\sum _{l=1}^{2M-1}{\left(i{K}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{l}iV\frac{\Delta t}{\hslash z}{\Psi }^0+O\left[(\Delta t{)}^{2M+1}\right].\end{array}$$

$$\left(1+i{\hat{K}}_y\frac{\Delta t}{2}+i\hat{V}\frac{\Delta t}{2}\right){\Psi }^{+}=\left(1-i{\hat{K}}_x\frac{\Delta t}{2}\right)\tilde{\Psi }.$$

$$\begin{array}{rl}\left(1-i{\hat{K}}_y\frac{\Delta t}{\hslash \bar{z}}-i\hat{V}\frac{\Delta t}{\hslash \bar{z}}\right){\Psi }^{\sigma +1}& =\left(1+i{\hat{K}}_x\frac{\Delta t}{\hslash z}\right)\tilde{\Psi }+\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{\hat{K}}_x\hat{V}{\Psi }^0\\ & -\sum _{l=0}^{2(M-1)}{\left(i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^l{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{\bar{z}}^2}{\Psi }^{\sigma }+\sum _{l=0}^{2(M-1)}{\left(i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^l{\hat{K}}_x{\hat{K}}_y\frac{(\Delta t{)}^2}{{\hslash }^2{z}^2}{\Psi }^0\\ & +\sum _{l=1}^{2M-1}{\left(i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{l}i\hat{V}\frac{\Delta t}{\hslash \bar{z}}{\Psi }^{\sigma }+\sum _{l=1}^{2M-1}{\left(i{\hat{K}}_x\frac{\Delta t}{\hslash \bar{z}}\right)}^{l}i\hat{V}\frac{\Delta t}{\hslash z}{\Psi }^0.\end{array}$$

$$f^{\prime \prime }(x)=\frac{1}{h^2}\sum _{k=-r}^{k=r}{c}_k^{(r)}f(x+kh)+O(h^{2r}),$$

$$\begin{array}{rl}{f}^{\prime \prime }(x)=& \sum _{\begin{array}{c}k=-r\\ (k\ne 0)\end{array}}^{k=r}{a}_k^{(r)}{f}^{\prime \prime }(x+kh)\\ & +\frac{1}{h^2}\sum _{k=-r}^{k=r}{b}_k^{(r)}f(x+kh)+O(h^{4r}).\end{array}$$

$$\begin{array}{rl}({\hat{K}}_x\Psi {)}_{j,m}=& \left(-\frac{{\hslash }^2}{2m}\right)\frac{1}{h_x^2}\sum _{k=-r}^{k=r}{c}_k^{(r)}{\Psi }_{j+k,m}\\ & \mathrm{for}0\le j+k\le \mathcal{J}.\end{array}$$

$$e=e^{(M)}+e^{(r)}=C_1(\Delta t{)}^{2M+1}+C_2{h}^{2r},$$

$$e_2^2={\int }_{x_0}^{x_{\mathcal{J}}}{\int }_{y_0}^{y_{\mathcal{M}}}|\Psi (x,y,T)-{\Psi }_{\mathrm{exact}}(x,y,T){|}^{2}dydx.$$

$$e_{\infty }=\underset{j,m}{\max }|{\Psi }_{j,m}-{\Psi }_{\mathrm{exact}}(x_j,y_m,t)|\mathrm{for}t=T.$$

$${\eta }_2^2={\int }_{x_0}^{x_{\mathcal{J}}}{\int }_{y_0}^{y_{\mathcal{M}}}|{\Psi }^{(M,r)}(x,y,T)-{\Psi }^{(M+1,r+1)}(x,y,T){|}^{2}dydx.$$

$$V(x,y)=-\frac{{\hslash }^2}{2m}(3-2{\tanh }^{2}x-2{\tanh }^{2}y)$$

$$\Psi (x,y,t)=\frac{i{e}^{it}}{2\cosh x\cosh y},$$

$$\mathrm{CPU}\mathrm{time}\propto {\left(e^{(r)}\right)}^{-1/(2r)}r.$$

$$e^{(M)}\propto (\Delta t/M^{\nu }{)}^{2M+1},$$

$$\mathrm{CPU}\mathrm{time}\propto N{M}^2\propto {\left(e^{(M)}\right)}^{-1/(2M+1)}{M}^{2+\nu }.$$

$$V(x,y)=\frac{1}{2}m({\omega }_x^2{x}^2+{\omega }_y^2{y}^2).$$

Full text

I Introduction

The determination of accurate numerical solutions of the time-dependent Schrödinger equation is an ongoing enterprise. The quantum wave equation is fundamental to the understanding of nonrelativistic atomic and subatomic systems and phenomena. Consequently it occurs in a diversity of physical systems. Ideally analytic solutions are available, but most realistic situations are too complex to yield such solutions.

In the last few years a number of improvements have been made to yield more accurate solutions with greater efficiency. The type of method often depends on the problem at hand, i.e., dimensionality, time dependence of the interaction, short- or long-time behaviour, etc. The “method of choice” for some years is the Chebyshev polynomial expansion of the time-evolution operator with (inverse) Fourier transformations to deal with the spatial development as time progresses Tal-Ezer and Kosloff (1984); Leforestier et al. (1991). More recently the Padé approximant representation of the time-evolution operator is exploited Puzynin et al. (1999); *puzynin00; van Dijk and Toyama (2007); van Dijk et al. (2011); van Dijk and Toyama (2014); van Dijk (2016). This approach is unitary, stable, and allows for systematic estimate of errors in terms of powers of the temporal and spatial step sizes. The two approaches have been shown to have comparable efficacy Formánek et al. (2010); van Dijk et al. (2011). Gusev et al. Gusev et al. (2014); Chuluunbaatar et al. (2008a, b) have recently given an improved and extended application of the method discussed by Puzynin et al. Puzynin et al. (1999); *puzynin00. They deal with the more general problem of a time-dependent Hamiltonian. Using a truncated Magnus expansion with additional transformations, they are able to obtain stable and efficient solutions which are accurate up to sixth-order in the time step.

Generally the various approaches involve time evolution and integration over space. Thus there are a number of ways of dealing with the time evolution. Crank-Nicolson approximates the exponential time-evolution operator by a Cayley form which retains unitarity, but is correct only to low order in time advance Crank and Nicolson (1947). The Chebyshev polynomial expansion can lead to high accuracy even over significant time intervals. It is not explicitly unitary. The generalized Crank-Nicolson approximates the evolution operator with a $[M/M]$ Padé approximant, factorized into $M$ factors of Cayley form. This form is unitary and has a truncation error of $O[(\Delta t)^{2M+1}]$, where $\Delta t$ is the temporal step size. This improves the precision rapidly with increasing $M$. Besides these three approaches there are other approximations of the time-evolution operator, e.g., the exponential split-operator method Feit et al. (1982) or the iterative Lanczos reduction Park and Light (1986). Like time development the spatial integration can be achieved in different ways, e.g., by different types of finite differencing or by the pseudospectral fast Fourier transform approach.

Since many of the calculations referred to have been done in one spatial dimension, in this paper we consider the generalized Crank-Nicolson with two spatial dimensions. A number of articles have appeared recently that describe methods of solving the two-dimensional time-dependent Schrödinger equation, including those with time-dependent potentials and nonlinear terms. See, for example, Refs. Tian and Yu (2010); Xu and Zhang (2012); Gao and Mei (2016); Symes et al. (2016); Wang et al. (2016); Zhang and Chen (2016). A number of these use the Cayley form for the time evolution operator. We wish to employ the higher order Padé form in order to enhance the efficiency of the approach. Given the two spatial dimensions, we pursue an alternating-direction implicit scheme which requires only solving one-dimensional implicit problems for each time step. Different approaches have been suggested, such as the use of multigrid partitioning Gaspar et al. (2014), but it is our intention to present one that provides the user with another efficient alternative. Clearly the method chosen will depend on the context.

In section II we formulate the time dependence of the problem. Section III is a description of the spatial integration. A number of applications are discussed in Sec. IV, and Sec. V presents conclusions and a discussion of the work.

II Accurate time-evolution scheme

We solve the two-dimensional time-dependent Schrödinger equation

[TABLE]

where

[TABLE]

starting with an initial wave function

[TABLE]

The time-evolution operator of the system gives an expression for the wave function at a time in terms of the wave function at an earlier time, i.e.,

[TABLE]

where $\Delta t$ is the time advance, and where we have suppressed the spatial coordinates $x$ and $y$ in the wave function. We will employ the factorized $[M/M]$ Padé approximant along with the alternating-direction implicit method Peaceman and H. H. Rachford (1955). In keeping with the expansion of the time-evolution operator discussed in Ref. van Dijk and Toyama (2007), the operator is written as

[TABLE]

where

[TABLE]

and $z_{s}^{(M)},s=1,\dots,M$ are the roots of the numerator of the $[M/M]$ Padé approximant of $e^{\textstyle z}$; the $\bar{z}_{s}^{(M)}$ are the corresponding complex conjugates. Since $\Psi^{(n+1)}=e^{-i\widehat{H}\Delta t/\hbar}\Psi^{(n)}$ ($n$ refers to the time $t_{n}=n\Delta t,n=0,1,\dots$), we write

[TABLE]

Defining $\Psi^{(n+s/M)}\equiv\widehat{\Theta}_{s}^{(M)}\Psi^{(n+(s-1)/M)}$, we can solve for $\Psi^{(n+1)}$ iteratively starting with $\Psi^{(n+1/M)}=\widehat{\Theta}_{1}^{(M)}\Psi^{(n)}$, then $\Psi^{(n+2/M)}=\widehat{\Theta}_{2}^{(M)}\Psi^{(n+1/M)}$, and so on.

Let us start with the basic substep of the procedure in going from $\Psi^{(n+(s-1)/M)}$ to $\Psi^{(n+s/M)}$, which we label below generically as $\Psi^{0}$ and $\Psi^{+}$, respectively. We then write

[TABLE]

or

[TABLE]

where $z$ is the generic $z^{(M)}_{s}$. Since $\widehat{H}=\widehat{K}_{x}+\widehat{K}_{y}+\widehat{V}$, we write

[TABLE]

so that

[TABLE]

In keeping with Peaceman and Rachford Peaceman and H. H. Rachford (1955), we define $\widetilde{\Psi}$ by the equation

[TABLE]

We insert this expression into Eq. (2.11) to obtain

[TABLE]

Operating on Eq. (2.13) with the inverse of $\left(1-i\widehat{K}_{x}\dfrac{\Delta t}{\hbar\bar{z}}\right)$, we get

[TABLE]

The inverse operators are expanded, but we must make sure that the expansions are correct to $O[(\Delta t)^{2M}]$ since the overall expansion (2.5) is of that order.

Simplifying and keeping terms up to $(\Delta t)^{2M}$ and assuming $M>1$, we obtain

[TABLE]

The $M=1$ case, for which $z=-2$, results in the equation

[TABLE]

This equation is a typical implicit equation with the Cayley form. We solve Eq. (2.15) iteratively so that $\Psi^{\sigma+1}\rightarrow\Psi^{+}$ as $\sigma=0,1,2,\dots$ increases in the equation

[TABLE]

We start the iteration with setting $\Psi^{\sigma=0}=\Psi^{0}$. When $\Psi^{\sigma+1}$ and $\Psi^{\sigma}$ are sufficiently close we stop. Note that we need to calculate $\widetilde{\Psi}$ only once for each sequence of iterations. We find that this approach can give accurate results; typically around six iterations are required for precise results. This process has to be repeated for each of the $M$ steps needed to achieve a full time step advance.

There is an alternative approach to solving Eq. (2.15) for $\Psi^{+}$. The terms on the right side involving $\Psi^{+}$ can be moved to the left side and one solves a linear system of equations upon the discretization of the spatial variables. However, as we show in the next section, the kinetic energy operators are banded diagonal matrices, and those operators raised to some power would result in matrices with the size of the bands increased. As a result the gains in efficiency of a banded matrix formulation are lost.

III Spatial integration

The numerical spatial integration of the partial differential equation (2.17) can be done in a number of ways. Two approaches often considered are the spectral decomposition of the spatial (kinetic energy) operator or the finite-difference representation of this operator. The relative merits are discussed by the authors of Ref. Cerjan and Kulander (1991). They point out that a“low-order differencing method is in principle faster than a spectral method since it scales as the ’bandedness’ times the size of the grid, $O(bN)$, rather than as $O(N\log N)$”. In the case of two-dimensional systems using the alternating-direction implicit approach $N$ is replaced by $N^{2}$, whereas $b$ is unchanged. For the purpose of this work we therefore use finite differences. One could choose the traditional three-point expression for the second-order partial derivative. There are however more precise methods. For instance the recent Numerov recent approach van Dijk (2016) gives much higher accuracy, as does the high-order compact finite difference approach in Refs. Tian and Yu (2010); Xu and Zhang (2012). The traditional approach is $O(h^{2})$, where $h$ is the spatial step size, whereas the high-order compact method is $O(h^{4})$, and the Numerov algorithm is $O(h^{5})$. The advantage of these approaches is that they lead to three-point formulas which may be convenient when crossing a discontinuity of the potential or considering an adaptive spatial grid van Dijk (2016).

As in earlier work van Dijk and Toyama (2007) we consider formulas which allows one to choose an arbitrary order of $h$. For a spatial grid (in one dimension) with step size $h$, the second derivative of $f(x)$ is expanded as

[TABLE]

where the $c_{k}^{(r)}$ are real constants, obtained from making series expansions of the functions $f(x\pm kh)$. A similar technique is used by Wang and Shao for the kinetic energy operator acting on the wave function of a two-dimensional stationary state problem Wang and Shao (2009). In another article the same authors suggest an expansion of the form Shao and Wang (2009)

[TABLE]

In one dimension the discretized kinetic energy is expressed as a banded diagonal matrix with bandwidth of $2r+1$, just like in the case of Eq. (3.1). Thus it seems that with virtually the same effort the calculation gives much more accurate results. A comparison of the two expansions van Dijk et al. (2011) shows that for smaller values of $r$ the calculation is indeed much more efficient, however for larger $r$ the accuracy decreases. The kinetic energy operator resulting from Eq. (3.1) can be made strictly diagonally dominant, whereas the diagonal dominance of the kinetic energy matrix from Eq. (3.2) becomes compromised when $r$ goes beyond ten. In this paper we use expansion (3.1) for the kinetic energy operators $\widehat{K}_{x}$ and $\widehat{K}_{y}$.

We consider a rectangular domain in space $[x_{0},x_{\mathscr{J}}]\times[y_{0},y_{\mathscr{M}}]\subset\mathbb{R}^{2}$, which we partition uniformly in each direction, so that with $h_{x}=(x_{\mathscr{J}}-x_{0})/\mathscr{J}$ and $h_{y}=(y_{\mathscr{M}}-y_{0})/\mathscr{M}$, $x_{j}=x_{0}+jh_{x},j=0,\dots,\mathscr{J}$ and $y_{m}=y_{0}+mh_{y},m=0,\dots,\mathscr{M}$. The time is also partitioned over the time interval from 0 to $T$ into $N$ subintervals, so that $\Delta t=T/N$ and the intermediate times are $t_{n}=n\Delta t$, where $n=0,1,\dots,N$. The equations we need to solve are typically of the type Eqs. (2.12) and (2.17). If we let $\Psi(x_{j},y_{m})=\Psi_{j,m}$, then $(\widehat{V}\Psi)_{j,m}=V_{j,m}\Psi_{j,m}$ and

[TABLE]

There is a similar relation for $\widehat{K}_{y}\Psi$ except that the summation is over the second index of $\Psi_{j,m+k}$. Thus in Eq. (2.17), for example, the right side is completely specified, but the $\Psi^{\sigma+1}$ on the left side needs to be found. This equation is really a linear system of equations with a banded diagonal coefficient matrix over the index $m$. It can be solved for each $j$ to obtain $\Psi^{\sigma+1}$. It is the strength of the alternating-direction implicit scheme that calculations are reduced to one-dimensional ones.

IV Implementation

In this section we consider four examples in which the method outlined previously is applied in order to investigate its accuracy and efficiency. We will also demonstrate the feasibility of calculating wave functions with more complex structure (several peaks and valleys) as they evolve in time.

IV.1 Errors

The truncation error of the series expansion of the wave function in time and space can be expressed as

[TABLE]

where we are considering the error of the real quantity $|\Psi(x,y,t)|$ and $C_{1}$ and $C_{2}$ are real positive numbers related to the $(2M+1)$th partial derivative with respect to $t$ and the $(2r)$th derivative with respect to $x$ or $y$, respectively. For simplicity we assume $h_{x}=h_{y}=h$ and a $xy$ symmetry of the wave function. If exact analytic solutions are available the error can be calculated by comparison. If that is not the case, a good estimate of the error can be made by comparing the solution for particular $M$ and $r$ to the one obtained when one or both of the $M$ and $r$ are increased by unity van Dijk et al. (2011).

To make comparisons of the numerically obtained solutions to analytic solutions in cases where the latter are known, we define the error $e_{2}$ such that

[TABLE]

The error $e_{2}$, which is a Euclidean/$\ell_{2}$-vector norm, is a measure of the accuracy of the wave function and its phase. Alternatively some authors have used the $\ell_{\infty}$-vector norm

[TABLE]

In the case that no exact solution is available, one can make an estimate of the error by comparing a solution obtained with particular values of $M$ and $r$ to the solution obtained with $M+1$ and $r+1$, e.g.,

[TABLE]

Since it turns out that $e_{2}$ and $\eta_{2}$ are very similar, evaluating $\eta_{2}$ provides a method to estimate the accuracy in the absence of an analytic solution van Dijk et al. (2011); van Dijk (2016).

Since in our applications the wave functions are zero near the boundary of the domain, we can use the simple rectangle rule for integration. The corrections to higher order polynomial approximations are all in terms of evaluations of the integrand near the end points, but since the wave function is zero there, one gets very accurate integrals with the simple quadrature Peters and Maley (1968).

IV.2 Example 1: solvable two-dimensional potential

Consider the potential

[TABLE]

with $\hbar=2m=1$, see Fig. 1.

This potential has been used by several authors as one which tests numerical methods $\breve{\rm C}$iegis et al. (2013); Tian and Yu (2010); Xu and Zhang (2012). We choose the domain $[-20,20]\times[-20,20]$ and solve the problem from $t=0$ to $t=T=1$ with the solution on the boundary equal to zero $\breve{\rm C}$iegis et al. (2013). The exact solution is

[TABLE]

which is also used to determine the initial wave function. It should be noted that wave function (4.6) is square integrable and is an energy eigenstate with energy $-\hbar^{2}/(2m)$. It describes a bound state at threshold; the energy spectrum at higher energies is a continuum and the corresponding wave functions are unbound.

In our numerical calculation we allow $M=1,\dots,6$ and $r=1,\dots,20$, with $\Delta t=0.01$ (or 100 time steps) and with ${\cal J}={\cal M}=200$. The $e_{2}$ are plotted as a function of $r$ for various values of $M$ in Fig. 2.

We compare this calculation to that of Ref. $\breve{\rm C}$iegis et al. (2013) since the other two calculations Tian and Yu (2010); Xu and Zhang (2012) are done over a much smaller spatial domain with much smaller time intervals. In our calculation we find generally that $e_{2}\gtrsim e_{\infty}$. In Ref. $\breve{\rm C}$iegis et al. (2013) the quoted errors are approximately $e_{\infty}\approx 10^{-4}$. Figure 2 shows that the error in the calculation is reduced significantly when one goes from $M=1$ to $M=2$ and 3. When $M>3$ the results are identical to those of $M=3$. Given that earlier calculations referred to are basically $M=1$ calculations with spatial errors of the order of $h^{3}$ or $h^{5}$, this method results in significant improvement in accuracy.

The errors of this example for $M\geq 3$ saturate at $e_{2}\approx 1.5\times 10^{-9}$. In the “gullies” of potential (4.5) the magnitude of the wave function is larger than elsewhere. In the gullies at the boundary of the computational space it is approximately $10^{-9}$. The numerical calculation assumes that the wave function is zero outside the computational domain. The discrepancy between the numerical wave function and the exact one outside the computational space is the source of the residual error.

In order to determine the efficiency of the calculation we obtain the CPU time when the error is close to, but less than, $10^{-6}$. For a particular $M$ we adjust the spatial step size by choosing $r$ until we reach a minimum CPU time. Similarly for a particular $r$, we chose $M$ to yield minimum CPU time. The results are listed in Tables 1 and 2. In Ref. van Dijk and Toyama (2007) we give an estimate for the CPU time as a function of $r$ for the one-dimensional calculation. In two-dimensions we expect the behaviour to be similar since errors in $x$ and $y$ integrations are similar and additive, especially in symmetric cases. Thus,

[TABLE]

To obtain a formula for the CPU time as a function of $M$ when $r$ and the error are constant, we take the relation of the error and $M$ to be

[TABLE]

where $\nu$ is a number less than unity. In the one-dimensional case we showed that $\nu\approx 1$, but Table 2 shows that a better relationship, especially for larger $M$, has $\nu<1$. Since $T=N\Delta t$ is a constant we write

[TABLE]

The CPU times as a function of $r$ and $M$ are shown in Fig. 3.

The estimates of the errors are shown as solid lines. For the CPU time as a function of $M$ we have estimated $\nu=1/2$. Such an estimate seems reasonable in light of the fact that the iterative part of the procedure increases the time, and the number of iterations vary with the value of $M$.

IV.3 Example 2: oscillating and pulsation harmonic oscillator wave functions

Consider the potential function for the two-dimensional anisotropic harmonic oscillator,

[TABLE]

An analytic solution for such a potential is van Dijk et al. (2014)

[TABLE]

where

[TABLE]

The various quantities in Eq. (4.12) are defined as follows:

[TABLE]

The wave function (4.12) is the pulsating and oscillating wave function of a particle subject to a one-dimensional harmonic oscillator characterized by $\omega$ or $\alpha$. The initial ($t=0$) wave function is the $n$th energy state of the particle subject to an oscillator characterized by $\beta$, rather than $\alpha$, displaced from the origin by amount $A$ and with a momentum $\hbar k$. The function $H_{n}(\xi)$ is the $n$th-order Hermite polynomial. This wave function provides a wave packet with more fluctuation than the traditional coherent wave packet; for instance, it has nodes which travel with the packet and whose occurrence spread and contract in time.

As an initial study of the accuracy of the method we consider the simplest case of an isotropic oscillator with the initial state the ground state. The values of the parameters are give in Table 3.

In Fig. 4 the error as a function of $M$, the order of the diagonal Padé approximant, is displayed. Note that the lower values of $M$, especially for larger $r$ give no results because the convergence of the iterative part of the calculation is not achieved. By decreasing $dt$ convergence can again be attained, but in Fig. 4 we keep $dt$ constant . The horizontal plateaux are not completed since there is no change in the error as $M$ is further increased. The calculations are done with double-precision floating-point arithmetic. We achieve an error less than $10^{-11}$ for $r=30$ and $M\geq 6$. The errors could be further reduced by increased computational precision.

For the same model we plot the error as a function of $r$ in Fig. 5 for a number of values of $M$. The lower pattern of dots is one that is obtained for each $M$ value up to a particular value of $r$, at which the error is constant as $r$ is increased further. These horizontal plateaux only extend to a certain point after which the iterative procedure becomes unstable. The plateaux are clearly visible for $M=4$ and 5, and the beginnings can be discerned for $M=1,2,3$. The instability of the calculation does not mean that we cannot obtain results in those regions. In this calculation $dt$ is the same in all cases; where instability sets in a smaller $dt$ will restore stability. The fact that the curves superimpose on the left can be seen from Fig. 4 where for each $M\geq 6$ the errors converge for $r$ sufficiently large.

On Fig. 5 we have plotted a solid line which is an estimate of the error obtained by considering the truncation error of the expansion in $x$ or $y$, i.e.,

[TABLE]

where $\xi$ is some value in the domain of $x$. We have assumed that $x$ dependence of the wave function of this example is a Gaussian and that the maximum value of it and its even-order derivatives occur when the argument is zero (Abramowitz and Stegun, 1965, p. 933), i.e.,

[TABLE]

The last factor on the right side of Eq. (4.14) is an adjustment to give reasonable agreement with the data. It amounts to an effective spatial step size which is smaller by a factor of 2.5. The shape of the solid curve is very sensitive to the form of this factor.

We plot the progression of the oscillating and pulsating wave packet as numerically determined in Fig. 6. The parameters used for this calculation are listed in Table 4. The error $e_{2}$ in the calculation

ranges from $3\times 10^{-7}$ after one time step to $5\times 10^{-3}$ after 557 steps.

Note that the oscillating frequency is $\omega_{x}=2\omega_{y}$. This means that when the motion has executed a complete cycle in the $x$ direction it has only gone through half a cycle in the $y$ direction. Thus the packet starts at (-5,0), travels along a quarter-elliptical path to (5,5), then back to (-5,0), to (5,-5), and to the initial point (-5,0). The pulsating frequencies in each direction are four times the corresponding oscillating frequencies.

IV.4 Example 3: free wave packet

For the free wave packet we consider the Hermite-Gaussian wave function of Ref. van Dijk et al. (2014),

[TABLE]

where

[TABLE]

with

[TABLE]

The travelling wave packet will have nodes whose distribution, if there is more than one node, spread in time. The model is similar to that of Galbraith et al. Galbraith et al. (1984) in whose calculation $n_{x}=n_{y}=0$ and $\alpha_{x}=\alpha_{y}$. The parameters we use are given in Table 5.

The free wave packet at times $t=10^{-3}$ and $t=6.3\times 10^{-3}$ is shown in Fig. 7. The separation of the peaks of the wave function as time progresses is clearly evident.

IV.5 Example 4: single-slit diffraction

The wave nature of electrons has been studied and observed in semiconductor nanostructures. Endoh et al. Endoh et al. (1992); *endoh99 have considered numerical simulations of the passage of such electrons through narrow constrictions. Recent experiments observed controlled electron diffraction for both single- and double-slit configurations Bach et al. (2013); Khatua et al. (2014).

The single slit in the barrier is obtained by introducing a potential

[TABLE]

where $f(x)$ is the difference of two Fermi functions

[TABLE]

In the calculation we choose $\mu=100$, $x_{w}=0.05$, and $V_{0}=1000$. Initially the free wave function (4.16) (with $n_{x}=n_{y}=0$ and $\alpha_{x}=\alpha_{y}=12.5\sqrt{2}$) impinges on the slit and diffracts. The parameters of the single-slit calculation are given in Table 6.

We plot the probability density $|\Psi(x_{0},y,t_{0})|^{2}$ as a function of $y$ when $x_{0}=1.253333$ and $t_{0}=0.0057$ in Fig. 8. The graph has a remarkable similarity to the Fraunhofer diffraction intensity. The slit is not of uniform width and hence we cannot compare parameters. The calculation does indicate that one can study different slit configurations and shapes Endoh et al. (1992); *endoh99 using this method.

In Fig. 9 we plot two snapshots of the wave packet passing through the slit. Since most of the packet is reflected, we multiply the amplitude of the diffracted packet in the figure by ten in order to make the packet’s shape in the region beyond the slit more visible.

V Discussion

We have shown that the generalized Crank-Nicolson method combined with the alternating-direction implicit procedure is a practical approach to the determination of numerical solutions of the two-dimensional Schrödinger equation. The method allows one to study the efficiency and accuracy systematically as a functions of powers of the temporal and spatial step sizes. Equation (2.15) is basic to the solution. It can be solved in different ways, but we choose to use an iterative approach which means one must find solutions of linear systems of equations whose coefficients form banded diagonal matrices. Since the number of iterations is low, the alternative noniterative approach leads to less sparse matrices and a correspondingly less efficient procedure.

A number of authors Tian and Yu (2010); Gao and Xie (2015); Xu and Zhang (2012) have considered alternating direction implicit compact finite difference schemes which give errors of order $O(h^{5}+\Delta t^{3})$. Although they include nonlinear equations in their analysis, they discuss, among others, example 1 of this paper as a test case. Our scheme gives errors $O(h^{2r}+\Delta t^{2M+1})$ where $M$ and $r$ are positive integers.

The examples demonstrate that this method is capable of accurate solutions even when there is significant fluctuation of the wave function. The methods described in this paper allow one to obtain a realistic theoretical analysis of the diffraction experiments that have been done recently. Given the recent attention to the two-dimensional nonlinear Schrödinger equation and the Pitaevskii equation, a future project is to expand the method of this paper to such systems, as well as those with time-dependent interactions or source terms. This in effect is a generalization of the work done earlier on one-dimensional systems van Dijk et al. (2014).

Furthermore it remains to systematically investigate the relative efficiency and accuracy of the approach of this paper to other methods that have been used or proposed. A generalization to three or higher spatial dimensions and the introduction of transparent boundary conditions are further natural extensions of this work.

Acknowledgements.

We are grateful to the Natural Sciences and Engineering Research Council of Canada for the Undergraduate Student Research Awards given to S.-J.P. (2013) and T.V. (2016).