Chinese Journal of Chemical Physics  2020, Vol. 33 Issue (2): 258-262

#### The article information

Yin Huang, Hai-lin Zhao, Syed Kazim Usman, Ganiyu Ayodele Ajibade, Zhi-gang Sun

High-Accurate Transparent Boundary Conditions for Time-Dependent Quantum Wave Packet Method

Chinese Journal of Chemical Physics, 2020, 33(2): 258-262

http://dx.doi.org/10.1063/1674-0068/cjcp1909167

### Article history

Accepted on: January 3, 2020
High-Accurate Transparent Boundary Conditions for Time-Dependent Quantum Wave Packet Method
Yin Huanga , Hai-lin Zhaoa , Syed Kazim Usmana,b , Ganiyu Ayodele Ajibadea,b , Zhi-gang Suna
Dated: Received on September 23, 2019; Accepted on January 3, 2020
a. State Key Laboratory of Molecular Reaction Dynamics and Center for Theoretical Computational Chemistry, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, Dalian 116023, China;
b. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Complex absorbing potential is usually required in a time-dependent wave packet method to accomplish the calculation in a truncated region. Usually it works effectively but becomes inefficient when the wave function involves translational energy of broad range, particularly involving ultra-low energy. In this work, a new transparent boundary condition (TBC) is proposed for the time-dependent wave packet method. It in principle is of spectral accuracy when typical discrete variable representations are applied. The prominent merit of the new TBC is that its accuracy is insensitive to the translational energy distribution of the wave function, in contrast with the complex absorbing potential. Application of the new TBC is given to one-dimensional particle wave packet scatterings from a barrier with a potential well, which supports resonances states.
Key words: Quantum wave packet    Absorbing potential    Transparent boundary conditions

Solving time-dependent Schrödinger equation is of vital importance in quantum molecular dynamics field nowadays, due to its advantages on the numerical scaling, the conceptual simplicity and the straightforward pictorial way over the time-independent method. As we know, in a wave packet calculation, the area of computation must be limited to appropriate size. In order to avoid affection of the undesirable spurious reflections at the boundaries in the long time period, in a modern time-dependent wave packet calculation, usually a complex absorbing potential (CAP) of the form $V(x)$ = $-if(x)$ or an equivalent damping function is introduced near the edge of the grid to attenuate the wave function gradually [1-3]. Due to its simplicity and local nature, the CAP gains general applications and usually is efficient.

However, when the distribution of the translational energy in the wave packet is of broad range, especially involving ultra-slow translational energy, the CAP becomes very inefficient. Since it must be strong enough to damp the part moving rapidly, at the same time, it must go up gently enough to avoid the reflection of the part moving slowly. Quite often, in order to keep the reflection and transmission less than 0.1%, the length of the grid range devoted to the CAP has to be several times of the largest de Broglie wavelength of the particle. That is, the CAP with length of tens even several hundred bohrs is required in a scattering processes involving H atom with translational energy less than 1.0 meV. The total computational area then becomes huge and difficult to handle with for describing a realistic reactive scattering process. Other methods for dealing with the open boundary conditions, such as perfect matched layers (PML), exterior complex scaling (ECS) and smooth exterior scaling (SES), suffer similar numerical difficulties to the CAP, even in principle they are more accurate [4].

Recent interest in the simulation of tunnelling electronic devices has prompted the development of a number of open-system boundary conditions [5, 6]. Exact boundary condition has been derived by the Laplace transform in a convolution form for the linear equation at the asymptotic region where the potential is a constant. Although this exact boundary condition in principle is global and continuous, its temporally semi-discretised form via the $Z$ transform and fully discretised form in both time and space have been treated [7, 8]. This boundary condition is insensitive to the energy distribution of the initial wave packet, only if the numerical parameters are capable of giving converged results in the interior region. As an alternative, the Fourier transform of the time-dependent Schrödinger equation may be adopted, which leads to boundary conditions being local both in space and time. This method results in designing an effective approximation for a square root $\sqrt{2\mu E -V(x)}$ at the boundary point $x_M$, which anyway can be easily transformed into a differential equation [9-13]; here $E$ denotes the energy of \newpage the wave packet, $\mu$ is the reduced mass of the system and $V$ is the potential energy. The accuracy of this kind of method depends strongly on the approximation of the square root $\sqrt{2\mu E -V(x)}$, i.e., depends on the energy distribution of the wave packet. The above two kinds of methods, named as transparent boundary conditions (TBC), often work with low-order finite difference method thus of low order accuracy, since their forms with high-order finite difference method are complicated [14].

