An improved method for reactive scatterings in ultra-cold conditions using the time-dependent approach

Abstract

A new efficient method for considering the long-range effect of reactive scattering processes in ultra-cold conditions has been developed using the time-dependent quantum wave packet theory, where the initial wave packet could be placed at a position near the interaction region. This is in contrast to previous methods, where the initial wave packet has to be placed far from the interaction region. The new method reduces the numerical effort significantly. Typical reactions, such as S(¹D) + H₂, D⁺ + H₂, and ⁷Li + ⁷Li₂ (v₀ = 1, j₀ = 0), under cold or ultra-cold conditions, are used to demonstrate the numerical efficiency of the new method.

---

1 Introduction

With the development of the molecular cross-beam technique, the collision energy range of a chemical reaction that can be studied has been significantly expanded; from the thermal collision energy^1,2 to several meV,^3–5 and now even to ultra-cold conditions,^6,7 many works have been reported. More and more interesting phenomena have been revealed, such as Feshbach resonance,⁸ resonance-induced quantum tunneling,³ and the statistical effects in ultra-cold chemistry.⁶

Experimentally, only the initial and final states are measurable. The intermediate process of the reaction remains elusive. Reactive scattering theory is important to uncover the dynamic mechanism of the observed interesting phenomena, such as using the quasi-classical trajectory (QCT) theory or the accurate quantum reactive scattering theory.⁹ The QCT theory is fast and applicable for large systems, which provides an intuitive understanding of the reaction.^10,11 However, for reactions involving quantum effects, the QCT theory fails, especially in the realm of cold and ultra-cold chemistry, where quantum effects are usually prominent. Here the quantum reactive scattering theory serves as a crucial tool for comprehending the reactive scattering processes.

The key to quantum reactive scattering theory is solving the Schrödinger equation. Generally, the Schrödinger equation can be solved in two different forms: one based on the time-independent framework,^12,13 and the other based on the time-dependent framework.^14,15 The time-independent method works well in cold and ultra-cold conditions, which usually requires a smaller basis than for thermal energy.^16–18 Anyway, the computational effort of the time-independent method scales as n³, where n represents the total basis functions excluding the ρ degree of freedom, which makes the time-independent method difficult for calculating tetratomic and polyatomic reactions.

It is well known that the computational effort of the time-dependent method scales as N², where N represents the total basis functions or grid points. If the fast Fourier transform technique is employed, the computational effort could be reduced to N log₂ N.¹⁹ However, usually the time-dependent quantum wave packet method is used to calculate reactions at thermal energy.^20–24 In the last several years, some attempts with the time-dependent method have been applied to calculate reactive scattering processes under cold conditions.^25,26 However, since the initial wave packet has to be placed in the true asymptotic region for the reaction under cold and ultra-cold conditions, which is far away from the interaction region in the time-dependent quantum wave packet methods, unlike for thermal conditions, the numerical efficiency is very low.

In this work, an efficient method was developed for reactions under ultra-cold conditions using the time-dependent quantum wave packet theory, avoiding placing the initial wave packet at the true asymptotic region. The method is a combined application of the time-independent method and the time-dependent, where the long range interaction in the asymptotic region is considered by the time-independent method, but the strong interaction region is considered by the efficient time-dependent method. In this way, the advantages of the time-dependent and time-independent method could be fully exploited and the numerical efficiency could be significantly improved. Some typical reactions, S(¹D) + H₂ and D⁺ + H₂ and ⁷Li + ⁷Li₂, under cold conditions, are taken as numerical examples to demonstrate the performance of the proposed new method. The rest of this paper is organized as follows: Section 2 presents the usual triatomic scattering theory and the new calculation strategy. The detailed results are provided in Section 3. The conclusions of this work are provided in Section 4.

2 Time-dependent quantum wave packet theory for reactive scatterings

Here we use a typical triatomic reaction to demonstrate the theoretical details of implementing the time-dependent quantum wave packet theory to calculate reactive scattering at the quantum state-resolved level.

