High Order Split Operators for the Time-Dependent Wavepacket method of Triatomic Reactive Scattering in Hyperspherical Coordinates

Abstract

Since the introduction of a series of methods for solving the time-dependent Schrödinger equation (TDSE) in the 80s of the last centry, such as the Fourier transform, the split operator (SO), the Chebyshev polynomial propagator, and complex absorbing potential, investigation of the molecular dynamics within quantum mechanics principle have become popular. In this paper, the application of the time-dependent wave packet (TDWP) method using high-order SO propagators in hyperspherical coordinates for solving triatomic reactive scattering was investigated. The fast sine transform was applied to calculate the derivatives of the wave function of the radial degree of freedom. These high-order SO propagators are examined in different forms, i.e., TVT (Kinetic–Potential–Kinetic) and VTV (Potential–Kinetic–Potential) forms with three typical triatomic reactions, H + H ${}_2$, O + O ${}_2$ and F + HD. A little difference has been observed among the performances of high-order SO propagators in the TVT and VTV representations in the hyperspherical coordinate. For obtaining total reaction probabilities with 1% error, some of the S class high-order SO propagators, which have symmetric forms, are more efficient than second order SO for reactions involving long lived intermediate states. High order SO propagators are very efficient for obtaining total reaction probabilities.

1. Introduction

Modern methods to solve the time-dependent Schrödinger equation play an important role in the description of atomic and molecular processes [1,2,3,4,5]. Especially, the time-dependent wave packet method for solving reactive scattering processes has become more and more popular due to its numerical scaling advantages. Usually, the Jacobi coordinate is applied for a reactive scattering process, since in it, the Schrödinger equation has a simple form. However, the time-dependent wave packet method using the Jacobi coordinate for a reactive scattering has two primary drawbacks [6]: First, one set Jacobi coordinate is only optimal to represent one of the arrangements. In order to extract state-to-state information, one suffers from the coordinate problem [7,8]. As far as the treatment of products with three free atoms is concerned, the Jacobi coordinates are not optimal choices. In contrast, the hyperspherical coordinate deals with all arrangement channels simultaneously and equally [6,9,10] which is capable of treating all the arrangement channels efficiently using only single propagation. Recently, Zhao et al. developed an efficient interaction-asymptotic region decomposition (IARD) method, where the adiabatically adjusting, principal axes hyperspherical (APH) coordinate presented by [9] was applied for the interaction region, but the corresponding Jacobi coordinates were applied for the asymptotic regions. The IARD method is very efficient for dealing with the state-to-state reactive scattering process using the time-dependent wave packet method [11].

In a numerical simulation of the quantum reactive scattering processes by TDWP-method, the efficiency strongly depends on the two main aspects; (i) the coordinate system and the corresponding grid representation and (ii) the time propagator to evolve the wave packet. Often, these two aspects are closely dependent and one needs to carefully design the whole numerical scheme.

Time-dependent methods are very easy to implement and they have many general applications, i.e., molecular reactive dynamics, prediction of laser atom or molecule interactions, photo-dissociation processes, ultra-cold reactions, etc. The dynamics on which we perform quantum calculations are very helpful for our understanding, through which we have accumulated a large amount of knowledge about micro-mechanisms about chemical reactions, such as the quantum bottleneck state over the transition state, chemical reaction and the geometric phase phenomenon etc.

The numerical grid methods by using the fast Fourier transform (FFT) and its mapped form, pioneered by Kosloff and his co-workers, have also proven to be very efficient and convenient for quantum molecular dynamics studies [12,13,14,15,16]. The other well known method for the solution of the molecular Schrödinger equation is the discrete variable representation (DVR) method [17,18,19,20]. Both of these two method are very effective and will be applied in the present work.

Many efficient wave packet propagation methods have been proposed in the past years, too. The most efficient and accurate propagator was the Chebyshev polynomial expansion, proposed by Tal-Ezer and Kosloff [21]. The second order SO (SSO) method [22,23] is another one of the most popular propagators. One of the interesting feature of the SSO is that it conserves the norm of the wave packet, even when a large time step is used. As a result, the propagation is exceedingly stable. The later introduced Chebyshev real wavepacket (CRWP) [24,25] is also very popular in quantum molecular dynamics field.

For reactive scattering processes, Sun et al. [26] have found that when wave packet is propagated on a flat Potential energy surface using Jacobi coordinate, the SSO appeared to be more efficient than the CRWP propagator, whereas when PES with deep potential is considered, the CRWP will give more efficient results than the SSO propagator. However, the CRWP requires a large number of iterations to obtain fully converged scattering informations for deep potentials. Recently, many groups [27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42] have investigated the application of high-order SOs for solving the Schrödinger equation. Most recently Sun et al. [43,44] draw significant numerical investigations with the high-order SOs and presented several typical tri-atomic reactive scattering calculations. The splitting integrator in either exponential VTV (potential-kinetic-potential splitting) or TVT (kinetic-potential-kinetic splitting) forms were employed in the Jacobi coordinates. They found that, generally, the high-order SO gives more efficient results in TVT form as compared to VTV version. The best high-order SOs in the Jacobi coordinate in most cases are more efficient than the SSO.

The Hamiltonian operators in the APH coordinates are very different from those in the Jacobi coordinates, thus the performances of the different high-order SOs in the APH coordinate should be very different. This motivated us to investigate how high-order SO for reactive scattering calculation works with the APH coordinates. In the literature, there are many different types of hyperspherical coordinates; for details see [45,46,47,48]. The main difference among the APH coordinate and other hyperspherical coordinates [45,46,47,48] lies in the selection of kinematic-angle and orientation of the body-fixed(BF) axes [9]. In the present work, we would only focus on the APH coordinate. The contents of the paper are further divided as follows: Section 2 contains a detailed description and derivation of the theoretical methods for the APH coordinates. Section 3 contains the numerical investigation on performance of the high-order SOs with three typical H + H ${}_2$, F + HD and O + O ${}_2$ reactions. Finally, Section 4 concludes our present work.

2. Theory: Coordinates System, Hamiltonian and Split Operators

2.1. Theory: Mass Scaled Jacobi Coordinate and Initial WavePacket

