The Finite Difference Time Domain (FDTD) Method to Determine Energies and Wave Functions of Two-Electron Quantum Dot

TL;DR

This paper demonstrates the application of the finite difference time domain (FDTD) method to accurately compute energies and wave functions of two electrons in a three-dimensional quantum dot modeled by a harmonic potential, using imaginary time evolution.

Contribution

The paper introduces the FDTD method for solving the time-dependent Schrödinger equation in imaginary time to find energies and wave functions of two-electron quantum dots, including excited states with orthogonality constraints.

Findings

FDTD provides accurate energies and wave functions for quantum dots.

The method successfully computes excited states with orthogonality constraints.

Results agree well with other numerical methods.

Abstract

The finite difference time domain (FDTD) method has been successfully applied to obtain energies and wave functions for two electrons in a quantum dot modeled by a three dimensional harmonic potential. The FDTD method uses the time-dependent Schr\"odinger equation (TDSE) in imaginary time. The TDSE is numerically solved with an initial random wave function and after enough simulation time, the wave function converges to the ground state wave function. The excited states are determined by using the same procedure for the ground state with additional constraints that the wave function must be orthogonal with all lower energy wave functions. The numerical results for energies and wave functions for different parameters of confinement potentials are given and compared with published results using other numerical methods. It is shown that the FDTD method gives accurate energies and wave functions.

Tables

Table 1. Table 1: Numerical ground and first excited state energies with its errors for two electrons in a box (side length = 1) computed by the FDTD method using various number of grids N 𝑁 N or grid spacing Δ ​ x Δ 𝑥 \Delta x . The exact ground dan first excited state energies are E 0 = 29.608813 subscript 𝐸 0 29.608813 E_{0}=29.608813 and E 1 = 44.413220 subscript 𝐸 1 44.413220 E_{1}=44.413220 . Numerical extrapolated energies are computed by extrapolating function Eq. ( 17 .

N	$\Delta x$	$E_{FDTD,0}$	$E_{FDTD,1}$	$E_0-E_{FDTD,0}$	$E_1-E_{FDTD,1}$
8	0.125000	29.230260	43.103716	0.378553	1.309503
10	0.100000	29.366090	43.570042	0.242723	0.843177
12	0.083333	29.440087	43.825748	0.168726	0.587471
14	0.071429	29.484776	43.980749	0.124037	0.432470
16	0.062500	29.513810	44.081681	0.095003	0.331538
Extrapolation (see text)		29.608196	44.406182	0.000618	0.007038

Table 2. Table 2: Numerical ground state energies computed using the FDTD method with extrapolation (see Fig.  3 ).

$\omega$	FDTD	ref [12]
0.25	1.092382	1.089262
1	3.734785	3.730120
4	13.536860	13.523214

Equations

$$i\hslash \frac{\partial }{\partial t}\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,t)=\hat{H}\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,t)=\left[-\sum _{j=1}^2\frac{{\hslash }^2}{2m}{{\nabla }_j}^2+V({\mathbf{r}}_1,{\mathbf{r}}_2)\right]\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,t)$$

$$V({\mathbf{r}}_1,{\mathbf{r}}_2)=V_{dot}({\mathbf{r}}_1)+V_{dot}({\mathbf{r}}_2)+\frac{e^2}{4\pi {ϵ}_0|{\mathbf{r}}_1-{\mathbf{r}}_2|}$$

$$\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,t)=\sum _{n=0}^{\infty }{c}_n{\varphi }_n({\mathbf{r}}_1,{\mathbf{r}}_2)\exp (-i{E}_{n}t/\hslash )$$

$${\varphi }_n({\mathbf{r}}_1,{\mathbf{r}}_2)=-{\varphi }_n({\mathbf{r}}_2,{\mathbf{r}}_1)$$

$$\frac{\partial }{\partial \tau }\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )=\frac{1}{2}\sum _{j=1}^2{{\nabla }_j}^2\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )-V({\mathbf{r}}_1,{\mathbf{r}}_2)\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )$$

$$\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )=\sum _{n=0}^{\infty }{c}_n{\varphi }_n({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )\exp (-E_n\tau )$$

$$\underset{\tau \to \infty }{\lim }\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )\approx [c_0\exp (-E_0\tau )]{\varphi }_0({\mathbf{r}}_1,{\mathbf{r}}_2)$$

