The article information
 Yin Huang, Hailin Zhao, Syed Kazim Usman, Ganiyu Ayodele Ajibade, Zhigang Sun
 黄寅, 赵海林, Usman Syed Kazim, Ajibade Ganiyu Ayodele, 孙志刚
 HighAccurate Transparent Boundary Conditions for TimeDependent Quantum Wave Packet Method
 应用于量子含时波包方法的高精度透明边界条件
 Chinese Journal of Chemical Physics, 2020, 33(2): 258262
 化学物理学报, 2020, 33(2): 258262
 http://dx.doi.org/10.1063/16740068/cjcp1909167

Article history
 Received on: September 23, 2019
 Accepted on: January 3, 2020
b. University of Chinese Academy of Sciences, Beijing 100049, China
Solving timedependent 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 timeindependent 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 timedependent wave packet calculation, usually a complex absorbing potential (CAP) of the form
However, when the distribution of the translational energy in the wave packet is of broad range, especially involving ultraslow 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 opensystem 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 semidiscretised form via the
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 pseudospectral 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 ultralow translational energy. At the same time, its implementation is simple and straightforward, and of highorder accuracy. We take the onedimensional motion of a particle to illustrate the proposed method.
The motion of a particle with the mass
$ \begin{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) \end{equation} $  (1) 
where atomic unit (a.u.) is applied. At time
$ \begin{equation} \Psi(x,0) = \exp\left[\frac{(xx_0)^2}{2\sigma_0}\right]\exp(ixk_0) \end{equation} $  (2) 
here
$ \begin{equation} 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] \end{equation} $  (3) 
which is presented in FIG. 1(B) with
The numerical solution of Eq.(1) can only be performed on a finite interval, which we choose as
Our new TBC to dealing the boundary conditions could be derived as the following. The wave packet in the boundary region
$ \begin{equation} \bar{\Psi}(x_i, t) = \sum\limits_{n = 0}^{N}a_n\exp(i\omega_n t)\exp(i k_n x_i) \end{equation} $  (4) 
where
$ \begin{equation} \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) \end{equation} $  (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
$ \begin{equation} \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)] \end{equation} $  (6) 
where
However, we find that when the ChebyshevGauss (CG) quadrature,
$ \begin{equation} \Psi(x_i, t) = \sum\limits_{n = 0}^{N}b_n\exp(i\omega_n t)\cos( k_n x_i) \end{equation} $  (7) 
Then the wave function
$ \begin{equation} \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)] \end{equation} $  (8) 
where
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
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
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
The resonance wave functions at
$ \begin{equation} {\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} \end{equation} $  (9) 
in logarithmic scale with base 10 is plotted against with time in FIG. 2(D). The standard
To strengthen the power of CGTDC, 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 ultralow energy, 0.6352 meV. The potential is given by Eq.(3) with
$ \begin{equation} V_{\rm CAP}(x) = iC\frac{(xx_{1/2})^2}{x_Mx_{1/2}} \end{equation} $  (10) 
with
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 CGTDC results. The oscillations come from the quality of the CAP, which is not long enough. However, the CGTDC results agree with the results using longer CAP (not displayed). Similar conclusion can be drawn with the real part of the
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,
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).
[1] 
R. Kosloff, and D. Kosloff, J. Comput. Phys.
63, 363(1986).
DOI:10.1016/00219991(86)901993 
[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/S00219991(03)001591 
[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 
b. 中国科学院大学, 北京 100049