Let A, B and C be system of three atoms with masses $m_{\alpha }$ where $\alpha$ = A, B, C and positions $P_{\alpha }$ w.r.t SF-axis. Total mass M, reduced mass $\mu$, and the scaling factor $d_{\alpha }$ for the three atoms are defined as follow:

$$M=\sum _{\alpha =A}^C{m}_{\alpha },$$

(1a)

$$\mu =\frac{m_A{m}_B{m}_C}{M},$$

(1b)

$$d_{\alpha }=\frac{m_{\alpha }}{\mu }\left(1-\frac{m_{\alpha }}{M}\right),$$

(1c)

With above three relations, we can now define the mass scaled Jacobi (MSJ) coordinates [45,46]

$$s_{\eta }=P_{\eta +2}-P_{\eta +1}$$

(2a)

$$S_{\eta }=P_{\eta }-\frac{m_{\eta +1}{P}_{\eta +1}+m_{\eta +2}{P}_{\eta +2}}{m_{\eta +1}+m_{\eta +2}}$$

(2b)

$$r_{\eta }=d_{\alpha }^{-1}{s}_{\eta }$$

(2c)

$$R_{\eta }=d_{\alpha }{S}_{\eta }$$

(2d)

where, subscripts $\eta$, $\eta +1$ and $\eta +2$ representing cyclic permutation of atoms A, B and C. $r_{\eta }$ represents the vector from particle B to C, $R_{\eta }$ position vector from center of BC to A. One advantage of the MSJ coordinates is simple orthogonal transformations among different sets are “kinematic-rotations” given by angle ${\chi }_{\eta +1,\eta }$

$$\left(\begin{array}{c}{R}_{\eta +1}\\ r_{\eta +1}\end{array}\right)=U({\chi }_{\eta +1,\eta })\left(\begin{array}{c}{R}_{\eta }\\ r_{\eta }\end{array}\right)$$

(3)

where the transformation matrix U is given by

$$\left(\begin{array}{c}{R}_{\eta +1}\\ r_{\eta +1}\end{array}\right)=\left(\begin{array}{cc}cos({\chi }_{\eta +1,\eta })& sin({\chi }_{\eta +1,\eta })\\ -sin({\chi }_{\eta +1,\eta })& cos({\chi }_{\eta +1,\eta })\end{array}\right)\left(\begin{array}{c}{R}_{\eta }\\ r_{\eta }\end{array}\right)$$

(4)

The kinematics angles are negative obtuse angles with following properties.

$${\chi }_{\eta +1,\eta }=-{\chi }_{\eta +1,\eta }$$

(5a)

$${\chi }_{\eta +1,\eta }=0$$

(5b)

$$tan({\chi }_{\eta +1,\eta })=-\frac{m_{\alpha }}{\mu }$$

(5c)

$$cos({\chi }_{\eta +1,\eta })=-\frac{\mu }{m_{\eta +2}{d}_{\eta }{d}_{\eta +1}}$$

(5d)

$$sin({\chi }_{\eta +1,\eta })=-\frac{m_{\alpha }}{\mu }$$

(5e)

These equations lead to the identity

$${\chi }_{AB}+{\chi }_{BC}+{\chi }_{CA}=2\pi$$

(6)

The MSJ coordinates are simply known as Jacobi coordinates due to extensive use in the literature. The SF or BF-axis set can be used to describe the positions of the three atoms. The six Jacobi SF coordinates consist of two orientation angles of each Jacobi vector and magnitude.

The Hamiltonian in the MSJ coordinates for triatomic reactive scattering can be written as

$$\widehat{\mathbf{H}}=-\frac{{\hslash }^2}{2\mu }\frac{1}{R_{\eta }}\frac{{\partial }^2}{\partial R_{\eta }^2}{R}_{\eta }-\frac{{\hslash }^2}{2\mu }\frac{1}{r_{\eta }}\frac{{\partial }^2}{\partial r_{\eta }^2}{r}_{\eta }+\frac{{\widehat{L_{\eta }}}^2}{2\mu R_{\eta }^2}+\frac{{\widehat{J_{\eta }}}^2}{2\mu r_{\eta }^2}+V,$$

(7)

where $\widehat{L}$ represents orbital angular momentum operator of atom A, $\widehat{J}$ is the rotational angular momentum operator of BC, ${\mu }_{R_{\alpha }}$ is the reduced mass between the center of mass of A and BC, and the total angular momentum is given by $J=\widehat{L}+\widehat{J}$.

For Wave packet calculations, the initial wave packet (IWP) was constructed using MSJ coordinates and then transformation of IWP was done in APH coordinates as described by [6]. For aforementioned problem, the initial wave packet in SF frame $(v_0,j_0,l_0)$ can be simply constructed as a wave function expanded in the SF MSJ coordinates as under:

$${\psi }_{\eta v_0{j}_0{l}_0}^{JM\in }(t=0)=\sum _{\eta v_0{j}_0{l}_0}\frac{1}{R_{\eta }{r}_{\eta }}G(R_{\eta }){\varphi }_{v_0{j}_0}(r_{\eta })|JM{j}_0{l}_0\in \rangle ,$$

(8)

where $|JM{j}_0{l}_0\in \rangle$ represents quantum numbers in SF-representation with the parity $\in ={(-1)}^{j_0+l_0}$, the ro-vibrational eigenfunction of diatom BC is ${\varphi }_{v_0{j}_0(r_{\eta })}$, and the shape of the initial wave function along the translational coordinate is Gaussian function $G(R_{\eta })$ and is given by

$$G(R_{\eta })={\left(\frac{2}{\pi {\sigma }^2}\right)}^{\frac{1}{4}}{\mathrm{e}}^{(-{\frac{(R_{\eta }-R_{\eta }^c)}{\sigma }}^2)}{\mathrm{e}}^{(-k_c{R}_{\eta })}.$$

(9)

Using Equation (8) we can easily transform coordinates system from the MSJ to the APH coordinates.

2.2. Hyperspherical Coordinate for Triatomic Reactive Scattering

Consider the kinematics rotation:

$$\left(\begin{array}{c}{K}_R\\ K_r\end{array}\right)=U({\chi }_{\eta })\left(\begin{array}{c}{R}_{\eta }\\ r_{\eta }\end{array}\right),$$

(10)

where ${\chi }_{\eta }$ continuous variable and its range is $[0,2\pi ]$. The ${\chi }_{\eta }$ differ only in origin for different choices of $\eta$ and are equivalent for different choices of $\eta$.

$${\chi }_{\eta }={\chi }_i-{\chi }_{{\eta }_i},$$

(11)

where ${\chi }_{{\eta }_i}$ are Jacobi kinematics angles given in Equation (5a). The kinematic-angle ${\chi }_{\eta }$ is selected to maximize magnitude of $K_R$, thus, a vector $K_R$ will move towards the vector $R_{\eta }$ for any atom $\eta$ which left other remaining two atom. To obtain maximum value of $K_R$ it can be obtained from:

$$tan(2{\chi }_{\eta })=\frac{2R_{\eta }.r_{\eta }}{{R_{\eta }}^2-{r_{\eta }}^2},$$

(12)

with $\chi$$ϵ$$[-\pi ,\pi ]$.

Transformations from one system $B{F}_{\eta }$ to anothere system $B{F}_{K_R}$ usually consist of rotations ${\beta }_{K_{\eta }}$ about their common y-axis. This ${\beta }_{K_{\eta }}$ is described by the following relation:

$$sin({\beta }_{K_{\eta }})=\frac{r_{\eta }sin({\chi }_{\eta })sin({\mathsf{\Theta }}_{\eta })}{K_R}$$

(13a)

and

$$cos({\beta }_{K_{\eta }})=\frac{R_{\eta }cos({\chi }_{\eta })+r_{\eta }sin({\chi }_{\eta })cos({\mathsf{\Theta }}_{\eta })}{K_R}$$

(13b)

here

$$K_R=\sqrt{r_{\eta }^2{sin({\chi }_{\eta })}^2.{sin({\mathsf{\Theta }}_{\eta })}^2+{(R_{\eta }cos({\chi }_{\eta })+r_{\eta }sin({\chi }_{\eta })cos({\mathsf{\Theta }}_{\eta }))}^2}$$

(13c)

The above system of equations are mapping between Jacobi and APH coordinates and will be used for transformation between these coordinates [9].

In the APH system, the hyperradial $\rho$ represents the radial part, and two angular parts $\theta$, ${\chi }_i$ are internal coordinates. The size of the triatomic system is described by hyperradial $\rho$ and shape of the system is given by hyperangles. The angle $\theta$ is a “bending” angle, it varies the triatomic shape i.e., from an equilateral triangle $\theta =0$ to a collinear geometry $\theta =\pi /2$, and for collinear geometry ${\chi }_i$ represents ratio of $r_{\eta }$ to $R_{\eta }$ for every fixed value $\rho$. These internal coordinates are defined as:

$$\rho =\sqrt{K_r^2+K_R^2}$$

(14a)

$$\theta =\frac{\pi }{2}-2{tan}^{-1}\left(\frac{K_r}{K_R}\right)$$

(14b)

${\chi }_i$ is defined in Equation (11). Here $\rho$$ϵ$$[0,\infty ]$ and $\theta$$ϵ$$[0,\frac{\pi }{2}]$. The internal coordinates deals with all arrangement channels equally.

The relation between internal coordinates of the MSJ coordinates in terms of APH coordinates is given as:

$$R_{\eta }=\frac{\rho }{\sqrt{2}}\sqrt{(1+sin\theta cos2({\chi }_i-{\chi }_{\eta i}))}$$

(15a)

$$r_{\eta }=\frac{\rho }{\sqrt{2}}\sqrt{(1-sin\theta cos2({\chi }_i-{\chi }_{\eta i}))}$$

(15b)

and

$$cos{\mathsf{\Theta }}_{\eta }=\frac{sin\theta sin2({\chi }_i-{\chi }_{\eta i})}{\sqrt{1-{sin}^2\theta {cos}^{2}2({\chi }_i-{\chi }_{\eta i})}}$$

(15c)

and also internal coordinates of APH system in term of the MSJ are as:

$$\rho =\sqrt{R_{\eta }^2+r_{\eta }^2}$$

(16a)

$$tan\theta =\frac{\sqrt{{(R_{\eta }^2-r_{\eta }^2)}^2+4(R_{\eta }^2{r}_{\eta }^2)}}{2R_{\eta }{r}_{\eta }sin{\mathsf{\Theta }}_{\eta }}$$

(16b)

and

$$tan\left[2({\chi }_i-{\chi }_{{\eta }_i})\right]=\frac{2R_{\eta }.r_{\eta }}{{R_{\eta }}^2-{r_{\eta }}^2},$$

(16c)

The Hamiltonian in APH coordinates is defined as:

$$\begin{array}{c}{\mathbf{H}}_{APH}=\frac{{\hslash }^2}{2\mu {\rho }^{5/2}}\frac{{\partial }^2}{\partial {\rho }^2}{\rho }^{5/2}-\frac{{\hslash }^2}{2\mu {\rho }^2}\left[\frac{4}{sin2\theta }\frac{\partial }{\partial \theta }sin2\theta \frac{\partial }{\partial \theta }+\frac{1}{{sin}^2\theta }\frac{{\partial }^2}{{\partial }^2{\chi }_i}\right]\\ +\frac{15{\hslash }^2}{8\mu {\rho }^2}+\frac{1}{\mu {\rho }^2}\left[\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\mathbf{J}+\left({\mathcal{C}}_{\theta }-\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\right){\mathbf{J}}_z^2\right]\\ +\frac{1}{2\mu {\rho }^2}\left[\frac{{\mathcal{A}}_{\theta }-{\mathcal{B}}_{\theta }}{2}({\mathbf{J}}_{+}^2+{\mathbf{J}}_{-}^2)+\hslash {\mathcal{D}}_{\theta }({\mathbf{J}}_{+}+{\mathbf{J}}_{-})\frac{\partial }{\partial {\chi }_i}\right]+V(\rho ,\theta ,\chi )\end{array}$$