The triatomic reaction scattering, A + BC, could be described using the atom–diatom Jacobi coordinate (R_v, r_v, θ_v, ω_v), where v = α, β, and γ which represent the reactant arrangement channel A + BC, the product arrangement channels B + AC, and C + AB, respectively. R_α(β,γ) is the length of the vector R⃑_α(β,γ) directing from A (B, C) to the center of mass of BC (AC, AB). r_α(β,γ) represents the length of the vector r⃑_α(β,γ) directing from B (C, A) to C (A, B). θ_α(β,γ) is the angle between the two vectors r⃑_α(β,γ) and R⃑_α(β,γ). ω_α(β,γ) expresses the Euler angles orienting R⃑_α(β,γ) in the space-fixed (SF) frame, which defines the scattering angle in the differential cross sections.

Using the reactant Jacobi coordinate (R_α, r_α, θ_α, ω_α), the time-dependent wave function in the body-fixed (BF) frame could be expressed as

We used the sine discrete variable representation (DVR) method to represent the wave function in the radial degree of freedom R_α,²⁷ and R_α,i denotes the grid points. φ_v(r_α) represents the vibrational wave function of diatomic molecule BC. The spherical harmonic functions y_jKα are used to describe the wave function in θ_α. ε is the parity of the system defined as ε = (−1)^j_α+l_α, where j_α represents the rotational angular momentum of diatomic molecule BC and l_α represents the orbital angular momentum. M and K_α represent the projection of total angular momentum J on the SF z-axis and BF z-axis, respectively. D represents the Wigner matrix,²⁸ which depends only on the Euler angles.

The initial packet is constructed in the SF representation in the reactant channel,^29–31 as follows:

Ψ^JMε_α = G(R_α)ϕ_v0j0(r_α)|JMj₀l₀ε〉 (2)

where ϕ_v0j0(r_α) is the eigenvector of the diatomic molecule BC with a specific rovibrational state (v₀, j₀). |JMj₀l₀ε〉 is the eigenvector with total angular momentum J in the SF representation, which can be transformed to the BF representation using Clebsch–Gordan (CG) coefficients.^29,31,32 Here, l₀ is the initial orbital angular momentum. G(R_α) is a Gaussian function in R_α provided bywhere . E₀ represents the central energy of the wave packet, and R_α,0 is the central position of the wave packet.

The collision energy-specified coefficient A(E) of the initial wave packet is usually calculated using

where , H_l is an outgoing Riccati–Hankel function, which is the solution of the following equationby assuming that the long range interaction potential could be ignored by placing the initial wave packet at a position which could be handled using the time-dependent wave packet method. This is usually true for a reaction under thermal reaction conditions, where the initial wave packet could be safely placed at a position with R_α no longer than 20 atomic units (a.u.).

After constructing the initial wave packet, the wave function is propagated from the reactant asymptotic region to the product asymptotic region using the split operator method.²³ In the product region, along with the increasing of R_β, the corresponding r_α in the reactant Jacobi coordinate also increases. This usually leads to a narrower effective region in θ_α, which requires lots of spherical harmonic functions y_jKα, to converge in the reactant Jacobi coordinate. Sometimes, the radial degree of freedom R_α would also increase, depending on the skew angle, which describes the angle between the reactant channel and product channel based on the Jacobi coordinate.³³

That is, an increase in the radial increment R_β in the product Jacobi coordinate would result in significant increase in the basis size and the computational effort. Therefore, for a reaction with product of low translational energy or long-range potential, a single reactant Jacobi coordinate for calculation is very inefficient.

One way to overcome this problem is the interaction-asymptotic region decomposition (IARD) method.^22,24 In this method, the entire region is divided into the interaction region, the reactant asymptotic region, and the product asymptotic region. The optimal Jacobi coordinate is used in the respective asymptotic region. The adjusting principal axes hyperspherical (APH) coordinate^34,35 or the Delves hyperspherical (DH) coordinate^36,37 is used in the interaction region. The IARD method has been verified as very efficient for the reaction F + HD, which involves product of very low translational energy. However, the IARD method is also not so efficient for calculation of the reaction under ultra-cold conditions, such as the ultra-cold Li + Li₂ reaction.

The difficulty for the time-dependent method in calculations of ultra-cold reactions could be conquered in the following way.

Instead of calculating A(E) using eqn (4), the new A(E), labelled as Ā(E), could be calculated as

