Nonstandard Finite Difference Time Domain (NSFDTD) Method for Solving the Schr\"odinger Equation

TL;DR

This paper introduces a nonstandard finite difference time domain (NSFDTD) method that improves the accuracy of solving the Schrödinger equation by iteratively refining eigen-energy estimates using a modified numerical scheme.

Contribution

The paper presents a novel NSFDTD approach that enhances eigen-energy calculations for quantum systems by replacing standard derivatives with non-standard schemes and iterative refinement.

Findings

Validated on infinite square well, harmonic oscillator, and Morse potentials.

Achieves more accurate eigen-energies compared to standard FDTD.

Demonstrates iterative improvement in numerical solutions.

Abstract

In this paper, an improvement of the finite difference time domain (FDTD) method using a non-standard finite difference scheme is presented. The standard numerical scheme for the second derivative in the spatial domain is replaced by a non-standard numerical scheme. In order to apply the non-standard FDTD (NSFDTD), first estimates of eigen-energies of a system are needed and computed by the standard FDTD method. These first eigen-energies are then used by the NSFDTD method to obtain improved eigen-energies. The NSFDTD method can be performed iteratively using the resulting eigen-energies to obtain accurate results. In this paper, the NS-FDTD method is validated for infinite square well, harmonic oscillator and Morse potentials.

Figures

Tables

Table 1. Table 1: Numerical eigen-energies for a particle in an infinite square well potential computed by the FDTD method (standard scheme), the NSFDTD method after first (1st), second (2nd), third (3rd) and final iteration are compared with exact values. The NSFDTD parameters used are Δ ​ x = 0.1 Δ 𝑥 0.1 \Delta x=0.1 dan the length of the square well or the computational domain is 1. The correct numbers are in bold.

n	FDTD	NSFDTD				Exact
		1st	2nd	3rd	final
1	4.894348	4.934469	4.934799	4.934802	4.934802	4.934802
2	19.098301	19.717997	19.738506	19.739186	19.739209	19.739209
3	41.221475	44.174164	44.395261	44.411870	44.413220	44.413220
4	69.098301	77.637465	78.778707	78.932758	78.956835	78.956835
5	100.000000	118.475509	122.325053	123.146042	123.370050	123.370055

Table 2. Table 2: Numerical eigen-energies for a particle in an harmonic oscillator potential computed by the FDTD method (standard scheme), the NSFDTD after first (1st), second (2nd), third (3rd) and final iteration are compared with the exact values. The length of the computational domains are 16 and 10 for Δ ​ x = 0.8 Δ 𝑥 0.8 \Delta x=0.8 and Δ ​ x = 0.5 Δ 𝑥 0.5 \Delta x=0.5 respectively.

No	FDTD	NSFDTD				Exact
		1st	2nd	3rd	final
Spatial Interval $\Delta x=0.8$
1	0.479077	0.497481	0.498004	0.498019	0.498019	0.5
2	1.391838	1.483898	1.491834	1.492521	1.492586	1.5
3	2.188108	2.435512	2.471583	2.476871	2.477781	2.5
4	2.954963	3.357062	3.439530	3.456553	3.460975	3.5
5	3.236360	4.050845	4.306467	4.384808	4.418670	4.5
Spatial Interval $\Delta x=0.5$
1	0.492062	0.499624	0.499705	0.499706	0.499706	0.5
2	1.459805	1.497690	1.498910	1.498949	1.498951	1.5
3	2.393936	2.491568	2.496825	2.497108	2.497125	2.5
4	3.292528	3.479435	3.493647	3.494727	3.494816	3.5
5	4.153118	4.459088	4.489421	4.492424	4.492753	4.5

Table 3. Table 3: Same as Table 2 except that this is for the sum of squared errors of wave functions.

No	FDTD	NSFDTD
		1st	2nd	3rd	final
Spatial Interval $\Delta x=0.8$
1	0.002640	0.001655	0.001631	0.001630	0.001630
2	0.015449	0.004473	0.003858	0.003807	0.003802
3	0.071309	0.011986	0.008113	0.007629	0.007548
4	0.134618	0.022206	0.011478	0.009849	0.009463
5	0.662104	0.118040	0.031514	0.017526	0.013284
Spatial Interval $\Delta x=0.5$
1	0.000611	0.000381	0.000379	0.000379	0.000379
2	0.003254	0.000878	0.000827	0.000826	0.000826
3	0.011101	0.001747	0.001501	0.001489	0.001488
4	0.029720	0.002703	0.001972	0.001925	0.001921
5	0.068005	0.003694	0.002016	0.001898	0.001886