$$E=\frac{<\psi |\hat{H}|\psi >}{<\psi |\psi >}=\frac{\int {\psi }^{\ast }\hat{H}\psi d^3{r}_1{d}^3{r}_2}{\int {\left|\psi \right|}^2{d}^3{r}_1{d}^3{r}_2}$$

$${\psi }^{n+1}(i_1,j_1,k_1,i_2,j_2,k_2)=\alpha {\psi }^n(i_1,j_1,k_1,i_2,j_2,k_2)+\beta {\phi }^n(i_1,j_1,k_1,i_2,j_2,k_2)$$

$$\varphi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2},k_{2})=\left[\begin{array}[]{l}\psi^{n}(i_{1}+1,j_{1},k_{1},i_{2},j_{2},k_{2})+\psi^{n}(i_{1}-1,j_{1},k_{1},i_{2},j_{2},k_{2})\\ +\psi^{n}(i_{1},j_{1}+1,k_{1},i_{2},j_{2},k_{2})+\psi^{n}(i_{1},j_{1}-1,k_{1},i_{2},j_{2},k_{2})\\ +\psi^{n}(i_{1},j_{1},k_{1}+1,i_{2},j_{2},k_{2})+\psi^{n}(i_{1},j_{1},k_{1}-1,i_{2},j_{2},k_{2})\\ +\psi^{n}(i_{1},j_{1},k_{1},i_{2}+1,j_{2},k_{2})+\psi^{n}(i_{1},j_{1},k_{1},i_{2}-1,j_{2},k_{2})\\ +\psi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2}+1,k_{2})+\psi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2}-1,k_{2})\\ +\psi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2},k_{2}+1)+\psi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2},k_{2}-1)\\ -12\psi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2},k_{2})\end{array}\right]$$

$$\alpha (i_1,j_1,k_1,i_2,j_2,k_2)=\frac{[1-\frac{\Delta \tau }{2}V(i_1,j_1,k_1,i_2,j_2,k_2)]}{[1+\frac{\Delta t}{2}V(i_1,j_1,k_1,i_2,j_2,k_2)]}$$

$$\beta i_1,j_1,k_1,i_2,j_2,k_2=\frac{\Delta \tau }{2(\Delta x{)}^2[1+\frac{\Delta t}{2}V(i_1,j_1,k_1,i_2,j_2,k_2)]}$$

$$\Delta \tau \le (\Delta x{)}^2/6$$