where H̄_l0 is found by solvingThis equation could be efficiently solved using the Numerov method in the time-independent framework.^30,38,39 Here, V_eff(R_α) represents the effective potential originating from the long-range interaction potential, which is given by V_eff(R_α) = 〈φ_v0j0|JMj₀l₀ε|ΔV|JMj₀l₀ε|φ_v0j0〉 where ΔV = V(R_α,r_α,θ_α) − V_BC(r_α). Similarly, H_lf(k_vfjfR_β/γ) in the product arrangement channel could also be replaced by H̄_lf(k_vfjfR_β/γ) to effectively account for the long-range interaction potential in the ultra-cold conditions.

This improvement of the time-dependent quantum wave packet method for reactive scattering, is easy to implement in the existing time-dependent code, such as with the αDH-IARD method.⁴⁰

For completeness, some details of the αDH-IARD method are briefly presented in the following. In the αDH-IARD method, the respective optimal Jacobi coordinate is used in different asymptotic regions. And, the αDH coordinate, based on the reactant Jacobi coordinate, is used in the interaction region. The internal coordinate (ρ, χ, θ) in this coordinate could be obtained using the following relationship

θ = θ_α (8c)

where d_α represents the mass-scaled factor given by and , where m_A, m_B, and m_C are the masses of atoms A, B, and C, respectively. The Euler angles in the αDH coordinate are the same as those in the reactant Jacobi coordinate. The Hamiltonian in the αDH coordinate could be expressed as

H = T_ρ + T_χ + T_J + V (9)

andwhere j_α represents the rotational angular momentum of the diatomic molecule BC. Correspondingly, the second-order split operator in the αDH coordinate is adopted as

e^−iHΔt ≈ S₂(Δt) = e^−iVΔt/2e^−iT_JΔt/2e^{−iT_χΔt/2}e^−iT_ρΔte^{−iT_χΔt/2}e^−iT_JΔt/2e^−iVΔt/2 (11)

since the second-order split operator is accurate enough for asymptotic region propagation of the wave function,^41,42 even for ultra-cold reactions.

For propagating the wave function in the interaction region, the higher-order split operator has to be included, since the shape of the potential energy surface is more complicated. The higher-order split operator with a symmetrical form is adopted and implemented as⁴²

S_n(Δt) = S₂(α_kΔt)⋯S₂(α₁Δt)⋯S₂(α_kΔt). (12)

In this work, the fourth-order split operator, named as 4S5 in our previous work, is used.⁴² The parameters are given as α₁ = 1/(4 − 4^1/3), α₂ = 1/(4 − 4^1/3), and α₃ = 1 − 4α₁.

Finally, the total reaction probability is obtained using the flux method

andhere, ψ^JMε_β(R⃑_β,r⃑_β;t) represents the wave function in the product asymptotic region for the B + AC channel.

We could see that, comparing with our previous implementation of the αDH-IARD method, the only difference of the new improved method is to replace A(E) with Ā(E) in eqn (14) and then obtain the total or product state-resolved reaction probability. Although it looks like the theoretical “patching” is very little, it could improve the numerical efficiency significantly.

In the time-dependent wave packet method for the ultra-cold reactions, the quality of the applied absorption potential is very important, since the translational energy is very low and usually very long range absorption potential has to be applied. In our calculations, the transmission free absorption potential proposed by Manolopoulos was applied, but in a modified form to avoid the unnecessary singular problem.⁴³ The transmission free absorption potential has the following form

where c ≈ 2.62206 and . R_α,s and R_α,e represent the starting point and end point of the absorbing region in the original equation.

In the original implementation of the transmission free absorption potential, R_α,e was simply set as the ending grid point in the calculation. However, this would lead to a singular potential in the calculation, which cannot be numerically accounted for and is meaningless. It could be understood from the relationship ΔEΔt < 2π and ΔpΔR_α < 2π. Here, ΔE and Δp represent the energy range of the Hamiltonian matrix and the momentum range of the system, respectively. To address this issue, in practical calculations, R_α,e is set as slightly greater than the actual end point of the grid, according to the chosen ΔR_α and Δt. This action not only avoids the singular potential problem, but also improves the absorption effect – the actual length of the applied absorption potential is longer. It is noted that the transmission free absorption potential is very convenient to use, since it is parameter-free and it does need to be optimized. In the following calculations, the typical length of the applied absorption potential is about 1000.0 bohr to achieve an effective damping of the wave function around the grid boundary.

3 Results and discussion

3.1 One-dimensional model

To understand the merits of the proposed new method more deeply, we first demonstrate it using a one-dimensional (1D) model.