(17)

where $V(\rho ,\theta ,\chi )$ is potential energy, $\mathbf{J}$ represents the Total Angular Momentum,

$${\mathbf{J}}_{\pm }={\mathbf{J}}_x\pm i{\mathbf{J}}_y$$

(18)

are the lowering and raising operators, and

$$\begin{array}{cccccccc}\hfill {\mathcal{A}}_{\theta }& =\frac{1}{1+sin\theta },\hfill & \hfill {\mathcal{B}}_{\theta }& =\frac{1}{2{sin}^2\theta },\hfill & \hfill {\mathcal{C}}_{\theta }& =\frac{1}{1-sin\theta },\hfill & \hfill {\mathcal{D}}_{\theta }& =\frac{cos\theta }{{sin}^2\theta }\hfill \end{array}$$

(19)

In calculations, for hyperradial “$\rho$” sine-DVR, for kinematics angle “$\chi$” Fourier DVR, and finite basis representation (FBR) of spherical harmonic basis “$y_{jK}({\theta }_{\alpha })$” used for angular “$\theta$” coordinate.

To propagate the wavepacket in APH coordinates $(\rho ,\theta ,{\chi }_i)$, the wavepacket can be expanded:

$${\mathsf{\Psi }}^{JM\in }(t)=\sum _K\frac{4}{{\rho }^{5/2}}{\overline{D}}_{MK}^{J\in *}({\overrightarrow{\mathsf{\Omega }}}_{\eta }){\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_0),$$

(20)

where ${\overline{D}}_{MK}^{J\in *}({\overrightarrow{\mathsf{\Omega }}}_{\eta })$ represents the parity-adapted normalized rotation matrix, which depend only on Euler angles ${\overrightarrow{\mathsf{\Omega }}}_{\eta }$,

$${\overline{D}}_{MK}^{J\in *}(\overrightarrow{\mathsf{\Omega }})=\sqrt{\frac{2J+1}{8{\pi }^2(1+{\delta }_{K0})}}\left[D_{KM}^{J*}({\overrightarrow{\mathsf{\Omega }}}_{\eta })+\in {(-1)}^{(J+K+\in )}{D}_{-KM}^{J*}({\overrightarrow{\mathsf{\Omega }}}_{\eta })\right],$$

(21)

where ∈ = ${(-1)}^{j+l}$ is the parity of the system, with l is the total orbital angular momentum quantum number and K is the projection of the total angular momentum J on the BF z axis and $D_{MK}^{Jϵ}({\overrightarrow{\mathsf{\Omega }}}_{\eta })$ is the Wigner rotation matrix. The wave function ${\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_o)$ only depends on the internal coordinates $(\rho ,\theta ,{\chi }_i)$ and APH Hamiltonian on the K part of the wave-packet is:

$$\begin{array}{c}{\mathbf{H}}_{APH}{\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_0)=(-\frac{{\hslash }^2}{2\mu }\frac{{\partial }^2}{\partial {\rho }^2}+\frac{15{\hslash }^2}{8\mu {\rho }^2}-\frac{{\hslash }^2}{2\mu {\rho }^2}\left[\frac{4}{sin2\theta }\frac{\partial }{\partial \theta }sin2\theta \frac{\partial }{\partial \theta }+\frac{1}{{sin}^2\theta }\frac{{\partial }^2}{{\partial }^2{\chi }_i}\right]\\ +V(\rho ,\theta ,\chi )+\frac{1}{\mu {\rho }^2}\left[\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\hslash J(J+1)+\left({\mathcal{C}}_{\theta }-\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\right){\hslash }^2{K}^2\right]\\ +\frac{1}{2\mu {\rho }^2}\sum _{K^{\prime }=0}^J〈{\widehat{D}}_{KM}^{J\in *}|\frac{{\mathcal{A}}_{\theta }-{\mathcal{B}}_{\theta }}{2}(J_{+}^2+J_{-}^2)+\hslash {\mathcal{D}}_{\theta }(J_{+}+J_{-})\frac{\partial }{\partial {\chi }_i}|{\widehat{D}}_{K^{\prime }M}^{J\in *}〉){\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_0)\end{array}$$

(22)

i.e., the wave function ${\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_0)$, can be expanded as:

$${\psi }^{JK\in }(\rho ,\theta ,{\chi }_i,t_0)=\sum _{n,m,j}{F}_{nmj}{u}_n(\rho ){\varphi }_m(\chi )y_{j{K}_{\alpha }}(\theta ),$$

(23)

here n, m represent basis-labels, $u(\rho )$ represents basis for hyperradius coordinate, $\varphi (\chi )$ basis for kinematics angles coordinates, and $y_{j{K}_{\eta }}=\sqrt{\frac{(2j+1)}{4\pi }}{d}_{K_{\eta 0}}^j$ is spherical harmonics. reduced Wigner rotational matrix [49] with $K_{\eta }=0$ is $d_{K_{\eta 0}}^j$. $A_{\theta }$, $B_{\theta }$, $C_{\theta }$ and $D_{\theta }$ all are given in the Equation (19). The matrix elements of the Equation (22) can be solved analytically [9]:

$$D_{mat}={\nu }_{J,K}^{+}{\nu }_{J,K+1}^{+}{\delta }_{K^{\prime },K+2}-{\nu }_{J,K}^{-}{\nu }_{J,K-1}^{-}{\delta }_{K^{\prime },K-1}+{(-1)}^{(J+K+\in )}{\nu }_{J,K}^{-}{\nu }_{J,K-1}^{-}{\delta }_{K^{\prime },2-K}$$

(24)

$$\begin{array}{c}\hfill \langle {\widehat{D}}_{KM}^{J\in *}|(J_{+}^2+J_{-}^2)|{\widehat{D}}_{K^{\prime }M}^{J\in *}\rangle ={\hslash }^2\frac{D_{mat}}{\sqrt{(1+{\delta }_{K,0})(1-{\delta }_{K^{\prime },0})}}\end{array}$$

(25)

and

