An efficient way to incorporate the geometric phase in the time-dependent wave packet calculations in a diabatic representation

The geometric phase (GP) effects arising from the crossing between the ground and exited electronic potential surfaces via a conical intersection (CI) have attracted great interest in recent years.^1–6 The effect was discovered independently by Pancharatnam in 1956 in crystal optics⁷ and by Longuet-Higgins et al. in 1958 in molecular systems⁸ and is now commonly referred to as the GP effect, or Berry’s phase effect, for its generalization by Berry to all adiabatic processes in 1984.⁹ Because the real-valued electronic wave function changes sign when the nuclei transport along a loop encircling a CI, the GP must be included in the nuclear quantum dynamics of a molecular system to make the total wave function single valued at each nuclear geometry, even when the total nuclear energy involved is much lower than that of the CI.

One way to include the GP in the adiabatic framework is to define a complex-valued phase angle determined by the hyperangle that encircles the CI. A vector potential acting on the complex-valued phase angle was then introduced to the kinetic energy operator of the Hamiltonian by Mead and Truhlar¹⁰ in a non-trivial manner. Dynamical calculations using the vector-potential-based time-independent (TID) coupled-channel (CC) approach has been successfully employed for the H + H₂ and its isotope reactions and has predicted GP effects both at thermal temperatures^1,11–14 and in the ultracold regime.^15,16 However, since the computational effort of the CC approach increases as N³ with the number N of the basis states, the TID CC calculations are computationally demanding when many partial waves contribute to the reaction.

In 2005, Althorpe and co-workers proposed to include the vector potential in the Jacobi-coordinate-based time-dependent wave packet (TDWP) calculations.^2,17 Because the TDWP method computationally scales much more favorably than the TID CC method, it succeeded in predicting the shift of the oscillation pattern in product differential cross sections (DCSs) for the H + H₂ reaction caused by the GP effects at collision energies above 1.8 eV. Recently, the efforts of GP on the differential cross section (DCS) were observed for the first time in the H + HD → H₂ + D reaction by Yuan et al.⁴ based on the comparison between high-resolution crossed molecular beam experiment measurements and TDWP calculations with/without the vector potential included, indicating that the TDWP approach is reliable on predicting the GP effects. However, incorporation of the vector potential yields a series of extra terms with singularities in the Hamiltonian, making the TDWP calculation much more expensive than that without the vector potential. As a result, it is challenging to use the TDWP with the vector potential to study reactions other than the H₃ system.

On the other hand, because the GP is implicitly included in the adiabatic-to-diabatic states (AtD) transformation, the GP effects can also be described in non-adiabatic dynamical calculations using the diabatic coupled PESs. Recently, Kendrick proposed a simple diabatization scheme for two-state systems with a CI by utilizing the information of ground and excited electronic PESs and the mixing angle that depends on the nuclear coordinates and can be obtained according to the location of the CI. This so-called 2 × 2 diabatic approach was successfully applied to study the GP effects on the ultracold H/D + HD (⁠$v$ = 4, j = 0) → H/D + HD (⁠$v$′, j′) reaction¹⁸ based on the TID CC approach and the APH coordinates. The reported rate coefficients agree perfectly with those obtained by using the vector potential approach, verifying the accuracy of the 2 × 2 diabatic representation in the description of the GP. It was also found that this 2 × 2 diabatic approach is more efficient than the vector-potential-based method because it eliminates the singularities in the Hamiltonian. Very recently, this approach was also used to study the non-adiabatic effects on ultracold H + H₂ (⁠$v$ = 4–8)¹⁹ and H + HD (⁠$v$ = 5–9) reactions.²⁰ 

In this paper, we present a new and efficient way to include the GP in TDWP calculations based on the 2 × 2 diabatic representation by only using the adiabatic ground electronic PES and the information of the CI. Test calculations on the benchmark H + H₂/HD reactions show that the new approach accurately reproduces the results computed using the vector potential approach and diabatic coupled PESs. Because a non-adiabatic TDWP dynamics calculation with N electronic states takes roughly N times of computational time of the corresponding adiabatic calculation, the new approach takes about two times of that on a single PES, considerably faster than the vector-potential-based TDWP calculation.