$${\psi }_{AS}({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )=\frac{1}{2}[\psi ({\mathbf{r}}_1,{\mathbf{r}}_2,\tau )-\psi ({\mathbf{r}}_2,{\mathbf{r}}_1,\tau )]$$

$${\psi }_{AS}^n(i_1,j_1,k_1,i_2,j_2,k_2)=\frac{1}{2}[{\psi }^n(i_1,j_1,k_1,i_2,j_2,k_2)-{\psi }^n(i_2,j_2,k_2,i_1,j_1,k_1)]$$

$$E=\frac{1}{\sum \psi (i_1,\dots ,k_2{)}^2}\sum \left\{V(i_1,\dots ,k_2)\psi {(i_1,\dots ,k_2)}^2-\psi (i_1,\dots ,k_2)\phi (i_1,\dots ,k_2)/2\Delta x^2\right\}$$

$$E(\Delta x)=E-m(\Delta x{)}^2$$

Full text

1 INTRODUCTION

The finite difference time domain (FDTD) method is a general numerical method that can be applied to solve many physical problems such as electromagnetic wave simulations [1], acoustic wave simulations [2] and quantum simulations [3, 4, 5, 6]. The FDTD method has many advantages. Beside its good accuracy, it is also easy to code and very flexible for systems of any geometrical shapes.

To solve the time-independence Schrödinger equation, the FDTD method can be devided into two methods: real-time FDTD (R-FDTD) method [3, 5] and imaginary-time FDTD (I-FDTD) method [6, 7, 8, 9]. The R-FDTD method uses evolution of a wavefunction by the discretized time-dependence Schrödinger equation and Fourier transformation procedure to obtain eigen energies and wavefunctions. The I-FDTD method uses a diffusion equation to evolve an initial wavefunction and after a large interval of time the wavefunction approaches the ground state wavefunction. Advantanges of using I-FDTD method are (1) it only uses real variables compared to complex variables in the R-FDTD method for non magnetic case, (2) it does not require Fourier transformation and an absorbing boundary condition to truncate the computational domain.

In this paper, we apply the I-FDTD method to compute energies and wavefunctions of two electrons in a quantum dot. In the two-electrons system we need to descritize a six dimensional differential equation that requires very large number of grids. One approach to overcome this problem is to use an approximate model such as Hartree-Fock (HF) model [3]. The HF model does not contain correlation term which is important to simulate correlated system. Therefore solving full numerical Schrödinger equation is desired. In this paper, in order to get accurate results, an extrapolation procedure as given in [10] is done using FDTD results with various grid spacing parameters.

2 THEORY

In this paper, we briefly explain the FDTD method. A detail explanation of the FDTD method can be found in [6]. To apply the FDTD method, we begin with the time-dependent Schrödinger equation (TDSE) for two electrons given by

[TABLE]

where $m$ is the effective mass of electron and $V(\mathbf{r}_{1},\mathbf{r}_{2})$ is the potential given by

[TABLE]

with an harmonic confining potential of a quantum dot $V_{dot}(\mathbf{r})=\frac{1}{2}m\left(\omega_{x}^{2}x^{2}+\omega_{y}^{2}y^{2}+\omega_{z}^{2}z^{2}\right)$.

The solution of Eq. (1) can be expressed in the form of expansion of eigen functions as

[TABLE]

where $c_{n}$ is the expansion coefficient, the eigen functions $\phi_{n}(\mathbf{r}_{1},\mathbf{r}_{2})$ and the energies $E_{n}$ are obtained by solving the time-independent Schrödinger equations (TISE), $\hat{H}\phi=E\phi$.

Beside the Schrödinger equation, the wavefunction $\phi_{n}(\mathbf{r}_{1},\mathbf{r}_{2})$ for electrons having the same spin must also satisfy the Pauli exclusion principle where the wave function must be anti-symmetric under exchange of two electrons [11], i.e.

[TABLE]

For simplicity in programming, in the remaining part of this paper we use atomic units with $\hbar=m=1$ and $e^{2}/4\pi\epsilon_{0}=1$. The first step in solving the TISE with the TDSE (Eq. (1)) using the FDTD method is to transform Eq. (1) into a diffusion equation by substituting $\tau=it$ in Eq. (1). The resulting diffusion equation is

[TABLE]

and the wavefunction is

[TABLE]

. It can be noted from Eq. (6) that for a large imaginary time $\tau$, because of the factor $\exp(-E_{n}\tau)$, wave functions $\phi_{n}({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}})$ with high energies diminish and the lowest energy state (or ground state) remains. Therefore after a large simulation time, the final wavefunction is

[TABLE]

where $[{{c}_{0}\exp(-{{E}_{0}}\tau)}]$ is just a constant factor which can be removed by normalizing the final wavefunction. Using an initial function $\psi({{\mathbf{r}}_{1}},{{\mathbf{r}}_{2}},\tau=0)$, we can obtain the ground state wave function by a simulation using Eq. (5) for long enough interval of $\tau$ provided that the initial function contains the ground state or ${{c}_{0}}\neq 0$. Excited states can be obtained using same procedure provided that all lower energy wavefunctions have been removed from the initial function.

After obtaining a wavefunction $\psi$, its energy can be computed by

[TABLE]

3 NUMERICAL METHOD