$$\langle {\widehat{D}}_{KM}^{J\in *}|(J_{+}+J_{-})|{\widehat{D}}_{K^{\prime }M}^{J\in *}\rangle =\frac{{\nu }_{J,K}^{+}{\delta }_{K^{\prime },K+1}-{\nu }_{J,K}^{-}{\delta }_{K^{\prime },K-1}+{(-1)}^{(J+K+\in )}{\nu }_{J,K}^{-}{\delta }_{K^{\prime },1-K}}{\sqrt{(1+{\delta }_{K,0})(1-{\delta }_{K^{\prime },0})}}$$

(26)

where ${\nu }_{J,K}^{\pm }$ is given by [9]: ${\nu }_{J,K}^{\pm }=\sqrt{(J\pm K+1)(J\mp K)}$. The results of these equations are asymmetric top and Coriolis coupling coefficients, respectively.

2.3. Split Operators

The SSO is very common nowadays in a molecular quantum dynamics calculation, For completeness, we make a brief introduction about it and its high-order forms.

2.3.1. Second Order Split Operator

The TD form of one-dimensional Schrödinger equation is given by:

$$i\hslash \frac{\partial }{\partial t}\mathsf{\Psi }(\overrightarrow{r},t)=\widehat{\mathbf{H}}\mathsf{\Psi }(\overrightarrow{r},t)=(\widehat{\mathbf{T}}+\widehat{\mathbf{V}})\mathsf{\Psi }(\overrightarrow{r},t),$$

(27)

where V is potential, which depends only on r. Using SSO, the solution of Equation (27) can be written as

$$\mathsf{\Psi }(\overrightarrow{r},t+{\Delta }_t)=e^{-i\widehat{\mathbf{H}}{\Delta }_t}\mathsf{\Psi }(\overrightarrow{r},t)=S_2({\Delta }_t)\mathsf{\Psi }(\overrightarrow{r},t)+O({\Delta }_t^3)$$

(28)

here

$$S_2({\Delta }_t)=e^{i\widehat{\mathbf{H}}{\Delta }_t}=e^{-i\frac{{\Delta }_t\mathbf{V}}{2}}{e}^{-i{\Delta }_t\mathbf{T}}{e}^{-i\frac{{\Delta }_t\mathbf{V}}{2}}$$

(29)

or

$$S_2({\Delta }_t)=e^{i\widehat{\mathbf{H}}{\Delta }_t}=e^{-i\frac{{\Delta }_t\mathbf{T}}{2}}{e}^{-i{\Delta }_t\mathbf{V}}{e}^{-i\frac{{\Delta }_t\mathbf{T}}{2}}$$

(30)

The SSO in Equation (29) is named as the VTV form and in Equation (30) is named as the TVT form, for facilitating the following discussions.

2.3.2. High Order Split Operator

There are a number of ways to develop the high-order SO as reported by [28,29,30,31,32,33,34,35,36,37,38,39,40,41,42], but the simplest and most straight forward way to develop a higher order SO is, product of the lower order integrators with different time steps as presented by [50,51,52] as following

$$S_{2k+2}({\Delta }_t)={\left[S_{2k}({\alpha }_1{\Delta }_t)\right]}^n{\left[S_{2k}({\alpha }_0{\Delta }_t)\right]}^m{\left[S_{2k}({\alpha }_1{\Delta }_t)\right]}^n$$

(31)

here ${\alpha }_1={\frac{1}{(2n-2n{m}^{2k})}}^{\frac{1}{(2k+1)}}$ and ${\alpha }_0=\frac{1-2n{\alpha }_1}{m}$ are $(2k+2)\mathrm{th}$ order split propagators. The efficiency of integrator in Equation (31) comes at higher cost [53], and only $4\mathrm{th}$ order operator have been studied in literature. An alternative way to construct the high-order integrator(s) is as

$$S_n({\Delta }_t)=S_2({\omega }_k{\Delta }_t)...S_2({\omega }_1{\Delta }_t)S_2({\omega }_0{\Delta }_t)S_2({\omega }_1{\Delta }_t)...S_2({\omega }_k{\Delta }_t)$$

(32)

The optimization of the coefficient ${\omega }_k$ is quite complicated, but SOs in Equation (32) is generally more efficient than those in Equation (31). Another way to obtain higher order split integrator is:

$$S_n({\Delta }_t)=e^{-i{\alpha }_{k+1}V}{e}^{-i{\beta }_{k}T}{e}^{-i{\alpha }_{k}V}{e}^{-i{\beta }_{k-1}T}...e^{-i{\alpha }_{2}V}{e}^{-i{\beta }_{1}T}{e}^{-i{\alpha }_{1}V}$$

(33)

Using this procedure, upto $8\mathrm{th}$ order splitting operator propagtors have been investigated [53,54,55,56]. Due to the flexibility of parameters in Equation (33), efficient results can be obtained by using Equation (33).

In our work, to facilitate the following discussion, the high-order SOs in Equations (32) and (33) are termed as A and S-Class/series, respectively, following the work of Sun et al. [43]. And all of the high-order SO propagators are denoted with the same names as given in the work of Sun et al. [43,44]. Numerical error function ${\Delta }_t^s$ is defined for comparing the efficiency of high-order SOs, where ${\Delta }_t^s$ is called effective time step or normalized time step as given by [43], and it is equal to time step (${\Delta }_t$) per stage:

$${\Delta }_t^s=\frac{{\Delta }_t}{N_s}$$

(34)

where $N_s$ indicates the number of the stages. The computational effort for each stage is the same as that for the SSO, thus clearly relates to numerical efficiency of high-order SO.

2.4. Split Operator in the APH Coordinate

For describing a reactive scattering in the APH coordinates using the wave-packet method, we follow our previous formalism [10,11]. The initial wavepacket was first construct in the MSJ coordinates as

$${\psi }_{\eta v_0{j}_0{l}_0}^{JM\in }(t=0)=\sum _{\eta v_0{j}_0{l}_0}\frac{1}{R_{\eta }{r}_{\eta }}G(R_{\eta }){\varphi }_{v_0{j}_0}({\mathbf{r}}_{\eta })|JM{j}_0{l}_0\in \rangle ,$$

(35)