Table 4. Table 4: Numerical eigen-energies for a particle in a Morse potential computed by the FDTD method (standard scheme), the NSFDTD after first (1st), second (2nd), third (3rd) and final iteration are compared with the exact values. The parameters used are Δ ​ x = 0.2 Δ 𝑥 0.2 \Delta x=0.2 , x 0 = 2 subscript 𝑥 0 2 x_{0}=2 and the computational domain is at an interval [0,8].

No	FDTD	NSFDTD				Exact
		1st	2nd	3rd	final
1	4.754956	4.865394	4.867215	4.867245	4.867246	4.875
2	13.365245	13.830904	13.851604	13.852525	13.852568	13.875
3	20.758056	21.762245	21.828617	21.832999	21.833308	21.875
4	27.087090	28.686514	28.818009	28.828801	28.829763	28.875
5	32.484206	34.633288	34.831976	34.850315	34.852177	34.875

Equations

$$\frac{\partial \psi (x,t)}{\partial t}=-\hat{H}\psi (x,t)=-\left[-\frac{1}{2}\frac{{\partial }^2}{\partial x^2}+V(x)\right]\psi (x,t)$$

$$\frac{{\partial }^2\psi (x,t)}{\partial x^2}+W(x)\psi (x,t)=0$$

$$E=\frac{\int {\psi }^{\ast }(x)\hat{H}\psi (x)dx}{\int |\psi (x){|}^{2}dx}$$

$$\frac{{\partial }^2\psi (x)}{\partial x^2}\approx \frac{\psi (i+1)-2\psi (i)+\psi (i-1)}{(\Delta x{)}^2}$$

$$\psi (i+1)-2\psi (i)+\psi (i-1)=$$

$$2\left[\frac{(\Delta x{)}^2}{2!}{\psi }^{\prime \prime }+\frac{(\Delta x{)}^4}{4!}{\psi }^{(4)}+\frac{(\Delta x{)}^6}{6!}{\psi }^{(6)}+\dots \right]$$

$$\psi (i+1)-2\psi (i)+\psi (i-1)=$$

$$2\frac{{\psi }^{\prime \prime }}{W}\left[\frac{(\Delta x\sqrt{W}{)}^2}{2!}-\frac{(\Delta x\sqrt{W}{)}^4}{4!}+\frac{(\Delta x\sqrt{W}{)}^6}{6!}-\dots \right]$$

$$\psi (i+1)-2\psi (i)+\psi (i-1)={\psi }^{\prime \prime }\left[\frac{4}{W}{\sin }^2\left(\sqrt{W}\Delta x/2\right)\right]$$

$$\frac{{\partial }^2\psi (x)}{\partial x^2}\approx \frac{\psi (i+1)-2\psi (i)+\psi (i-1)}{g(i)}$$