The nuclear Hamiltonian for a two-state system in the diabatic representation can be written as

where $T^$ is the kinetic energy operator and the real symmetric diabatic potential matrix W can be calculated as

where λ is the so-called mixing angle and V_g and V_ex represent the potential energy of ground state and excited state, respectively. It was found by Herzberg and Longuet-Higgins²¹ that the phase change of the adiabatic electronic wave function around an intersection of two electronic states can be described as

where α is a suitable polar angle encircling the CI. Interestingly, the definition of α is equivalent to η used in the vector potential approach, so we can directly extract the explicit expression of α from the Jacobi coordinates following the previous theoretical treatment designed for η.^3,14,17,18

Based on this idea, Kendrick constructed diabatic PESs for H₃¹⁹ and its isotopic analogs^18,20 based on the accurate V_g and V_ex. However, the use of accurate V_ex not only requires additional efforts to construct an excited electronic PES for a system but also renders the calculation to include other possible non-adiabatic effects associated with dynamics on V_ex, in addition to the GP effects. In fact, in the vector potential approach, the system is restricted on the ground electronic PES with only the GP effects included in quantum dynamics calculations. Similarly, if we set V_ex sufficiently high in the diabatic treatment, we can also restrict the system in the ground electronic state with only the GP included properly as in the vector potential approach. This is the very idea of the present study, and we call this new method as the diabatic version of vector potential (DV-VP) approach. When the total energy of reaction systems is below that of the CI point, the DV-VP approach should be accurate because the excited electronic state is inaccessible anyway. In addition, the new method makes possible the study of pure GP effects when the total energy of the reaction system is higher than the CI. Since the vector potential approach has a singularity at the CI and is hard to be employed when the system energy is above the CI, the DV-VP approach has the potential to be widely used in the study of pure GP effects in reaction systems containing an energetically accessible CI, such as Li + LiNa²² and K + KRb.²³ 

We take the H + HD/H₂ reaction as an illustrative example to show the accuracy of the DV-VP method, as the hydrogen exchange reaction is the prototype chemical reaction for the study of the GP effect. The TDWP dynamical calculations were performed based on the diabatic potential matrix constructed using the BKMP2 PES.²⁴ The calculations for the H + HD → H₂ + D reaction at high collision energies follows closely the previous diabatic dynamics calculations performed in the body-fixed Jacobi coordinates (R, r, θ).^4,25 We applied a total of 161 sine basis functions (including 127 for the interaction region) in the R range of 0.1–15.0 a₀, a total of 127 sine basis basis functions (including 47 for the asymptotic region) in the r range of 0.2–14.0 a₀, and 110 rotational basis functions to converge calculation results after 800 iterations with Δt = 10. The range of the total angular momentum was taken as 0 ≤ J ≤ 50, and all helicity channels were included. Similar damping functions were employed to prevent the wave packet from reflecting back from the boundaries as in Ref. 25. The state-to-state reaction probabilities were extracted by the coordinate transformation reactant coordinate based wave packet method.^25–27 

The calculations at low temperatures for the H + H₂ (⁠$v$ = 4, j = 0) reaction were performed using our newly developed TDWP method for low energy reactive collisions.²⁸ We applied a total of 8191 sine basis functions (including 63 for the interaction region and 511 for the asymptotic reaction) in the R range of 0.1–819.0 a₀, a total of 50 vibrational basis functions (including 10 for the asymptotic region) in the r range of 0.5–6.0 a₀, and 35 rotational basis functions to converge the calculation results at the lowest collision energy after 10 × 10⁶ iterations with Δt = 5.