Fig. 1(a) displays the 1D model. The potential energy curves of this 1D model are taken as the effective potential of ⁷Li + ⁷Li₂ (v₀ = 1, j₀ = 0) as a function of R when R is greater than 10 bohr. For R less than 10 bohr, the potential is set as a constant. These two regions are connected with a smooth transition function. It is noted that the distance with negative values simply indicates that the region with R < 0 has no realistic meaning. This is just for numerical convenience to obtain convergent results. Fig. 1(b) displays the reaction probability at collision energy 0.0002 meV, i.e., the probabilities going from the right side to the left side of the potential, as a function of the propagation time, calculated using the new method with R₀ = 40 bohr and the previous time-dependent method with R₀ = 40 and 1235 bohr. The result calculated using the time-independent method is also present (the propagation time is meaningless), with the Numerov method.³⁸

Fig. 1 (a) The 1D model. The green line is the potential energy curve, the red line represents the initial wave packet located at R₀ = 40 bohr in the calculation using the new method, and the blue line represents the initial wave packet located at R₀ = 1235 bohr in the calculation using the previous time-dependent method. (b) The reaction probability at collision energy 0.0002 meV as a function of the propagation time, obtained using the new method with R₀ = 40 bohr and the previous time-dependent method with R₀ = 40 and 1235 bohr. The green line represents the accurate results obtained using the time-independent method (the propagation time is meaningless here).

It is seen that the reaction probability obtained using the previous time-dependent method with R₀ = 1235 bohr agrees well with that of the new method with R₀ = 40.0 bohr, and agrees well with that of the time-independent method (Fig. 1(b)). However, the results obtained using the previous method with R₀ = 40.0 bohr, differ significantly.

The reason is, that in the previous method, when the initial wave packet was placed not far enough away, the weak long-range interaction potential could not be ignored in ultra-cold conditions. It could only be eliminated after placing it in a position very far away in true asymptotic region; the blue line shows the effect on obtaining accurate results, i.e., with R₀ = 1235 bohr.

However, by employing the new method with R₀ = 40.0 bohr, the long-range interaction potential could be incorporated in eqn (6) and (7) by Ā(E). Consequently, the convergent results could be obtained by placing the initial wave packet at a much closer position.

We know that the total propagation time is determined by the complete evolution time of the initial wave packet from the initial position to the final position. In the new proposed method, since the initial wave packet could be located in a much closer position, the total propagation time could be significantly reduced and the computational effort could also be significantly reduced correspondingly.

It is seen in Fig. 1(b) that the reaction probability oscillates as a function of total propagation time, this is mainly due to the remaining evolution of the wave function. The amplitude decreases with the increasing propagation time and finally disappears. The period of the oscillations satisfies the uncertainty principle: ΔT ≈ 2π/E. Therefore, the oscillations in the results by the new proposed method with R₀ = 40.0 bohr disappear much more rapidly than the results by the previous method with R₀ = 1235.0 bohr.

3.2 The S(¹D) + H₂ and D⁺ + H₂ reaction under cold conditions

Since the consideration of the long-range interaction potential in the S(¹D) + H₂ and D⁺ + H₂ reactions is important for obtaining accurate results at relatively low collision energy, we took these two reactions to illustrate the advantages of the proposed new method for reactions under cold conditions.^40,44 The S(¹D) + H₂ is a typical insertion reaction, in the following calculations using the IARD method for this reaction, the αDH coordinate is chosen to describe the reaction in the interaction region, our previous numerical test indicated that the αDH coordinate is more effective than the APH coordinate.⁴⁰ In the calculations, the latest potential energy surface (PES) constructed by Yuan and coworkers is used.⁴⁵ For more details of the calculations, one may refer to our previous work.⁴⁰

Fig. 2 displays the total reaction probabilities of the S(¹D) + H₂ reaction for J = 0 and 2, calculated with the initial wave packet located at different positions of R_α,0 using the previous time-dependent method. It is observed that the total reaction probabilities with the initial wave packet located at R_α,0 = 20 a.u. differ significantly from the results located at R_α,0 = 100 a.u. In particular, the total reaction probabilities with the initial wave packet located at R_α,0 = 20 a.u. are much larger than 1 at low collision energy, as shown in Fig. 2(b). These facts indicate that the long-range effect has to be considered in the S(¹D) + H₂ system, and the central position of the initial wave packet must be placed far away to obtain accurate results using the previous time-dependent method, i.e., at least R_α,0 = 100 bohr. This increases a lot of the computational effort, as compared with the calculations at thermal energy.