and then was transformed into the APH coordinates, where $|JM{j}_0{l}_0\in \rangle$ represents quantum numbers in SF-representation with the parity $\in ={(-1)}^{j_0+l_0}$, the ro-vibrational eigenfunction of diatom BC is ${\varphi }_{v_0{j}_0(r_{\eta })}$, and the shape of the initial wave function along the translational coordinate is Gaussian function $G(R_{\eta })$ and is given by Equation (9).

Wave packet propagation involves the implementation of sine transform on hyperradius in order to calculate the kinetic energy operators effect on the wavepacket. The interaction of an angular kinetic energy operator on wavepacket can be evaluated in the FBR using associated Legendre polynomials, and for evalution of potential energy, the DVR technique was applied [20,26].

Split operator in TVT form in APH coordinate are written as

$$\begin{array}{cc}\hfill \mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t+{\Delta }_t)& =S({\Delta }_t)\mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t)\hfill \\ & \equiv e^{-i\mathbf{T}{\Delta }_t/2\hslash }{e}^{-i\hslash \mathbf{V}{\Delta }_t/2\hslash }{e}^{-i\mathbf{T}{\Delta }_t/2\hslash }\mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t)\hfill \end{array}$$

(36)

or the VTV form as

$$\begin{array}{cc}\hfill \mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t+{\Delta }_t)& =S({\Delta }_t)\mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t)\hfill \\ & \equiv e^{-i\mathbf{V}{\Delta }_t/2\hslash }{e}^{-i\hslash \mathbf{T}{\Delta }_t/2\hslash }{e}^{-i\mathbf{V}{\Delta }_t/2\hslash }\mathsf{\Psi }(\rho ,\theta ,{\chi }_i,t)\hfill \end{array}$$

(37)

where

$$e^{-i\mathbf{T}{\Delta }_t/2\hslash }=S_2^{\mathbf{T}}({\Delta }_t)\approx e^{-i{\mathbf{T}}_{\chi }{\Delta }_t/2\hslash }{e}^{-i{\mathbf{T}}_{\theta }{\Delta }_t/2\hslash }{e}^{-i{\mathbf{T}}_{\rho }{\Delta }_t\hslash }{e}^{-i{\mathbf{T}}_{\theta }{\Delta }_t/2\hslash }{e}^{-i{\mathbf{T}}_{\chi }{\Delta }_t/2\hslash }$$

(38)

$$\begin{array}{c}\hfill {\mathbf{T}}_{\chi }=\frac{{\hslash }^2}{2\mu {\rho }^2}\left[\frac{1}{{sin}^2\theta }\frac{{\partial }^2}{{\partial }^2{\chi }_i}\right]+\frac{1}{\mu {\rho }^2}\left[\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\mathbf{J}+\left({\mathcal{C}}_{\theta }-\frac{{\mathcal{A}}_{\theta }+{\mathcal{B}}_{\theta }}{2}\right){\mathbf{J}}_z^2\right]\\ \hfill +\frac{1}{2\mu {\rho }^2}\left[\frac{{\mathcal{A}}_{\theta }-{\mathcal{B}}_{\theta }}{2}({\mathbf{J}}_{+}^2+{\mathbf{J}}_{-}^2)+\hslash {\mathcal{D}}_{\theta }({\mathbf{J}}_{+}+{\mathbf{J}}_{-})\frac{\partial }{\partial {\chi }_i}\right],\end{array}$$

(39)

$$\begin{array}{c}\hfill {\mathbf{T}}_{\theta }=\frac{{\hslash }^2}{2\mu {\rho }^2}\left[\frac{4}{sin2\theta }\frac{\partial }{\partial \theta }sin2\theta \frac{\partial }{\partial \theta }\right],\end{array}$$

(40)

$$\begin{array}{c}\hfill {\mathbf{T}}_{\rho }=\frac{{\hslash }^2}{2\mu }\frac{{\partial }^2}{\partial {\rho }^2}+\frac{15{\hslash }^2}{8\mu {\rho }^2}.\end{array}$$

(41)

Higher order SO can be easily implemented by combining series of SSO in both VTV and TVT form. The important point is to carefully design local time step as discussed by [43] for the higher order SO. Due to the energy transformation, time step have upper limit for higher order SOs, which is defined by $(E_{\max }-E_{\min })<{\displaystyle \frac{2\pi }{\Delta t}}$, where $(E_{\max }-E_{\min })$ is the maximal energy distribution of wavepacket.

Finally, we use an absorption potential of the following form [43] to avoid the wave function reaching towards the grid boundary along the $\rho$ degree of freedom:

$$D(R)=i{\Delta }_{t}C\left(\frac{\rho -{\rho }_a}{{\rho }_b-{\rho }_a}\right),{\rho }_a⩽\rho ⩽{\rho }_b={\rho }_{\mathrm{end}},$$

(42)

${\rho }_a$ is the starting point of the absorbing potential region so ${\rho }_a⩽\rho ⩽{\rho }_b$ and C defines the strength of the absorbing potential.

3. Results and Discussion

The numerical error was estimated from the total reaction probabilities, which was calculated by the flux formalism method as

$$P(E)=\frac{1}{\mu }\mathrm{Im}\left[\langle {\mathsf{\Psi }}_f^{JK\in }(E)|\frac{\partial }{\partial \rho }|{\mathsf{\Psi }}_{f^{\prime }}^{JK\in }(E)\rangle \right],$$

(43)

where ${\mathsf{\Psi }}^{JK\in }(E)$ is given by

$${\mathsf{\Psi }}^{JK\in }(E)=\frac{1}{a(E)}\int e^{(iEt)}{\mathsf{\Psi }}^{JK\in }(t)dt,$$

(44)

where $a(E)$ is determined by the initial wavepacket [8,26]. The error was defined as

$$\sigma =\frac{1}{M}\sum _{k=1}^M|P(E)-P^0(E)|/P^0(E),$$

(45)

where M is the number of the collision energies $E_k$ and $P^0(E_k)$ was the converged results calculated with very small time step. For different reactions the value of M and range of energy are given in

Table 1. Reactions, Number of Collision Energies M and Range of energy. Collision energies are evenly distributed in the given range.