In a modern quantum wave packet calculation, usually the discrete variable representation method is taken as the basis functions to realise the transformation between coordinate (or grid) and basis representation (DVR), which is a pseudo-spectral method. Thus the TBC, which is designed with low order finite different method, is not easy to apply in a wave packet calculation.

In this work, to crack the open boundary conditions difficulty, a new TBC, which is adapted to the DVR method, is proposed. It is expected to be particularly useful in a calculation where energy distribution is broad and involves ultra-low translational energy. At the same time, its implementation is simple and straightforward, and of high-order accuracy. We take the one-dimensional motion of a particle to illustrate the proposed method.

The motion of a particle with the mass $\mu$ in the interval $x $$\in$$ [-x_M, x_M]$ may be described by the one-dimensional Schrödinger equation

 $$$i\frac{\partial }{\partial t} \Psi(x, t) = -\frac{1}{2\mu}\frac{\partial^2}{\partial x^2} \Psi(x, t) +V(x) \Psi(x, t)$$$ (1)

where atomic unit (a.u.) is applied. At time $t$ = 0, the initial wave packet $\Psi$($x_0$, $t$ = 0) moves from right to left, starting from central position $x_0$ with central energy $E_0$, which is defined as

 $$$\Psi(x,0) = \exp\left[\frac{-(x-x_0)^2}{2\sigma_0}\right]\exp(ixk_0)$$$ (2)

here $\sigma_0$ defines the width and thus the energy distribution of the initial wave packet, and $k_0$ = $\sqrt{2\mu E_0}$. The potential energy is defined by

 $$$V(x) = V_0\left[\frac{ \rm{e}^{\beta x}}{(1+ \rm{e}^{\beta x})^2}-\alpha\frac{ \rm{e}^{5\beta x}}{(1+ \rm{e}^{5\beta x})^2}\right]$$$ (3)

which is presented in FIG. 1(B) with $V_0$ = 2.0 eV, $\alpha$ = 1.0 a.u. and $\beta$ = 0.3 a.u.

 FIG. 1 (A) The wave function (navy lines) at a series of time and (B) the potential energy, traveling from right to left. The red lines in the boundary regions are calculated by the CG-TBC, which are used to update the wave function around boundaries. The grey area in panel (B) indicates the CAP.

The numerical solution of Eq.(1) can only be performed on a finite interval, which we choose as $x $$\in$$ [-x_M, x_M]$, discretised as $-x_M, -x_{M-1} \dots, x_i, \dots$, $x_{M-1}$, $x_M$ with equal spacing $\Delta x$ when sincDVR is applied. One of the appropriate boundary conditions to chose $\Psi(\pm x_M)$ = 0 for the case that the particle is far from the boundaries. However, reflections occur for a long time propagation when the particle could approach to the boundaries.

Our new TBC to dealing the boundary conditions could be derived as the following. The wave packet in the boundary region $[-x_M, -x_1]$ (see FIG. 1), where usually the potential energy becomes a constant, may be expanded as a series of plane waves of the form around $-x_M$ (boundary conditions around $x_M$ can be similarly derived)

 $$$\bar{\Psi}(x_i, t) = \sum\limits_{n = 0}^{N}a_n\exp(-i\omega_n t)\exp(-i k_n x_i)$$$ (4)