Fig. 2 (a) Total reaction probabilities of the S(¹D) + H₂ reaction for J = 0 with the initial wave packet located at R_α,0 = 20.0 bohr and 100.0 bohr using the previous time-dependent method, and with the initial wave packet located at R_α,0 = 20.0 bohr using the proposed new method. (b) The results on the right hand side show the enlarged parts of the left panel. (c) Total reaction probabilities of the S(¹D) + H₂ reaction for J = 2 with the initial wave packet located at R_α,0 = 20.0 bohr and 100.0 bohr using the previous time-dependent method, and with the initial wave packet located at R_α,0 = 20.0 bohr using the proposed new method. (d) The results on the right hand side show the enlarged parts of the left panel.

However, with the implementation of the proposed new method, accurate results could be obtained with R_α,0 = 20 bohr. This demonstrates the power of the proposed new method.

For calculation of the D⁺ + H₂ reaction, the PES of the electronic ground state, constructed by Viegas and colleagues, is used.⁴⁶ This PES analytically describes the long-range interaction potential. The numerical parameters in the asymptotic regions could be found in our previous work.⁴⁴ For the interaction region, different from our previous work using the APH coordinate, the αDH coordinate is utilized here: 95 grid points of the sine DVR method are employed for the hyper-radial ρ degree of freedom in the range [0.01, 8.0] a.u.; 79 grid points of the sine DVR method are used for the hyper-angle degree of freedom χ in the range [0, π/2]; 31 spherical harmonic basis functions with even number (j_max = 60) are used in the angular degree of freedom. For efficient and accurate propagation of the wave function in the interaction region, a high-order split operator is utilized.⁴²

Fig. 3(a) displays the total reaction probability of the D⁺ + H₂ reaction with the initial wave packet located at different positions of R_α,0. It is evident that the total reaction probabilities calculated using the previous time-dependent method with the initial wave packet located at R_α,0 = 20 bohr are quite different from those with the initial wave packet located at R_α,0 = 300 bohr. The results with the initial wave packet located at R_α,0 = 300 bohr converge quite well with respect to the position of the initial wave packet. However, there are still some slight oscillations in the results with R_α,0 = 300 bohr due to the limited propagation time 2 × 10⁸ a.u.

Fig. 3 (a) Total reaction probabilities of the D⁺ + H₂ reaction with initial wave packets located at R_α,0 = 20.0 bohr and 300.0 bohr using the previous time-dependent method, and with the initial wave packet located at R_α,0 = 20.0 bohr using the proposed new method, along with the results calculated using the time-independent ABC code with ρ_max = 302. (b) Total reaction probabilities of the D⁺ + H₂ reaction calculated using the time-independent ABC code with different ρ_max.

However, by applying the proposed new method, the results with the initial wave packet located at R_α,0 = 20 bohr, agree with those of R_α,0 = 300 bohr, very well by ignoring the slight oscillations. This might be surprising but is easily understandable.

For a more comprehensive comparison, the total reaction probabilities obtained by the time-independent ABC code are depicted in Fig. 3(b) with varying ρ_max.⁴⁷ In the calculations, Δρ is set to 0.00625 bohr. The maximum energy E_max and angular basis function j_max, are set as 4.2 eV and 50, respectively. It is observed that the total reaction probabilities vary with different ρ_max. When ρ_max is larger than 202 bohr, the results converge. This is consistent with the time-dependent method. For comparison, the converged results obtained by the time-independent ABC code are also displayed in Fig. 3(a), which agrees with the converged time-dependent method results.

In Fig. 3(b), we observe that the peak at a collision energy of about 0.045 meV varies significantly with different ρ_max. To understand this phenomenon, the time-independent scattering wave function at 0.045 meV is displayed in Fig. 4, calculated using the time-dependent wave packet method. It is observed that the wave function varies drastically with R_α less than 50.0 a.u., and with R_α larger than about 200.0 a.u., the wave function gradually exhibits asymptotic behaviour. Therefore, it is not surprising that ρ_max significantly influences the results at a collision energy of about 0.045 meV.