Reactions	M	Range of Energies
H + H ${}_2$	700	[0.3, 1.0]
O + O ${}_2$	701	[0.02, 0.16]
F + HD	341	[0.01, 0.035]

.

The Hamiltonian representation in APH coordinates is more complex than in Jacobi coordinates, hence, it requires smart execution of the SO to achieve efficiency and accuracy. In the Jacobi coordinate, the two radial kinetic energy operators can commute with each other, but in the APH coordinates, none of the kinetic energy operator can freely commute with others, which makes its efficient implementation quite difficult. To apply the high-order splitting operator in the APH coordinates for reactive scattering processes, the arrangement of the kinetic energy operators for high-order SO in both the TVT and VTV form is very important. In our strategy, the following forms are used to keep the symmetry and to avoid unnecessary calculations as in Equation (38): $T_{\chi }$ at outer most position, $T_{\theta }$ at central position and then $T_{\rho }$ at inner most position. Using this arrangement for the kinetic operators, the high-order SO in both the TVT and VTV form are implemented in the following calculations.

In order to demonstrate the efficiency of the high-order SOs in either the TVT or the VTV form for a reactive scattering process, H + H ${}_2$, O + O ${}_2$, and F + HD reactions are considered. The aforementioned chemical reactions possess variant dynamic characteristics. H + H ${}_2$ reaction is the easiest direct chemical reaction, while, with lasting resonance states and a deep potential well about 1.1 eV the O + O ${}_2$ chemical reaction is most complex. However, F + HD is direct but have a Fashbach resonances which makes it interesting. Therefore, these reactions gives a comprehensive view for understanding the performance of the high-order SO propagator for triatomic reactions in the APH coordinate.

We plotted total reaction probabilities of the O + O ${}_2$ reaction and F + HD reaction, in Figure 1, calculated using high-order SOs with very small time steps and with time steps which are capable of giving results with errors of about 1%. We could see that in principle, the difference between them is very small and hard to discern by eye. It is also seen that there are many peaks in the probabilities for the O + O ${}_2$ reaction, which implies that there are many long lived resonance intermediate states in the potential well. There is only one peak in the probabilities for the F + HD reaction, which is the ground Feshbach resonance state in the adiabatic potential energy curve D⋯HF ($v^{\prime }$ = 3). This state resides in the van der Waals well in the product channel. The total reaction probabilities of H + H ${}_2$ reaction can be found in many studies and we do not plot them here.

3.1. H + H ${}_2$ Reaction

The adiabatic BKMP2 PES [57] is employed for this reaction in our calculations to examine the performance of the high-order SOs in the APH coordinate.

In Figure 2, the numerical convergences of high-order split propagators in the APH coordinate in the TVT form (A and S class) are presented. It is seen that most of the operators are less efficient than the SSO for the whole error range, except for the S class (4S5b, 4S7) operators for results of high accuracy with small time steps. We note that the propagator named as 4S5 indicates that it is a 4th order S class propagator with four stages. Sometimes, there are more than one such propagators reported in the literature, then we further add a, b and c ... to discern.

The high-order SOs with effective time steps that are smaller than 4 or 5 a.u converge in a way of the SSO, which may suggest that the error mainly comes from the splitting of the kinetic operators. As expected, the SSO remains to be the best choice for this reaction due to the smoothness of the potential, which is similar to the conclusion in the Jacobi coordinate [43,44]. Comparison of the effective time step of the most efficient TVT form propagators in the APH coordinates and the Jacobi coordinates [43,44] at error of 1% are given in

Table 2. Comparison of the most efficient propagators in the TVT form in the Jacobi and the APH coordinates.

	Jacobi	APH	Jacobi	APH
Reactions	A Class/Time Step (a.u)	A Class/Time Step (a.u)	S Class/Time Step (a.u)	S Class/Time Step (a.u)
H + H ${}_2$	4A4b/5.0, 4A6b/3.3	4A4b/11.9, 4A6b/13.1, 6A6/14.91	4S5b/4.5	4S5a/12.1, 4S5b/8.7
O + O ${}_2$	-	4A4b/5.25	-	4S9/16.41
F + HD	4A6a/9.0	4A6a/4.08	4S7, 4S9/5.1	4S5b/5.99

. It is interesting to see that the high-order SOs in the APH coordinate in the TVT form for this reaction are more efficient than those in the Jacobi coordinate.

Similarly, in Figure 3, the results obtained using a high-order SO in the VTV form are presented. The performance of the high-order SOs in the VTV form is very similar to those in the TVT form, except for lower accuracy. There is little difference between the performance of the SSO in the TVT and the VTV form. Thus again, similar to that in the Jacobi coordinate, the SSO is the best choice among all of the SOs in the VTV form for this simple reaction.

3.2. O + O ${}_2$ Reaction

For the O + O ${}_2$ calculations, the SSB PES [58] was applied. In Figure 4, the numerical convergences of high-order split propagators in the APH coordinate in the VTV form of both A and S class are plotted. The results in panel (A) demonstrate that the 4th-order A class propagators are not efficient for the O + O ${}_2$ reaction and their efficiency is lower than the SSO, even for results of high accuracy. They all converge in a way of the SSO, which suggests that the splitting is not so successful for these operators.

In comparison with these 4th-order A class propagators, the 4th-order S class propagators in VTV-form are very efficient and are able to give good results with large effective time steps, as shown in panel (B) of Figure 4. Among them, in order to obtain results with error about 1%, the 4S5b, 4S7, 4S9 and 4S11 are optimal choices to perform quantum calculations with effective time step $\Delta t$ = 15.6 a.u, $\Delta t$ = 16.6 a.u, $\Delta t$ = 16.55 a.u and $\Delta t$ = 16.96 a.u respectively. In contrast, the A class operators were proved to be most suitable choice to perform quantum calculations in the Jacobi coordinates for this reaction [43].