$$g(x)=\left\{\begin{array}[]{ll}(\Delta x)^{2}&\text{for}\ W(x)=0\\ \frac{4}{W(x)}\sin^{2}\left(\sqrt{W(x)}\Delta x/2\right)&\text{for}\ W(x)>0\\ \frac{4}{-W(x)}\sinh^{2}\left(\sqrt{-W(x)}\Delta x/2\right)&\text{for}\ W(x)<0\end{array}\right.$$

$$g(x)\approx (\Delta x{)}^2-\frac{W(x)}{12}(\Delta x{)}^4+\frac{[W(x){]}^2}{360}(\Delta x{)}^6$$

$${\psi }^{n+1}(i)=\alpha (i){\psi }^n(i)+\frac{\beta (i)\Delta t}{2g(i)}[{\psi }^n(i+1)-2{\psi }^n(i)+{\psi }^n(i-1)]$$

$$E$$

$$V(x)=D[1-\exp (-(x-x_0)){]}^2$$

Full text

\corres

wayan.sudiarta@unram.ac.id \msinfo3 January 20183 January 20183 January 2018

\doinum

??/?? \artcitid#### \volnum? 2018 \pgrange1–7

\lp7

1 Introduction

The finite difference time domain (FDTD) method has been used extensively to simulate various electromagnetic wave radiations and interactions [1, 2]. Because of its usefulness, easy in programming, good accuracy, and flexibility, the FDTD method has been applied to other fields such as sound wave propagations [3] and simulations of quantum systems [4, 5].

To simulate a quantum system in real time using the FDTD method, one can use split method where the time-dependent Schrödinger equation (TDSE) is splitted into two equations for real and imaginary parts [6] or one can also directly use complex computation without splitting the TDSE [7]. Eigen values and eigen functions can not be directly obtained from the real-time FDTD method. It requires a Fourier transformation procedure to obtain eigenvalues and an additional FDTD simulation to compute the eigen functions. Beside that for infinite spatial domain, the FDTD method requires also an absorbing boundary condition (ABC) to remove unwanted reflection from the computational boundary [8, 9].

Sudiarta and Geldart [5] showed that the Fourier transformation procedure and the ABC are not needed when the TDSE in imaginary time is used to compute eigen energies and eigen functions. The TDSE in real time is transformed into the TDSE in imaginary time (or a diffusion equation) by replacing the real time variable with an imaginary time variable. The diffusion equation is then used to simulate a decaying of a wave function to the ground state wave function. The ground state energy and the ground state wave function can be computed directly in single simulation run. Beside that another advantage of the FDTD method in imaginary time is that it uses only real values of the wave function therefore it requires less computer memory.

The FDTD method in imaginary time (called simply the FDTD method for the rest of this paper) has been applied for various quantum systems. Sudiarta and Geldart [10] have used the FDTD method for an electron in a quantum dot and in a magnetic field. Dumitru et al. [11] have studied quarkornium states in QCD plasma. Strickland and Yager-Elorriaga [12] have used the FDTD method with parallel algorithm for a cumputer cluster to simulate a large quantum system. Margotta et al [13] have used the FDTD method to study quarkonium binding energies using a realistic complex-valued potential. The FDTD method has also been used by Alford and Strickland [14] to study quarkonium states in an external magnetic field. Sudiarta and Geldart [15] have applied the FDTD method to compute single particle density matrices. Sudiarta and Angraini (2016) have also applied the FDTD method with the super symmetry technique to obtain ground and excited states of a particle in one dimensional potentials.

For a quite small spatial grid spacing, the FDTD method has been found to give accurate results for one to three dimensional quantum systems [5]. For larger systems the FDTD method with small spacing is not possible due to limited computer memory. One can overcome this problem by using a cluster of computers and doing computation in parallel which can be easily done for the FDTD method as presented by Strickland and Yager-Elorriaga [12].

Another method to overcome this limitation is to increase its accuracy by using a nonstandard finite difference (NSFD) scheme that was found to be successfully applied for various differential equations [20, 17, 19]. The NSFD scheme for the Schrödinger equation was applied by Mickens and Ramadhani [20]. It is shown by Mickens and Ramadhani [20] and Chen et al. [21] that the nonstandard finite difference scheme is performing better than the standard method and it is also found to be exact for a constant potential. Application of the FDTD method with the NSFD scheme has not been done previously. In this paper, the NSFD scheme is used to modify the standard FDTD method such that the nonstandard FDTD (or NSFDTD) method can be used for large spatial grid spacing.

The NSFD scheme is generally applied for a system with one spatial dimension. Extension to two or more spatial dimensions is currently in progress and it will be published in our future publication. One possible extension is given by Cole [22, 23] for electromagnetic waves. In this paper only the theory and numerical method of the NSFDTD method for one spatial dimension are presented.

This paper is organized as follows: the theory of the NSFDTD method is introduced in the next section, numerical results of energies and wavefunctions for three potential wells are given in section three and the conclusion is given in section four.

2 Theory

The time-dependent Schrödinger equation (TDSE) in imaginary time for a system of one particle in a one dimensional potensial $V(x)$ using atomic units ($\hbar=m=1$) is given by [5]

[TABLE]

The corresponding time-independent Schrödinger equation (TISE), which is an eigen value equation, is given by

[TABLE]

where $W(x)=2[E-V(x)]$. It is noted that $E$ is the eigen energy.

A solution of the TISE can be obtained iteratively by using Eq. (1). Using an arbitrary initial wavefunction $\psi(x,0)$ and a suitable discretization scheme, a simulation with the discretized Eq. (1) can be done to evolve the initial wavefunction. After a long enough simulation time, the wavefunction converges to the ground state of the system. Excited states of the system can be obtained in turn by evolving the initial wavefunction with the component of lower energy wavefunctions removed from the wavefunction [5]. For computational purpose, the computational domain is limited and the boundary of computation domain is terminated by a Dirichlet boundary condition, $\psi(boundary)=0$. This boundary condition is equivalent to the infinite potential wall. After the wavefunction is obtained, the corresponding energy is computed with

[TABLE]

In the FDTD method [5], the standard discritization for the second derivative of Eq. (1) is the central finite difference scheme [24] in the form of

[TABLE]

where we have used a notation $\psi(i)\equiv\psi(i\Delta x)$ and $\Delta x$ is the spatial grid spacing.

For certain applications, the NSFD scheme developed by Mickens [16, 18] has been shown to give exact numerical results. A nonstandard scheme for the Schrödinger equation given by Mickens and Ramadhani [20] can be found by using a Taylor series for finite difference of the second derivative as the following:

[TABLE]

where $\psi^{\prime\prime}=\partial^{2}\psi/\partial x^{2}$, $\psi^{(4)}=\partial^{4}\psi/\partial x^{4}$, and $\psi^{(6)}=\partial^{6}\psi/\partial x^{6}$. Equation (2) with a constant potential, $W(x)=W$, can be used to simplify Eq. (5). It can be shown that

[TABLE]

Equation (6) can be simplified further using Maclaurin series of $2\sin^{2}(\theta/2)$, we obtain

[TABLE]

Therefore, the second derivative can be accurately approximated by

[TABLE]

where the denominator function $g(i)\equiv g(i\Delta x)=g(x)$ is given by

[TABLE]

. $g(x)$ can be rewritten in a simpler form using a Taylor series as

[TABLE]

It is noted in Eq. (10) that the NSFD scheme reduces to the standard scheme when the value $W(x)$ is closed to zero or the spatial interval $\Delta x$ is small. The NSFD scheme of Eq. (8) is found to give exact results for Eq. (2) when $W(x)$ has no spatial dependence or a constant as indicated by Eq. (7). This is in agreement with Mickens and Ramadani [20] and Chen et al. [21]. When $W(x)$ is not a constant, the nonstandard scheme has accuracy of the order $(\Delta x)^{2}$ and known to perform better than other methods such as the standard finite difference method and the Numerov method for large values of $\Delta x$ [21].

In order to apply the NSFDTD scheme, initial values of eigen energies $E$ are required as input parameters. The first estimates of energies can be obtained by using the standard FDTD method or by using $E=0$. The values of eigen energies are then improved by using the NSFDTD method. Further improvement can be done by reapplying the NSFDTD method. This proses can be iterated until converged energies are obtained.

3 Numerical Methods

An explicit iterative procedure for the NSFDTD method is following Sudiarta and Geldart [5]. The temporal derivative in Eq. (1) is approximated by the forward difference equation and the spatial second derivative is approximated using Eq. (8). For better numerical stability and accuracy, $\psi(x,t)$ in Eq. (1) is equal to an average value $(\psi^{n+1}(i)+\psi^{n}(i))/2$. After some manipulations, it is found that an explicit equation is given by

[TABLE]

where $\alpha(i)=[1-\Delta tV(i)]/[1+\Delta tV(i)]$ and $\beta(i)=1/[1+\Delta tV(i)]$.

Using initial random values for an initial wavefunction $\psi(i)$, we then iterate with Eq. (11) and after enough number of iterations, the wavefunction converges to the groundstate wavefunction. Excited states are obtained by redoing the iteration with lower states removed from the wavefunction (see Sudiarta and Geldart [5] for detail).

The energy is computed numerically by

[TABLE]

.

As previously mentioned that in order to use the NSFDTD method, initial eigen energies are needed and in this paper we obtain them by performing iteration with $g(i)=(\Delta x)^{2}$ or setting $W(x)=0$ in Eq. (9). The resulting energies and wavefunctions for this case are the same as the numerical results found using the standard FDTD method. Alternatively we can also use $E_{n}=0$ as the initial eigen energies and it is found that same final results of the NSFDTD method are obtained regardless the values of the initial eigen energies. In this paper, only numerical results using $W(x)=0$ are given.

4 Results and Discussion

4.1 Infinite Square Well Potential

To demonstrate significant improvements when using the NSFDTD in particular for a constant potential $V(x)=V_{0}$ or simply $V(x)=0$, we first consider a free particle in an infinite square well given by $V(x)=0$ for $0<x<1$ and $V(x)=\infty$ for $x<0$ and $x>1$. The eigen energies are known to be $E_{n}=n^{2}\pi^{2}/2$ and the eigen states are given by $\psi_{n}(x)=\sqrt{2}\sin(n\pi x)$. The computational domain is terminated using $\psi(x)=0$ at $x=0$ and $x=1$ which corresponds to an infinite potential wall.

The numerical parameters used for computation are $\Delta x=0.1$ and $\Delta t=(\Delta x)^{2}/10$. The numerical results for eigen energies and state functions are shown in Table 1 and Fig. 1.

It is noted that after the first iteration the NSFDTD results have shown better accuracy than the standard FDTD method. Further reapplication of the NSFDTD method produces even better eigen energies. After about 10 iteration, the final eigen energies converges to the exact values. These results confirm that the NSFDTD method is exact method when the potential is constant. Furthermore, the resulting wavefunctions are in good agreement with the analytical wavefunctions and the same as the standard FDTD results. This is a consequence of the fact that the values of $g(i)$ are constants. Discritization of the TISE using the standard and nonstandard schemes yields difference equations that differ only in an energy scaling factor.

4.2 Harmonic Oscillator Potential

Second application of the NSFDTD method is for a particle in an harmonic oscillator potential $V(x)=\frac{1}{2}\omega x^{2}$. The eigen energies are $E_{n}=(n+\frac{1}{2})\omega$. The wavefunctions are products of $\exp(-\omega x^{2}/2)$ and Hermite polynomials (details can be found in [25, 26]). Numerical results of eigen energies for an angular frequency $\omega=1$, grid spacing $\Delta x=0.8$ and $\Delta x=0.5$ with computational length of 16 and 10 are given in Table 2. Numerical wave functions for the first five states are shown in Fig. 2. To compare accuracy of the FDTD and the NSFDTD results, absolute errors of numerical wavefunctions for states ($n=3$ and $n=4$) are shown in Fig. 3.

In this application, the potential $V(x)$ and consequently the $W(x)$ are dependent on $x$. It is clearly shown, especially for excited states, the first iteration of the NSFDTD method is also found to improve significantly the standard FDTD method. Although the NSFDTD method is not exact in this case, its accuracy is still higher than the standard FDTD results. It is obviously shown that the numerical results become more accurate when the spacing is reduced. Because the wavefunction of the harmonic oscillator has exponential dependence, in order to get accurate results using the NSFDTD method we need to verify that the wavefunction has decayed to approximately zero near the computational boundary. This is to ensure that the boundary does not affect wavefunctions.

Comparison of numerical and analytical wavefunctions are shown in Fig. 2 and the sum of squared errors (SSE) of the wavefunctions are given in Table 3. The numerical wavefunctions are shown to be in agreement with the analytical wavefunctions. It is noted that the SSE of wavefunctions is decreasing as we increase the number of iteration. It is also shown that the wavefunctions are very small near the computational boundaries. Numerical wavefunctions obtained by the NSFDTD method is more accurate than by the FDTD method which is clearly shown by smaller absolute errors in Fig. 3.

4.3 Morse Potential

For application in atomic vibrations, the Morse potential well is usually chosen as the model potential. The Morse potential as an example for this paper is in the form of

[TABLE]

. The analytical eigen energies are given by $E_{n}=\omega(n+\frac{1}{2})[1-\omega(n+\frac{1}{2})/4D]$ with $\omega^{2}=2D$.

For numerical computations, the NSFDTD parameters are $x_{0}=2$, $\Delta x=0.2$, the computational domain is $[0,8]$ and the Morse parameter $D=50$ is used for simplicity. The resulting numerical eigenvalues are given in Table 4. In this case, it is also found that the NSFDTD method is more accurate that the standard FDTD method.

5 Conclusion

A numerical method known as the nonstandard (NS) FDTD method has been presented and validated using three potential wells. The numerical results of eigen energies and wavefunctions are compared with the standard FDTD and the analytical results. It has been shown that the NSFDTD method is an exact method for the case of a constant potential. The NSFDTD method has been also shown to perform better than the standard FDTD method for all cases.

Acknowledgement

This work is partially supported by Directorate Research and Community Service - Directorate General of Strengthening for Research and Development - Ministry of Research, Technology, and Higher Education Republic of Indonesia.