Fig. 4 Time-independent scattering wave function for the D⁺ + H₂ reaction at a collision energy of 0.045 meV as a function of r_α and R_α where θ_α has been integrated out (a), and as a function of R_α where r_α and θ_α have been integrated out (b).

3.3 The ⁷Li + ⁷Li₂ reaction under ultra-cold conditions

The efficiency and accuracy of the proposed new method for the S(¹D) + H₂ and D⁺ + H₂ reactions suggest that the time-dependent wave packet method may be applied for a realistic ultra-cold reaction of alkali metals. The ⁷Li + ⁷Li₂ reaction plays an important role in cold and ultra-cold conditions. Using magnetically tunable Feshbach resonances, long-lived Bose–Einstein condensates can be created.⁴⁸ In the following, we will show the power of the proposed new method for this reaction under ultra-cold conditions. For the calculations of this reaction, the αDH-IARD method was chosen, and the PES constructed by Soldán and coworkers was used.⁴⁹

The numerical details are as follows. In the reactant asymptotic region, 16 383 grid points are utilized for the sine DVR method in R_α in the range [17.0, 1617.0] bohr. The length of the absorbing potential is 1395.5 a.u. to effectively dampen the wave function with extremely low translation energy. 6 vibrational wave functions of ⁷Li₂ are used in r_α, and 16 spherical harmonic basis functions with even number (j_max = 30) are used in θ_α. In the product asymptotic region, the parameters for R_β and r_β are the same as those for R_α and r_α, and 31 spherical harmonic basis functions are used in θ_β. The second-order split operator is used in the asymptotic region with a time step of Δt = 100 a.u. In the interaction region, 191 grid points for the sine DVR are employed for ρ in the range [5.0, 25.0] a.u., 79 grid points for the sine DVR are used for χ in the range of [0, π/2], and 36 spherical harmonic basis functions with an even number (j_max = 70) are used for the angular degree of freedom. The higher-order split operator 4S5 is adopted in the interaction region with a time step of Δt_i = 100 a.u. for obtaining accurate results.⁴² The total propagation time is 1.0 × 10⁹ a.u. under these cold conditions.

Fig. 5(a) displays the cross section for the ⁷Li + ⁷Li₂ (v₀ = 1, j₀ = 0) reaction with different total angular momentum from 0 to 10 with the proposed new method. The collision energy range is from 0.0002 to 2 meV, which corresponds to 2.32 mK to 23.2 K. For a more in-depth comparison, the results obtained by Cvitaš et al. using the time-independent method are also shown in Fig. 5(b).⁵⁰ It is observed that these two results agree well with each other. We note that the nuclear spin is considered in the time-independent method,⁵¹ while it is not considered in our time-dependent framework. Therefore, these two results cannot be exactly the same.

Fig. 5 (a) The cross sections for the ⁷Li + ⁷Li₂ (v₀ = 1, j₀ = 0) reaction under ultra-cold conditions with total angular momentum J from 0 to 10 using the proposed new method. (b) The cross sections in the same energy range with total angular momentum J from 0 to 10 using the time-independent method obtained by Cvitaš et al.⁵⁰

One may worry about the accuracy of our results, especially for collision energy at 0.0002 meV. It requires a very long propagation time and a very long range absorption potential. However, as in our calculations, indeed a very long range absorption potential (about 1400 a.u.) was applied. At the same time, from Fig. 1, we could see that, with a total propagation time of 1.0 × 10⁹ a.u., the results could have an uncertainty of about 1%. Therefore, the results in Fig. 5 are credible, even though the collision energy looks so low for the time-dependent wave packet method.

4 Conclusions

In this study, we have proposed a strategy, combining the time-independent and time-dependent methods, to improve the numerical accuracy and efficiency for a time-dependent wave packet method considering the long-range interaction potential for reactions under cold or ultra-cold conditions. Compared to the calculations using the usual time-dependent wave packet method where the initial wave packet must be placed at a position very far from the interaction region for a reaction in cold conditions, in the proposed new method, the initial wave packet could be placed much closer. Two typical reactions S(¹D) + H₂ and D⁺ + H₂ under cold conditions, which involves long range interaction potential, and the ⁷Li + ⁷Li₂ (v₀ = 1, j₀ = 0) reaction under near ultra-cold conditions, are taken as the numerical examples to examine the performance of the new method, which clearly demonstrate that the improved new method not only reduces the total propagation time but also saves quite a lot of grid range.