where $\omega_n$ and $k_n$ are the corresponding energy and momentum, with $k_n$ = $\pm\sqrt{[2\mu(\omega_n-V)]}$. The sign of $k_n$ determines the traveling direction of the waves. The wave function $\Psi(x_i, t+\Delta t)$ then usually is written as

 $$$\bar{\Psi}(x_i, t+\Delta t) = \sum\limits_{n = 0}^{N}a_n\exp[-i\omega_n( t +\Delta t)]\exp(-i k_n x_i)$$$ (5)

which leads to spurious reflections at the grid boundaries when the wave function closes. Physically, there is only left outgoing wave function in the boundary region $[-x_M$, $-x_1]$, we may effectively write the wave function $\Psi$($x_i$, $t$+$\Delta t$) in this way

 $$$\bar{\Psi}(x_i, t+\Delta t) = \sum\limits_{n = 0}^{N}a_n\exp(-i\omega_n t)\exp[-i k_n (x_i+v_n\Delta t)]$$$ (6)

where $v_n$ = $k_n/\mu$ is the velocity of the particle. The problem is that $a_n$ is unknown and is not easy to calculate accurately with the Fourier function, due to the discontinuity of the wave functions at the ends of the boundary regions.

However, we find that when the Chebyshev-Gauss (CG) quadrature, $Q_k$ = $\cos[(k$+0.5)$\pi$/($N$+1)], is adopted, we actually can extract $a_n$ accurately with $N$ no less than 200. For convenience, $\arccos Q_k$ can be simply chosen as the grid points $x_i$ in the boundary region $[-x_M, -x_1]$. It is noted that in order to resolve the component of the lowest translational energy, the length of boundary region should be comparable or longer than the largest wave length of the particle. Then the wave function in the boundary region $[-x_M, -x_1]$ may be written as

 $$$\Psi(x_i, t) = \sum\limits_{n = 0}^{N}b_n\exp(-i\omega_n t)\cos( k_n x_i)$$$ (7)

Then the wave function $\Psi$($x_i,$ $t$+$\Delta t)$ at the points around $-x_M$ may be written according to

 $$$\Psi(x_i, t) = \sum\limits_{n = 0}^{N}b_n\exp(-i\omega_n t)\cos[k_n (x_i+v_n\Delta t)]$$$ (8)

where $v_n$ = $k_n$/(2$\mu$) is the velocity of the particle. In a practical calculation, we do not update the wave function for all $x_i$(0$\le $$i$$ \le$$N$) in the boundary region, since around the ends of boundary region the wave function can not be accurately calculated due to the discontinuity. We only update about 3/5 grid points in the boundary region around $-x_M$ with the wave function given by Eq.(8). Although the wave function on the grid points quite close $-x_M$ is not that accurate, they are pure outgoing and would not go back to interior region. Thus we may expect that they have little influence on the calculation, which will be verified by the following numerical examples. We named this new TDC as the CG-TDC.

In FIG. 1(A), snapshots of the wave function at a series of times are displayed as navy lines, where the boundary conditions are satisfied by using Eq.(8). The initial wave function is constructed with $x_0$ as 75.0 a.u., $\sigma_0$ as 0.5 a.u. and $E_0$ as 1.1 eV. The calculated interval is [$-$130.0, 130.0] a.u. with total 1535 grid points. Second-order split operator with time step 10.0 a.u. was applied to propagate the initial wave function. The total propagation time is 8$\times$10$^5$ a.u. The red lines in the boundary region are the wave functions calculated by Eq.(8).

It is seen clearly that the wave functions at the points around the ends of red lines do not agree well with the navy lines. After updating the propagated wave function (navy lines) on the points around $-x_M$ by the values calculated using Eq.(8) (left part of the red lines, about 3/5), correct behaviours of the wave function during the propagation were guaranteed. In the calculation, the number of the grid points is taken as 240 and the wave function on the outside 140 grid points among them was used to replace the propagated wave function to satisfy the boundary condition.