Panels (C) and (D) of Figure 4 presents the numerical convergence of 6th- and 8th-order propagators. Again, we see that the S class SOs are very efficient. Among them, the most efficient 6th-order method is the 6S7 and the most efficient 8th-order is the 8S19. With effective time step $\Delta t$ = 11.3 a.u and 15.7 a.u, they would be able to give reaction probabilities with error less than 1%. We would like to note that, since in the current calculations for the O + O ${}_2$ reaction only collision energy in a very small range is considered, which allows huge total time step, the 8S19 propagator becomes effective. When collision energy in a large range is considered, the total time step becomes smaller and the 8S19 may not be the best choice anymore. From the discussion above, we see that the S class propagators in the VTV form in the APH coordinates are the optimal choices and are much more efficient than those of A class.

The numerical convergence of high-order SO in the TVT form is given in Figure 5. By comparing the results in panel (A) of Figure 5 and Figure 4, it can be seen that the 4th-order A class SOs in the TVT form exhibits a better convergence than those in the VTV form. Anyway, they are still not more efficient than the SSO.

The results in panel (B) of Figure 5 demonstrate that all the examined S class 4th-order propagators in the TVT form are more efficient than the SSO, except the 4S5a, whose efficiency is a little less than the SSO. The best one of them is 4S9, which with effective time step $\Delta t=16.41$ a.u can give results with error less than 1%.

In panels (C) and (D) of Figure 5, the numerical results of the 6th- and 8th-order SO in the TVT form are given. It is seen that the 6S7 is the most efficient one among the 6th-order SOs and the 8S19 is the best one among the 8th-order SOs. With effective time step $\Delta t$ = 9.1 a.u, the 6S7 is able to give probabilities with error less than 1%. In Table 2 and

Table 3. Comparison of the most efficient propagators in the VTV form in the Jacobi and the APH coordinates.

	Jacobi	APH	Jacobi	APH
Reactions	A Class/Time Step (a.u)	A Class/Time Step (a.u)	S Class/Time Step (a.u)	S Class/Time Step (a.u)
O + O ${}_2$	4A6a/20.0	4A4b/5.05	-	4S11/16.96
F + HD	6A8/10.0	4A6a/3.41	4S7/4.5	4S5b/5.94

optimal propagators with the effective time step for giving error of 1% are listed for this reaction, both in the TVT and VTV form, along with the corresponding results in the Jacobi coordinates. It can be concluded that the S class propagators in the VTV form for O + O ${}_2$ reaction in the APH coordinates are very efficient, which is of comparable performance of the best A class SO in the Jacobi coordinate.

3.3. F + HD→ HF + D Reaction

The FXZ PES [59,60] describing the F + HD reaction was used in our following calculations. Figure 6 presents the numerical convergence of high-order SO in the TVT of both A and S class. The results in panel (A) indicate that the examined 4th-order A class SO are less efficient than the SSO, and converges in a second order way. However, the results in panel (B) of Figure 6 indicate that all S class operator in TVT form shows faster convergence as a function of the effective time step. Among them, the 4S5b is the best one. With effective time step as 5.99 a.u, it is able to give reaction probabilities with errors of less than 1%. Comparison of the effective time steps of the best high-order SO in the TVT form for giving error less than 1% in the APH-coordinates and the Jacobi coordinates are presented in Table 2.

Figure 7 gives the numerical convergence for high-order A and S class SOs in the VTV form. It can be seen that these A class operators are inefficient, and their efficiency is only comparable with the SSO for very small time step per stage. They are not as good as the A class operators in the TVT form as shown in panel (A) of Figure 6.

As seen from panel (B) of Figure 7, all examined high-order S class operators in the VTV form are clearly more efficient. The 4S5b is the best one, and with effective time step ${\Delta }_t^s$ = 5.94 a.u, it can give reaction probabilities with errors less than 1%. Comparison of the effective time steps of the best high-order SO in the VTV form for giving error less than 1% in the APH coordinates and the Jacobi coordinates are presented in Table 3.

The 6th-order SOs only works with small time step, and similar behaviour is observed with 8th-order SOs. Numerical convergence are given in Figure 6 and Figure 7.

From the discourse above, we can conclude that the S class high-order SOs in the TVT form are a little more efficient than those in the VTV form, whereas for the A class propagators, the high-order SOs in both TVT and VTV forms have similar numerical convergence behaviours in APH coordinates for the F + HD reaction. The most important observation is that most of the examined high-order SOs for the F + HD reaction exhibit high-order numerical convergence, which is quite unusual. This fact may suggest that most high-order SOs are more efficient than the SSO for calculating the results with high accuracy for the F + HD reaction. This might result from the fact that the Feshbach resonance resides in the product channel, which is the most important region for the reaction. However, for the other two reactions, the most important region for the reaction should be the interaction region. Thus, the error comes from the kinetic operators splitting and being quite different from that in the O + O ${}_2$ reaction.

4. Conclusions

In this study, the performance of a series of high-order SOs for a triatomic reactive scattering process in the APH coordinate is examined. Since the kinetic energy operator in the APH coordinate is more complicated than the kinetic energy operator in the Jacobi coordinate, and the high-order SOs are derived using the one-dimensional model, it is not clear (a priori) if the high-order SOs are still effective. The numerical investigation suggests that the S class high-order SOs are very effective in the APH coordinate for reactions involving resonance, such as O + O ${}_2$ and F + HD reaction. This is different from that in the Jacobi coordinate, where the most effective high-order SOs are the ones of A class. At the same time, we notice that in the APH coordinate, the performance of the high-order SOs in the TVT and VTV form are almost the same. However, in the Jacobi coordinate, the operators in the TVT form are a little better. For the simple direct reaction H + H ${}_2$, the most efficient propagator is proved to be the 2nd-order SO. This is consistent with the results given in the Jacobi coordinate.

In this work, an interesting thing is noticed for the F + HD reaction. Usually, the convergence rate of higher order SOs is of 2nd order at small time steps due to the complicated kinetic operators of the Hamiltonian in the APH coordinate, where their splitting only has a second order convergence. However, irrespective of the time step range, the high-order convergence of the SOs of the F + HD reaction is kept. This suggests that for obtaining results of high accuracy for this reaction, most of the high-order SOs are more efficient than the 2nd order SO.

Currently, there is focused interest on studying ultra-cold reactions. We can expect that the high accuracy property of high-order SOs is attractive and high-order SOs would gain more attention, where high precision results are required.