The proposed method shares similarities with that of the time-independent method in the work by Lara and coworkers,^52,53 in which the whole reactive region is split into the inner region and outer region. In their work, the hyperspherical coordinate is used in the inner propagation region, and the Jacobi coordinate is used in the outer propagation region. In this way, the computational efficiency could be significantly improved.

Although the method is proposed and demonstrated with triatomic calculations, it is straightforward to extend to tetratomic and polyatomic reactions. In tetratomic and polyatomic reactions, the main issue with the time-independent method is that the numerical scaling scales as n³, and the computational effort could be substantial. By employing our proposed method, the initial wave packet could be placed much closer, which saves computational effort for the time-dependent method for calculating ultra-cold reactions. We expect that the ultra-cold tetratomic reactions or polyatomic reactions could effectively deal with the time-dependent method in the near future, and the combination of the time-dependent and time-independent theory may become popular for investigating reactive scatterings in cold or ultra-cold conditions.

Data availability

The authors confirm that all data supporting the findings of this study are available within the article.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 21825303). H. Z. thanks the Special Research Assistant Funding Project of Chinese Academy of Science and China National Postdoctoral Program for Innovative Talents.

Notes and references

1. D. M. Neumark, A. M. Wodtke, G. N. Robinson, C. C. Hayden and Y. T. Lee, J. Chem. Phys., 1985, 82, 3045–3066.
2. T. G. Yang, J. Chen, L. Huang, T. Wang, C. L. Xiao, Z. G. Sun, D. X. Dai, X. M. Yang and D. H. Zhang, Science, 2015, 347, 60–63.
3. T. G. Yang, L. Huang, C. L. Xiao, J. Chen, T. Wang, D. X. Dai, F. Lique, M. H. Alexander, Z. G. Sun, D. H. Zhang, X. M. Yang and D. M. Neumark, Nat. Chem., 2019, 11, 744–749.
4. M. Lara, S. Chefdeville, K. M. Hickson, A. Bergeat, C. Naulin, J. M. Launay and M. Costes, Phys. Rev. Lett., 2012, 109, 133201.
5. R. Wild, M. Nötzold, M. Simpson, T. D. Tran and R. Wester, Nature, 2023, 615, 425–429.
6. Y. Liu, M. G. Hu, M. A. Nichols, D. Z. Yang, D. Q. Xie, H. Guo and K. K. Ni, Nature, 2021, 593, 379–384.
7. H. Yang, X. Y. Wang, Z. Su, J. Cao, D. C. Zhang, J. Rui, B. Zhao, C. L. Bai and J. W. Pan, Nature, 2022, 602, 229–233.
8. T. Wang, J. Chen, T. Yang, C. Xiao, Z. G. Sun, L. Huang, D. Dai, X. Yang and D. H. Zhang, Science, 2013, 342, 1499–1502.
9. J. M. Bowman and G. C. Schatz, Annu. Rev. Phys. Chem., 1995, 46, 169–196.
10. J. C. Juanes-Marcos, S. C. Althorpe and E. Wrede, Science, 2005, 309, 1227–1230.
11. P. G. Jambrina, D. Herraez-Aguilar, F. J. Aoiz, M. Sneha, J. Jankunas and R. N. Zare, Nat. Chem., 2015, 7, 661–667.
12. R. Kosloff, Annu. Rev. Phys. Chem., 1994, 45, 145–178.
13. B. Jackson, Annu. Rev. Phys. Chem., 1995, 46, 251–274.
14. D. E. Manolopoulos and D. C. Clary, Annu. Rep. Prog. Chem., Sect. C, 1989, 86, 95–118.
15. P. Honvault and J.-M. Launay, Quantum Dynamics of Insertion Reactions, 2004, pp. 187–215.
16. J. M. Hutson and P. Soldan, Int. Rev. Phys. Chem., 2007, 26, 1–28.
17. G. Quemener and P. S. Julienne, Chem. Rev., 2012, 112, 4949–5011.
18. N. Balakrishnan, J. Chem. Phys., 2016, 145, 150901.
19. G. G. Balint-Kurti, R. N. Dixon and C. C. Marston, Int. Rev. Phys. Chem., 1992, 11, 317–344.
20. S. C. Althorpe, J. Chem. Phys., 2001, 114, 1601–1616.
21. S. C. Althorpe and D. C. Clary, Annu. Rev. Phys. Chem., 2003, 54, 493–529.
22. H. L. Zhao, U. Umer, X. X. Hu, D. Q. Xie and Z. G. Sun, J. Chem. Phys., 2019, 150, 134105.
23. H. Zhao and Z. Sun, Numerical Methods for Solving the Time-Dependent Schrödinger Equation for a Molecular Dynamics Process, 2023, ch. 6, pp. 271–348.
24. H. Zhao and Z. Sun, J. Chem. Theory Comput., 2024, 20, 1802–1810.
25. J. Y. Huang, S. Liu, D. H. Zhang and R. V. Krems, Phys. Rev. Lett., 2018, 120, 143401.
26. B. Buren, M. D. Chen, Z. Sun and H. Guo, J. Phys. Chem. A, 2021, 125, 10111–10120.
27. D. T. Colbert and W. H. Miller, J. Chem. Phys., 1992, 96, 1982–1991.
28. R. N. Zare, Angular Momentum, Wiley, New York, 1988.
29. S. Y. Lin and H. Guo, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 74, 022703.
30. Z. G. Sun, X. Lin, S.-Y. Lee and D. H. Zhang, J. Phys. Chem. A, 2009, 113, 4145–4154.
31. Z. G. Sun, H. Guo and D. H. Zhang, J. Chem. Phys., 2010, 132, 084112.
32. J. Z. H. Zhang and W. H. Miller, J. Chem. Phys., 1989, 91, 1528–1547.
33. S. Gómez-Carrasco and O. Roncero, J. Chem. Phys., 2006, 125, 054102.
34. J. Crawford and G. A. Parker, J. Chem. Phys., 2013, 138, 054313.
35. H. Zhao, X. Hu, D. Xie and Z. G. Sun, J. Chem. Phys., 2018, 149, 174103.
36. L. M. Delves, Nucl. Phys., 1958, 9, 391–399.
37. L. M. Delves, Nucl. Phys., 1960, 20, 275–308.
38. J. L. M. Q. Gonzalez and D. Thompson, Comput. Phys., 1997, 11, 514–515.
39. G. Quéméner, P. Honvault, J. M. Launay, P. Soldán, D. E. Potter and J. M. Hutson, Phys. Rev. A: At., Mol., Opt. Phys., 2005, 71, 032722.
40. H. Zhao, D. Xie and Z. G. Sun, J. Phys. Chem. A, 2021, 125, 2007–2018.
41. Z. G. Sun, S.-Y. Lee, H. Guo and D. H. Zhang, J. Chem. Phys., 2009, 130, 174102.
42. Z. G. Sun, W. Yang and D. H. Zhang, Phys. Chem. Chem. Phys., 2012, 14, 1827–1845.
43. D. E. Manolopoulos, J. Chem. Phys., 2002, 117, 9552–9559.
44. H. Zhao, D. Xie and Z. G. Sun, J. Phys. Chem. A, 2021, 125, 2460–2471.
45. J. C. Yuan, D. He and M. D. Chen, Sci. Rep., 2015, 5, 14594.
46. L. P. Viegas, A. Alijah and A. J. C. Varandas, J. Chem. Phys., 2007, 126, 074309.
47. D. Skouteris, J. F. Castillo and D. E. Manolopoulos, Comput. Phys. Commun., 2000, 133, 128.
48. S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. Riedl, C. Chin, J. H. Denschlag and R. Grimm, Science, 2003, 302, 2101–2103.
49. P. Soldán, M. T. Cvitaš and J. M. Hutson, Phys. Rev. A: At., Mol., Opt. Phys., 2003, 67, 054702.
50. M. T. Cvitaš, P. Soldán, J. M. Hutson, P. Honvault and J. M. Launay, Phys. Rev. Lett., 2005, 94, 033201.
51. G. Quéméner, P. Honvault and J. M. Launay, Eur. Phys. J. D, 2004, 30, 201–207.
52. M. Lara, P. G. Jambrina, A. J. C. Varandas, J. M. Launay and F. J. Aoiz, J. Chem. Phys., 2011, 135, 134313.
53. M. Lara, P. G. Jambrina, F. J. Aoiz and J. M. Launay, J. Chem. Phys., 2015, 143, 204305.

---