In FIG. 2(A), total transmission probabilities calculated by the flux method (M1), the amplitude method (M2) and the correlation function method (M3) at position $x_f$ = 75.0 are displaced [15], together with the one calculated using the CAP with the flux method. In the CAP calculation, the same parameters as those in the CG-TBC calculation are adopted, except that the boundary conditions are satisfied with a CAP with length of 50.0 a.u. It is seen that the four sets of the results agree with each other excellently. The inset in panel (B) displays the enlarged plots of that in FIG. 2(B). The results calculated using the CG-TDC indeed are very accurate, and the results calculated using the M2 method are inferior to the other two methods.

 FIG. 2 (A) and (C) Total transmission probabilities as a function of collision energy, calculated using the flux method (M1), amplitude method (M2) and correlation function method (M3) with the CG-TBC and transmission probabilities calculated using the M3 method with the CAP. (B) Resonance wave functions at translational energy 0.36912 eV calculated using the CG-TBC and CAP. (D) The relative $L^2$ error in log scale as a function of propagation time, defined in Eq.(9).

The resonance wave functions at $E_ \rm{R}$ = 0.36912 eV obtained by the CG-TDC and the CAP are compared in FIG. 2(C), which is calculated by a time-energy Fourier transform. Except the parts in the boundary regions, the two wave functions agree with each other very well. In the subset, it is seen that starting from $x$ = 80.0 a.u., where the CAP arises, the wave function calculated with the CAP gradually damps to zero, whereas the wave function with the CG-TDC keeps its oscillations, due to its transparency nature. The relative $L^2$ error, which is defined as

 $$${\rm error} = \frac{ {\sum\limits_i}|\Psi(t,x_i)-\Psi^0(t,x_i)|^2}{ {\sum\limits_i}|\Psi^0(t,x_i)|^2}$$$ (9)

in logarithmic scale with base 10 is plotted against with time in FIG. 2(D). The standard $\Psi^0(t,x_i)$ was calculated with much larger grid interval [$-$520.0, 520.0] a.u. but with the same grid spacing. The error increases along with time slowly. At time 4$\times$10$^4$ a.u. the residual population only is about 0.03, the error still less than 5$\times10^{-5}$. This fact clearly demonstrates the high accuracy of the CG-TBC.

To strengthen the power of CG-TDC, in FIG. 3(A) total transmission probabilities from translational energy of 0.01 meV to 2.0 eV were presented using a potential, which has a resonance state at ultra-low energy, 0.6352 meV. The potential is given by Eq.(3) with $V_0$ = 0.01 eV, $\alpha$ = 2.0 a.u., and $\beta$ = 0.3 a.u. The initial wave function is set up with $x_0$ as 80.0 a.u., $\sigma_0$ as 0.35 a.u. and $E_0$ as 0.8 eV. The calculated interval is [$-$150.0, 150.0] a.u. with total 1535 grid points. There are 240 grid points (amount to interval 47.0 a.u.) in the CG-TBC, among them 160 points are used to update the wave function. The same calculation was carried out, with the CAP of similar length to 50.0 a.u., as the following form

 FIG. 3 Comparison of the total transmission probabilities (A), the real (B) and imaginary part (C) of the $S$-matrix, as a function of collision energy calculated using the CG-TBC and the CAP with boundary regions of length of 50 a.u. The enlarged plots in the subsets of panel (A) clearly demonstrate the rationality of the CAP strength.
 $$$V_{\rm CAP}(x) = -iC\frac{(x-x_{1/2})^2}{x_M-x_{1/2}}$$$ (10)

with $C$ = 0.002 a.u. which has been fully optimised for a minimal reflection of the wave function with small translational energy, at the same time, for a minimal transmission of the wave function with large translational energy, which finally would be reflected back into the interior region. Since the wavelength of the particle at 0.6 meV is about 29 a.u. with $\mu$ as 0.58 amu, we may expect that the CAP could not work well in this case. In order to resolve the lowest resonance peak, the initial wave packet was propagated for a long time as 10$^7$ a.u.