To verify that the DV-VP calculation is not sensitive to the energy of the exited electronic state used, we calculated the DCS for the H + HD → D + H₂ (⁠$v$′ = 2, j′ = 3) reaction, with the V_ex set as the exact excited potential energy PES, a fixed inaccessible value of 3.5 eV and 5.5 eV. Figure 1 presents the comparison between the calculated products DCSs as a function of collision energy in the backward direction [panel (a)] and as a function of angle at a collision energy of 2.5 eV [panel (b)]. Since the total energy of the reaction system lies below the CI point when E_c < 2.53 eV, the upper electronic state is energetically inaccessible in the energy regime considered here. It can be clearly seen that the calculations with different V_ex agree perfectly with each other, indicating that the value of V_ex has no influence on the calculations when the upper electronic state is energetically inaccessible. For convenience, the value of V_ex is fixed to 5.5 eV in the following calculations.

In Fig. 2, we compare the DCSs of the H₂ (⁠$v$′ = 1, j′ = 9 and $v$′ = 2, j′ = 3) products from the H + HD → D + H₂ reaction at the collision energy of 1.8 eV calculated by the new method with the diabatic results computed using the newly constructed coupled diabatic PESs, which has been verified to describe the dynamics of the reaction accurately. As can be seen from Fig. 2, the angular distribution patterns from the adiabatic calculations with no GP do not agree with the diabatic results. In contrast, the calculations based on the DV-VP approach accurately reproduce the oscillation pattern of the DCSs computed through diabatic theory, indicating that the DV-VP method accurately incorporated the GP into the TDWP calculations.

Figure 3 presents the comparison of backward scattering DCSs for the H + HD → D + H₂ (⁠$v$′ = 2, j′ = 3) reaction as a function of collision energy with/without the GP included. The GP calculations were carried out using the vector potential approach,⁵ the DV-VP approach, and the newly constructed coupled diabatic PESs. It is clear that the DV-VP results are in good agreement with the calculations using the diabatic PESs and the vector-potential-based approach, verifying that the DV-VP approach is capable of describing the GP in the energy regime up to 2.8 eV. It is not surprising that our theoretical treatment for the excited electronic state still remains valid at a collision energy of 2.8 eV because a previous study had shown that the excited electronic state has a slight influence on the DCSs of the reaction at a collision energy of 2.77 eV.⁴ 

The GP in the hydrogen exchange reaction H + H₂ and its isotopic reactions is also reported to effectively control the reaction rates at cold (T < 1 K) and ultracold temperatures (T < 1 mK). In Fig. 4, we compare the total reaction probabilities for the H + H₂ (⁠$v$ = 4, j = 0) reaction for the total angular momentum J = 0 in a wide range of energies extending to the low energy regime obtained using the DV-VP approach and the adiabatic results without the GP effects. It can be seen that the GP modifies the reaction probabilities of the reaction at all energies considered here, with the GP results considerably larger than the NGP ones. In order to verify the accuracy of the DV-VP approach for low energy collisions, we also show in Fig. 4 the diabatic results of Kendrick calculated using the time-independent quantum mechanical (TIQM) approach.¹⁹ As can be seen, the DV-VP calculations reproduce quite well the TIQM results at all energies, confirming further that the new DV-VP approach is capable of studying reliably the GP effects at low temperatures.

In a brief summary, we proposed a diabatic version of the vector potential approach to include the GP in the TDWP study of a two-state system connected via a CI. The test calculations on the DCS for the H + HD → H₂ + D reaction and on the total reaction probability for the H + H₂ (⁠$v$ = 4, j = 0) reaction at low temperatures showed that the DV-VP approach is very accurate on dealing with the GP effects. The new approach only requires the adiabatic PES of the ground electronic state and the location of the CI as the vector potential approach. Furthermore, it merely increases the computational cost by a factor of 2 compared to the adiabatic dynamics calculation, considerably more efficient than the vector potential approach. Due to its superiority in numerical efficiency, this method is expected to be of great use in studying the GP effects in more complex reaction systems, such as the HO₂ reaction system in which the GP effects has been reported to play an important role.¹⁴ 

We thank B.K. Kendrick for sending us the diabatic results computed using the time-independent method. This work was supported by the National Natural Science Foundation of China (Grant Nos. 21590804 and 21688102) and the Chinese Academy of Sciences (Grant No. XDB17010200).