For numerical computation, we use a notation ${{\psi}^{n}}({{i}_{1}},{{j}_{1}},{{k}_{1}},{{i}_{2}},{{j}_{2}},{{k}_{2}})=\psi({{i}_{1}}\Delta x,{{j}_{1}}\Delta x,{{k}_{1}}\Delta x,{{i}_{2}}\Delta x,{{j}_{2}}\Delta x,{{k}_{2}}\Delta x,n\Delta\tau)$ (where $\Delta\tau$ and$\Delta{{x}}$ are the temporal and spatial spacing. Following derivation in [6], it can be obtained that an explicit numerical iterative scheme for simulating evolution of a wavefunction is given by

[TABLE]

where $\varphi^{n}(i_{1},j_{1},k_{1},i_{2},j_{2},k_{2})$ is

[TABLE]

and the coefficients $\alpha$ and $\beta$ are

[TABLE]

[TABLE]

In order to have a stable iteration, the time step $\Delta\tau$ must satisfy a stability condition given by

[TABLE]

.

In this paper to initiate the iteration, a random wavefunction is used as the initial wavefunction and the computational boundary is set to be zero or $\psi(\mathbf{r}_{1},\mathbf{r}_{2})_{boundary}=0$. The wavefunction for same spin electrons must also satisfy the Pauli exclusion principle which can be accomplished by updating the wavefunction at every step using

[TABLE]

or in a discretized form given by

[TABLE]

The energy is computed numerically by an approximation of Eq.  (8) given by

[TABLE]

4 RESULTS AND DISCUSSIONS

4.1 Non Interacting Electrons In a Box

To validate the FDTD method, we perform FDTD simulations for two non-interacting electrons in a cubic box with side length of one. The potential of the box is given by $V_{box}(\mathbf{r})=0$ inside the box ($|x|<1/2$, $|y|<1/2$, $|z|<1/2$) and $V_{box}(\mathbf{r}_{1},\mathbf{r}_{2})=\infty$ outside the box. Numerical results for ground state and first excited state energies using various grid parameters are given in Table 1. It is noted that the numerical FDTD results approach exact energies as the grid spacing becomes smaller.

The grid spacing is equal to the length of computation domain $L$ divided by number of grids $N$, that is $\Delta x=L/N$. The FDTD computation for two electrons requires $(N+1)^{6}$ grid points. This needs more than $40\times(N+1)^{6}$ bytes of computer memory when using double precision computation. Using a regular PC with four gigabites of memory can do computations with maximum number of grids per dimension $N=17$.

To get accurate results we perform extrapolation using FDTD results for various spatial intervals $\Delta x$. Because the FDTD method uses central finite difference (CFD) scheme for the second derivative in the Schrödinger equation, the numerical errors of the CFD scheme are proportional to $(\Delta x)^{2}$. Numerical errors for ground state and excited state energies are plotted in Fig. 1. The numerical errors are perfectly matched by a linear curve fitting. Therefore a fitting function used in this paper in the form of

[TABLE]

is suitable for extrapolation. Other extrapolating functions such as Padé approximating function and the Richardson method may also be used [10].

Ground state and excited state energies obtained by extrapolation are given in Table 1. It is shown clearly that a significant improvement is resulted when using this extrapolation method.

As an example of comparison between numerical and analytical wavefunctions are shown in Fig. 2 for values at the $x$ axis, that is the values of $\psi(x,0,0,0,0,0)$. The numerical wavefunction is in a good agreement with the analytical wavefunction.

4.2 Two Electrons in an Isotropic Quantum Dot

To demonstrate further application of the FDTD method, we consider two interacting electrons with different spins in an isotropic quantum dot ($\omega_{x}=\omega_{y}=\omega_{z}$ with three harmonic confining potentials ($\omega=0.25,1,4$). Numerical results for ground state energies for four grid spacings are given in Fig.  3. Similar to previous section, linear variation of energies as a function of $(\Delta x)^{2}$ are also shown for this case. After extrapolation procedure (see Fig.  3), numerical ground state energies for the three confining potentials can be found in Table 2. The FDTD results are in a good agreement with results of the discrete variable representation (DVR) method given in [12]. A variation of numerical wavefunction of two electrons in quantum dot along x-axis is given in Fig. 4. It is expected due to harmonic confinement and singlet state, the maximum wavefunction is located in the middle of the quantum dot.

5 CONCLUSION

A numerical method known as the FDTD method to obtain eigen energies and wavefunctions for two electrons in a confining potential has been presented. The FDTD method for two electrons requires large amount of computer memory, therefore limited number of grids with large grid spacings are used. One method to overcome this limitation is by performing FDTD simulations with various grid spacings and then extrapolating the FDTD results. It has been shown that errors of the FDTD results are proportional to the square of the grid spacing. The FDTD method with extrapolation gives accurate eigen energies and wavefunctions that have been validated by comparing with results from other analytical and numerical methods.

6 ACKNOWLEDGMENTS

This work is partially supported by the University of Mataram under PNBP research scheme.