As seen in subsets in FIG. 3(A), there are considerable oscillations in the low and high energy part for the CAP results, which deviate from the CG-TDC results. The oscillations come from the quality of the CAP, which is not long enough. However, the CG-TDC results agree with the results using longer CAP (not displayed). Similar conclusion can be drawn with the real part of the $S$-matrix in FIG. 3(B) and the imaginary part of the $S$-matrix in FIG. 3(C). The $S$-matrix was calculated using the M3 method. The results in FIG. 3 demonstrate the numerical superiority of the CG-TDC over the CAP for limiting the computational interval.

To conclude this work, we would like to note that previous TDC in a similar spirit was accomplished by an effective approximation for satisfying the dispersion relation, $k$ = $\pm\sqrt{[2\mu(\omega-V)]}$, using the low order finite difference method thus of low order accuracy. Thus those TDC strongly depends on the energy distribution of the wave function. In the work by Hu and Zheng [16], although the wave function on the end grid point was calculated using FFT method, which is of spectral accuracy, the spatial derivative of the equation was calculated using second-order finite difference method. Our new proposed CG-TDC in principle is of global spectral accuracy, i.e., exhibits exponential convergence with respective to the number of grid points for both spatial derivative and boundary conditions and in the future work, one may find better ways to resolve $b_n$ to further improve the numerical efficiency of the TDC with similar ideas.

Acknowledgements: This work was supported by the National Natural Science Foundation of China (No.21733006, No.21825303 and No.21688102), the Strategic Priority Research Program of Chinese Academy of Sciences (No.XDB17010200).

Reference
 [1] R. Kosloff, and D. Kosloff, J. Comput. Phys. 63, 363(1986). DOI:10.1016/0021-9991(86)90199-3 [2] M. S. Child, Mol. Phys. 72, 89(1991). DOI:10.1080/00268979100100041 [3] D. Neuhauser, and M. Baer, J. Chem. Phys. 90, 4351(1989). DOI:10.1063/1.456646 [4] A. Nissen, H. O. Karlsson, and G. Kreiss, J. Chem. Phys. 133, 054306(2010). DOI:10.1063/1.3458888 [5] X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, and A. Schädle, Comm. Phys. Comp. 4, 729(2008). [6] G. Pang, L. Bian, and S. Q. Tang, Phys. Rev. E 86, 066709(2012). DOI:10.1103/PhysRevE.86.066709 [7] X. Antoine, and C. Besse, J. Comp. Phys. 188, 157(2003). DOI:10.1016/S0021-9991(03)00159-1 [8] A. Arnold, and M. Ehrhardt, J. Comp. Phys. 145, 611(1998). DOI:10.1006/jcph.1998.6043 [9] J. R. Hellums, and W. R. Frensley, Phys. Rev. B 49, 2904(1994). DOI:10.1103/PhysRevB.49.2904 [10] T. Shibata, Phys. Rev. B 43, 6760(1991). DOI:10.1103/PhysRevB.43.6760 [11] J. P. Kuska, Phys. Rev. B 46, 5000(1992). DOI:10.1103/PhysRevB.46.5000 [12] T. Fevens, and H. Jiang, SIAM J. Sci. Comput. 21, 255(1999). DOI:10.1137/S1064827594277053 [13] D. Ruprecht, A. Schädle, F. Schmidt, and L. Zschiedrich, SIAM J. Sci. Comput. 30, 2358(2008). DOI:10.1137/070692637 [14] J. W. Zhang, Z. Z. Sun, X. N. Wu, and D. S. Wang, Commun. Comput. Phys. 10, 742(2011). DOI:10.4208/cicp.280610.161110a [15] J. Q. Dai, and J. Z. H. Zhang, J. Phys. Chem. 100, 6898(1996). DOI:10.1021/jp9536662 [16] J. S. Hu, and C. X. Zheng, J. Comp. App. Math. 326, 116(2017). DOI:10.1016/j.cam.2017.05.018

a. 中国科学院大连化学物理研究所, 分子反应动力学国家重点实验室和理论计算中心, 大连 116023;
b. 中国科学院大学, 北京 100049