J. Phys. Chem. A 2003, 107, 7197-7203
7197
Full-Dimensional Wave Packet Studies of Collisional Vibrational Relaxation of Both p- and o-H2† Shi Ying Lin and Hua Guo* Department of Chemistry, UniVersity of New Mexico, Albuquerque, New Mexico 87131 ReceiVed: February 17, 2003; In Final Form: April 10, 2003
We report converged full-dimensional quantum dynamical calculations of the vibrational relaxation in the collision: H2 (V1 ) 1, j1 ) 0,1) + H2 (V2 ) 0, j2 ) 0,1) f H2 (V′1 ) 0, j′1) + H2 (V′2 ) 0, j′2), employing a recent global potential energy surface fitted to a large number of high-level ab initio points. The scattering dynamics is characterized by a time-independent wave packet approach based on the Chebyshev polynomial expansion of Green’s operator, which requires repetitive calculations of the action of the system Hamiltonian onto the propagating wave packet. The full-dimensional Hamiltonian within the coupled-states approximation is discretized in a mixed grid/basis representation with the adaptation of the parity and diatomic exchange symmetry, and its action is efficiently computed in the appropriate representation facilitated by a series of one-dimensional pseudospectral transformations. Scattering involving both p- and o-H2 are studied. Rate constants up to a high temperature (3500 K) are obtained from S-matrix elements and compared with available experimental measurements as well as with previous theoretical results.
I. Introduction As the most abundant molecule in the universe, the hydrogen molecule plays an important role in many areas of astrophysics and astrochemistry. For example, collision-induced energy transfer between H2 molecules and between H2 and other atoms/ molecules is believed to be closely related to various astrophysical phenomena, such as cooling of primordial gas and shock wave-induced heating in interstellar media.1-3 To accurately model the thermal balance and kinetics of such systems, stateto-state rate constants become highly desirable.4 Experimentally, it is rather difficult to measure such quantum-state resolved quantities in these systems, and only a limited set of data exists.5-10 On the other hand, accurate theoretical data require both a reliable potential energy surface (PES) and accurate dynamical treatment of the scattering event. The calculation of the electronic energy of the H2 + H2 system is considered to be straightforward as it contains a small number of electrons and is amenable to high-level ab initio methods with a large basis set. However, it is only recently that such calculations have been carried out on a sufficiently large scale to cover the relevant configuration space.11-13 Analytical fits of these ab initio points have since been developed by Aguado, Suarez, and Paniagua (ASP)14 and by Boothroyd, Martin, Keogh, and Peterson (BMKP).13 The BMKP PES is much more reliable than previous versions as a large number (48 180) of MRD-CI points were used in the fitting of the global potential energy function and it has an estimated uncertainty less than 1 kcal/mol. Hence, this PES provides a valuable testing ground for dynamical models. Dynamical studies of the H2 + H2 system have been rather extensive. Classical15 and semiclassical treatments16-19 of the H2 + H2 collision have been reported, but there is always the risk of missing out on important quantum effects in this simplest diatom-diatom system. The quantum mechanical framework †
Part of the special issue “Donald J. Kouri Festschrift”.
for characterizing collisional dynamics between two diatomic molecules is well-established.20-26 However, earlier applications to the H2 + H2 system have resorted to various approximations, such as the distorted wave,20,27 rigid-rotor,20,22,24,26-28 effective potential,23 and two-state approximations,1,22 as well as asymmetric treatments of the two colliding hydrogen molecules.29-33 Although these approximations appear to be reasonable, little has been done until recently to verify their validity by performing exact full-dimensional calculations. The lack of fulldimensional dynamic studies can probably be attributed to the involvement of six internal degrees of freedom, which represents a significant numerical challenge. Moreover, many previous dynamical studies of this system have used empirical21,23 or low-quality ab initio PESs.34-36 The uncertainties arising from both the dynamical approximations and the inaccurate PES make it difficult for a rigorous comparison with the experimental data or an unambiguous evaluation of the validity of either the PES or the dynamical method used in the calculation. The recent emergence of high-quality global PESs of the H4 system has stimulated some interests in dynamical calculations by quantum mechanical methods with as few approximations as possible. This is important not only for checking the validity of the PES used in the calculation but also for establishing an accurate benchmark for more approximate dynamical models. To this end, Pogrebnya and Clary37 reported extensive fulldimensional close-coupling (CC) calculations of rovibrational inelastic collisions between p-H2 molecules as well as between o-H2 molecules employing the BMKP PES. The calculated rate constants over the 20∼300 K temperature range were found to overestimate experimental measurements. They attributed this disagreement to the strong anisotropy of the BMKP PES and modified the PES to lower the anisotropy. This modification was shown to improve the agreement between theory and experiment. In the meantime, our group has studied the pure rotational38 as well as rovibrational inelasticity39 induced by collisions between two p-H2 molecules by using a time-
10.1021/jp0302156 CCC: $25.00 © 2003 American Chemical Society Published on Web 06/07/2003
7198 J. Phys. Chem. A, Vol. 107, No. 37, 2003
Lin and Guo
independent quantum wave packet approach, also with full dimensionality. Rate constants over a wide range of temperatures (0∼3500 K) were calculated. The trend of our results is in qualitative agreement with experimental measurements and the agreement improves at high temperatures (g500 K). However, significant disagreement was found at low temperatures where the theoretical results significantly overestimate experimental measurements. This observation is consistent with the conclusion of Pogrebnya and Clary.37 In this contribution, we extend our previous studies of the para-para H2 collision to include both ortho-ortho and orthopara collisions. In particular, we include scatterings with nonzero Ω, which were avoided in our earlier work on the para-para H2 collision. The discretization scheme thus needs to be modified to accommodate the change. But otherwise, the dynamical methods used here are quite similar to our earlier work. The calculations yield transition probabilities and crosssections up to a high energy (2.2 eV), which allow the calculation of rate constants up to a rather high temperature (3500 K). This paper is organized as follows. In Section II, the discretization scheme and time-independent Chebyshev wave packet method are briefly reviewed in the context of diatomdiatom collisions. In Section III, results from the dynamical calculations are presented and discussed. A short summary is given in Section IV. II. Theory A. Discretization and Evaluation of H ˆ ψ. For the molecular system studied here, it is convenient to use the diatom-diatom Jacobi coordinates (r0, r1, r2, θ1, θ2, and φ). Here, the first three radial coordinates denote the intermolecular and diatomic internuclear distances. The Jacobi angles (θ1 and θ2) are defined between the diatomic and the intermolecular vectors, while φ is the dihedral angle. The Hamiltonian in these coordinates can be written as below (p ) 1):40
H ˆ )-
1
∂2
2µ0 ∂r
2
+ 2
0
2
hˆ i + ∑ ∑ i)1 i)1
ji2
2µiri2
+
(Jˆ - ˆj1 - ˆj2)2
+
2µ0r02
(Jˆ ) and ˆj12 ) ˆj1 + ˆj2 are given by Ω. The total wave packet is expanded to a parity-adapted basis:
ΨJp )
ψiJpi i j j mΩ|i0〉|i1〉|i2〉|j1j2mΩ:JMp〉 ∑ i i i j j mΩ 01212
where the total angular momentum (J), its projection onto the space-fixed (SF) z-axis (M), and the parity (p ) (1) are all good quantum numbers, j1 and j2 are rotational quantum numbers of the diatoms, and m is the projection of rotational angular momentum j1 onto the BF z-axis. As there is no external field, the quantum number M can be arbitrarily chosen. Thus, we assume M ) 0 and drop it from subsequent equations for simplicity. Our choice is equivalent to the E2 frame discussed by Gatti et al.41,42 The discretization strategy is similar to the recent work of Goldfield and Gray,43 although differences exist in implementation. The treatment of the three radial coordinates is straightforward. As shown above, the radial basis is given in a threedimensional direct product grid representation, indexed by i0i1i2. The action of the first radial kinetic energy operator (KEO) in eq 1 is evaluated by fast (sine) Fourier transform on an equidistant grid.44 The two vibrational reference Hamiltonians in eq 2 are represented on potential-optimized discrete variable representation (PODVR) grids.45,46 The use of PODVR minimizes the grid size and improves computational efficiency. Partial sum is used to achieve quasi-linear scaling in the matrixvector multiplication. In contrast to the radial degrees of freedom, the angular part is complicated and deserves a more thorough discussion. The parity-adapted angular basis used in the expansion in eq 3 is given below:
|j1j2mΩ:Jp〉 ) (2 + 2δΩ,0δm,0)-1/2 × [|JΩ〉|j1mj2Ω - m〉 + p(-1)J|J - Ω〉|j1 - mj2 - Ω + m〉] (4) where the quantum number Ω is restricted to nonnegative values to avoid linear dependency, and in the case of Ω ) 0, m is nonnegative. The basis for the overall rotation is the Wigner J rotational matrix (DΩ,M ) defined in terms of the Euler angles 47 (R, β, γ):
2
V(r0, r1, r2, θ1, θ2, φ) -
Vi(ri) ∑ i)1
(1)
where µ0, µ1, and µ2 are the appropriate reduced masses, Jˆ and ˆji are respectively angular momentum operators corresponding to the overall and diatomic rotations, V is the PES of the system, and Vi(ri) (i ) 1, 2) is the potential energy function of free diatomic molecules. The vibrational reference Hamiltonians hˆ i are given by
hˆ i ) -
1 ∂2 + Vi(ri) 2µi ∂r 2
(i ) 1, 2)
(2)
i
As shown below, the major numerical task in propagating the wave packet is the evaluation of the action of the Hamiltonian on to the propagation state, namely, H ˆ ψ. This is done in a mixed grid/basis representation that minimizes the dimensionality of the wave function and renders efficient computation of H ˆ ψ. In particular, we choose the body-fixed (BF) frame with the intermolecular vector (r b0) as the reference z-axis. This is convenient as projections of both the total angular momentum
(3)
01212
|JΩ〉 )
x
2J + 1 J* D Ω,0(R, β, γ) 8π2
(5)
The internal angular basis is, on the other hand, given as a product of uncoupled spherical harmonics:
|j1mj2Ω - m〉 ) Θj1m(θ1)Θj2Ω-m(θ2)x(1/2π) exp(imφ) (6) where Θjm(θ) are normalized associated Legendre functions with the Condon-Shortley phase convention.48 Note that this is not a direct product basis. The use of uncoupled angular basis |j1j2mΩ:Jp〉, instead of the conventional coupled one |j1j2j12Ω:Jp〉,25,49 is because of the pseudospectral transforms used in calculating the action of PES to the wave packet, as discussed below. These two bases can be readily interconverted. In the angular finite basis representation (FBR) outlined above, the matrix of KEOs is very sparse, albeit not diagonal. In particular, the first two angular KEOs in eq 1 are diagonal with the following matrix elements:
〈j1′j2′m′Ω′:Jp|jˆi2|j1j2mΩ:Jp〉 ) ji(ji + 1)δj1′,j1δj2′,j2δm′,mδΩ′,Ω (7)
Collisional Vibrational Relaxation of H2
J. Phys. Chem. A, Vol. 107, No. 37, 2003 7199
For convenience, the third angular KEO is broken into four parts:
(Jˆ - ˆj1 - ˆj2) ) O ˆd + O ˆ 12 + O ˆ1 + O ˆ2 2
forward direction in the following equation: L(θ1:m)
(8a)
(8b)
O ˆ 12 ) ˆj1+ˆj2- + ˆj1-ˆj2+
(8c)
O ˆ i ) -Jˆ +ˆji- - Jˆ -ˆji+
i ) 1, 2
LR,j(θ:m) ) xwRΘj,m(θR)
(8d)
Lβ,m(φ) ) x1/nφΦm(φβ)
The matrix elements of these operators are given as eq 9a-d ( ) [j(j + 1) - m(m ( 1)]1/2. Similar (below), where λj,m expressions have been derived before by Gatti et al.42 and by Goldfield and Gray.43 However, neither provided explicit matrix elements for a parity-adapted basis. Our results also differ from those of Goldfield and Gray by a phase factor due to the Condon-Shortley convention used in this work. Because of the restriction Ω, Ω′ g 0, the terms associated with p(-1)J in eq 9 survive only for special values of Ω and Ω′ (for example, Ω ) Ω′ ) 0 in eq 9b). Despite the formidable appearance of these matrix elements, the action of KEOs on to the wave packet is quite easy to evaluate. In this work, we impose the coupled-state (CS) approximation, which ignores the Coriolis coupling between different Ω channels.50,51 The CS approximation renders Ω a good quantum number and neglects the off-diagonal Coriolis terms in eq 9c,d. Previous experience has indicated that the CS approximation is a reasonable one when applied to the H2 + H2 system.26 Finally, we discuss the calculation of the action of the potential operator on to the wave packet. It can be shown that at a given radial grid point labeled by i0i1i2, the potential matrix element has the form of eq 10 (below), where
Φ(l) m (φ)
) [π(1 + δm,0)]
-1/2
L(φ)
(11)
where R1, R2, and β respectively denote the indices of the angular grids in the θ1, θ2, and φ coordinates. The transformation matrixes are defined as
where
O ˆ d ) Jˆ 2 + ˆj12 + ˆj22 - 2Jˆ z(jˆ1z + ˆj2z) + 2jˆ1zˆj2z
L(θ2:Ω - m)
i0i1i2Ω ψji10ji21mi2Ω 798 ψRi01i1j2i2mΩ 798 ψRi01i1Ri22Ω m 798 ψR1R2β
j ) |m|, ..., jmax; R ) 1, ..., nθ (12a) β ) 1, ..., nφ
(12b)
where nθ and nφ are the numbers of grid points in the polar and dihedral angles, respectively; φβ ) (β - 1/2)π/nφ, and θR and wR are respectively the abscissas and weights of the GaussLegendre quadrature.52 The transformed wave packet in the direct product grid is then multiplied by the diagonal potential matrix. After that, the resulting wave packet is transformed back to the original FBR by three sequential transformations in the backward direction in eq 11. The above pseudospectral transform strategy for the angular coordinates can be traced back to earlier work by several authors43,53-55 including us.56,57 Due to different nuclear spins, H2 molecules exist in either the para or ortho forms, which are practically noninterconvertable by collisions. To maintain total symmetry, p-H2 can only populate even j states, while o-H2 can only populate odd j states. For para-para or ortho-ortho collisions, one can further simplify the calculation by adapting the diatomic exchange symmetry that recognizes the indistinguishability of the H2 molecules.20,21,23-25 As discussed in detail earlier,57 the existence of such symmetry renders nearly half of the grid redundant. Consequently, one only needs to store the wave function at a grid point (i0, i1, i2, j1, j2, m) for i1 g i2. When needed, values at other grid points can be generated by a simple mapping:
{
ψix0,i2,i1,j2,j1,m ) xψix0,i1,i2,j1,j2,m
sin (mφ) (l ) 1) × cos(mφ) (l = 2)
(13)
where x ) (1. In this scheme, it is the wave function, rather than the Hamiltonian matrix, that is symmetrized, which simplifies the coding. Since the exchange and parity operators commute, there are four symmetry species for these H4 systems. For collisions between p- and o-H2, however, the H2 molecules should be treated as distinguishable moieties with no exchange symmetry.
As the potential matrix depends parametrically on three radial indices, a large amount of memory is required to store the matrix elements. This problem can be efficiently circumvented by the pseudospectral method: the wave packet (ψji10ji21mi2Ω) originally in the angular FBR is transformed to the grid representation (ψRi01i1Ri22Ω β ) by three sequential transformations indicated by the
〈j1′j2′m′Ω′:Jp|O ˆ d|j1j2mΩ:Jp〉 ) [J(J + 1) + j1(j1 + 1) + j2(j2 + 1) - 2Ω2 + 2m(Ω - m)]δj1′,j1δj2′,j2δm′,mδΩ′,Ω
(9a)
ˆ 12|j1j2mΩ:Jp〉 ) [(1 + δΩ′,0δm′,0)(1 + δΩ,0δm,0)]-1/2 [λj+1,mλj-2,Ω-mδj1′,j1δj2′,j2(δm′,m+1δΩ′,Ω + 〈j1′j2′m′Ω′:Jp|O p(-1)Jδ-m′,m+1δ-Ω′,Ω) + λj-1,mλj+2,Ω-mδj1′,j1δj2′,j2(δm′,m-1δΩ′,Ω + p(-1)Jδ-m′,m-1δ-Ω′,Ω)] (9b) - 〈j1′j2′m′Ω′:Jp|O ˆ 1|j1j2mΩ:Jp〉 ) -[(1 + δΩ′,0δm′,0)(1 + δΩ,0δm,0)]-1/2 [λJ,Ω λj1,mδj1′,j1δj2′,j2(δm′,m-1δΩ′,Ω-1 + + + λj1,mδj1′,j1δj2′,j2(δm′,m+1δΩ′,Ω+1 + p(-1)Jδ-m′,m+1δ-Ω′,Ω+1)] (9c) p(-1)Jδ-m′,m-1δ-Ω′,Ω-1) + λJ,K - 〈j1′j2′m′Ω′:Jp|O ˆ 2|j1j2mΩ:Jp〉 ) -[(1 + δΩ′,0δm′,0)(1 + δΩ,0δm,0)]-1/2 × [λJ,Ω λj2,Ω-mδj1′,j1δj2′,j2(δm′,mδΩ′,Ω-1 +
〈j1′j2′m′Ω′:Jp|Vi0i1i2|j1j2mΩ:Jp〉 )
[
+ + λj2,Ω-mδj1′,j1δj2′,j2(δm′,mδΩ′,Ω+1 + p(-1)Jδ-m′,mδ-Ω′,Ω+1)] (9d) p(-1)Jδ-m′,mδ-Ω′,Ω-1) + λJ,K
(1 + δm′,0)(1 + δm,0)
]∑
(1 + δΩ′,0δm′,0)(1 + δΩ,0δm,0)
1/2 2
{[δΩ′,Ω +
l)1
(l) (φ)|Vi0i1i2|Θj1,m(θ1)Θj2,Ω-m(θ2)Φ(l) p(-1)J+l δ-Ω′,Ω]/2〈Θj1′,m′(θ1)Θj2′,Ω′-m′(θ2)Φm′ m (φ)〉} (10)
7200 J. Phys. Chem. A, Vol. 107, No. 37, 2003
Lin and Guo
B. S-Matrix and Chebyshev Propagation. In this work, we follow the correlation function formulation proposed independently by Tannor and Weeks58 and by Kouri et al.59 in calculating the S-matrix elements. To this end, the S-matrix element for a transition from initial (i) to final (f) state at energy E is given as
Sf ri(E) )
i 〈χ |G+(E)|χi〉 2πai(E)af*(E) f
(14)
where the initial (final) state wave packet |χi〉(|χf〉) is normally expressed as a product of an internal state (i, f) eigenfunction and a translational wave packet localized in the asymptotic region. The energy amplitude a(E) ) 〈E|χ〉 depends on the choice of the translational wave packet and can be readily evaluated. The major numerical task in eq 14 is to calculate the action of Green’s operator G+(E) ) (E - H ˆ )-1 on to the initial wave packet. A common approach is to expand the operator in terms of the time propagator.58 Alternatively, a Chebyshev polynomial expansion of G+ can be used.60-64 As discussed in our earlier work,38,39 the S-matrix elements are calculated as a Fourier transform of the Chebyshev correlation functions:
Sf ri(E) )
x
e-ikfr0
f
2kf
πµ0
2H
x1 -
-
∑(2 - δk,0)e
-ikarccos(Enorm)
Ck
comparison, the time propagator cannot be evaluated directly and approximations are inevitable.67 In addition, the Chebyshev wave packet can be propagated exclusively in real arithmetic provided the Hamiltonian is real-symmetric and the initial wave packet is real. In contrast to the necessarily complex wave packet in time propagation, this is sometimes referred as a “real wave packet” propagation.68 C. Cross-Sections and Rate Constants. For given Ω, p, and x, the state-resolved integral cross-section for the transition from an initial state V1j1V2j2j12 to a final state V1′j1′V2′j2′j12′ is given by
σVΩpx (E)) 1′j1′V2′j2′j12′rV1j1V2j2j12
πℵ
(E) (2J + 1)PVJΩpx 1′j1′V2′j2′j12′rV1j1V2j2j12 )
∑ σVΩpx′j ′V ′j ′j
1 1 2 2 12′rV1j1V2j2j12
JgΩ
where kf ) x2µ0(E-Ef). To avoid divergence, it is necessary to normalize both energy and its corresponding operator by upper and lower spectral bounds (Hmax and Hmin) of the Hamiltonian:
Enorm ) (E - H+)/H-
H ˆ norm ) (H ˆ - H+)/H-
(16)
with H( ) (Hmax ( Hmin)/2. The correlation function in eq 15 is the overlap between the final internal state eigenfunction (|φf〉) and the Chebyshev wave packet (|ψk〉 ) Tk|ψ0〉) at an asymptotic intermolecular distance (r0f): Ck ) 〈φf|ψk (r0 ) r0f)〉. Noting that the Chebyshev ˆ norm) ) cos kΘ ˆ ],65 polynomial is a cosine-type propagator [Tk(H this expression bears strong resemblance to the time correlation formulism advanced by Tannor and Weeks.58 Like the timedependent method,66 the Chebyshev approach is energyglobal: a single propagation with the initial wave packet |χi〉 produces a column of the S-matrix at all desired energies, provided the corresponding value of ai(E) is not very small. In practice, the propagation of the Chebyshev wave packet |ψk〉 is carried out with asymptotic damping to enforce the outgoing wave boundary condition:62,63
ˆ norm|ψk〉 - D|ψk-1〉) |ψk+1〉 ) D(2H
(17)
with |ψ0〉 ) |χi〉 and |ψ1〉 ) DH ˆ norm|ψ0〉. The coordinatedependent damping function (D) is related to an energydependent negative imaginary potential and should be chosen to damp the wave packets smoothly at the edges of grid. As in our earlier work,38,39 a Gaussian shaped damping function in the intermolecular coordinate is used at r0 > r0d: D ) e-γ(r0-r0d)2. An advantage of the Chebyshev propagation is that it can be realized efficiently and accurately by a three-term recursion (eq 17), which involves mainly matrix-vector multiplication. In
(J,E) (18)
where kV1j1V2j2 ) [2µ0(E - EV1j1V2j2 Depending on the existence of exchange symmetry, ℵ equals to (1 + δV1V2δj1j2)(1 denotes the transition + δV1′V2′δj1′j2′) or 1. PVJΩpx 1′j1′V2′j2′j12′rV1j1V2j2j12 probability for a given J, Ω, p, and x and is expressed in terms of the S-matrix element: )]1/2.
PVJΩpx (E) ) |δV1V1′δj1j1′δV2V2′δj2j2′ 1′j1′V2′j2′j12′rV1j1V2j2j12 (E)|2 (19) SVJΩpx 1′j1′V2′j2′j12′rV1j1V2j2j12
k)0 E2normai(E)
(15)
∑×
kV21j1V2j2 JgΩ
We note that in the above equations the initial or final rotational state is labeled in terms of coupled angular basis with the index j1j2j12 in which |j1 - j2| e j12 e j1 + j2. It can be readily transformed to the uncoupled angular basis representation used in propagation. Finally, the state-resolved thermal rate constant is obtained by a Boltzmann average of the corresponding integral crosssection over the collision energy (Ec ) E - EV1j1V2j2):
(T) ) kVΩpx 1′j1′V2′j2′j12′rV1j1V2j2j12
(
1 8 kBT πµ0kBT
∫0∞σVΩpx′j ′V ′j ′j
)
1/2
×
(Ec)e-Ec/kBTEc dEc (20)
1 1 2 2 12′rV1j1V2j2j12
where T is temperature and kB is the Boltzmann constant. Usually experimental measurements provide more averaged quantities, for example, cross-sections or rate constants corresponding to the V1′j1′V2′j2′ r V1j1V2j2 transition without the resolution of Ωpx. These quantities can easily be obtained by averaging over all allowed values of Ωpxj12 and summing all allowed j12′. Taking the integral cross-section in eq 18 as an example (similarly for the rate constant in eq 20), the averaged quantity is given as
σV1′j1′V2′j2′rV1j1V2j2(E) ) 1 (2j1 + 1)(2j2
∑ + 1)Ωpxj j
12 12′
wxσ VΩpx (21) 1′j1′V2′j2′j12′rV1j1V2j2j12
where w+ ) 3 and w- ) 3 for collisions between o-H2, and w+ ) 1 and w- ) 0 for collisions between p-H2.24 In other cases where two diatomic molecules are distinguishable, e.g., collisions between p- and o-H2, the sum over x collapses. 2/
1/
III. Results For the collision between p-H2 with the lowest rotational quantum numbers (j1,2 ) 0), some simplifications are possible
Collisional Vibrational Relaxation of H2 as demonstrated in our earlier work.38,39 The most significant one is that all the Ω * 0 components vanish even for J > 0. However, when the collision involves o-H2, whose lowest rotational state is j ) 1, such simplifications do not exist and the calculation is much more demanding. For the ortho-ortho collision with j1 ) j2 ) 1, for example, there are three possible values of j12 ()0, 1, 2). For each j12, there are multiple allowed channels labeled by Ω ) 0, ..., j12. Therefore, propagation of multiple initial wave packets corresponding to different j12 and Ω is required. For the ortho-para collision, the number of propagations needed is less than that of the ortho-ortho collision, but the exchange symmetry is lost. We note that there are two cases for the collision between p- and o-H2: The vibrationally excited o-H2 can be relaxed by the collision with the ground-state p-H2 and vice versa. We differentiate these two cases as ortho-para for the former and para-ortho for the latter. Following Gray and Balint-Kurti,68 we used in our calculations a real initial wave packet: |χi〉 ) Ne-[(r0-r0i)/σ]2/2 cos kir0|φi〉, where |φi〉 represents the internal state eigenfunction and N is the normalization constant. The outgoing part of the initial wave packet does not impact the calculation as long as it is effectively eliminated by the damping. The initial wave packets were launched at r0i ) 12.0 bohr with a Gaussian width σ ) 0.3 bohr and central kinetic energy of 0.5 eV. The final state projection was made at r0f ) 12.0 bohr. The onset of the damping was placed at r0d ) 16.0 bohr with γ ) 0.01 bohr-2. To test the convergence, several calculations with different grid/basis sizes and grid ranges have been performed. The numerical parameters used to obtain the results reported in this work are as follows. For the scattering coordinate (r0), 128 equally spaced grid points extending from 2.0 to 21.0 bohr were used. For each of the two vibrational coordinates ri (i ) 1, 2), four PODVR points were found to be sufficient. The rovibrational eigenvalues and eigenfunctions of the free H2 molecule were calculated variationally and represented in the PODVR (j ) 0) grid. For the angular FBR, we have used j1max ) j2max ) 7 for the ortho-ortho collision and j1max ) 7, j2max ) 8 for the ortho-para collision. Only j ) odd (even) for o-H2 (p-H2) were included in the angular basis. For the angular grid, nθ1 ) nθ2 ) nφ ) 12 points were used for the three angular coordinates. The PES was cut off at 5.4 eV to minimize the spectral range of the Hamiltonian. Typically, 3000 steps of Chebyshev propagation were found to be sufficient. Because the existing experimental data are concerned with vibrational relaxation [(V1 ) 1,V2 ) 0) f (V1′ ) 0, V2′ ) 0)], only the S-matrix elements in this final vibrational channel were calculated. Thus, the results reported below are exclusively related to the vibrational relaxation in the collision. To obtain accurate integral cross-sections, one has to calculate numerous partial wave contributions (see eq 18). The J dependence of the final-state-summed partial cross-section Σf σfri(J, E) is displayed in the upper panel of Figure 1 for the transitions from j1 ) j2 ) j12 ) 1 with Ω ) 1, p ) +1, and x ) +1 at two collision energies. Here, f refers to the rotational states in the V1′ ) V2′ ) 0 channel. Typically, the partial wave contribution rises initially and decays to zero at very large J. The number of partial waves needed to converge the crosssection increases with the collision energy. At the collision energy of 1.6 eV, for example, about 70 partial waves are necessary. It can also be noted from the figure that the dependence of the cross-section on the total angular momentum quantum number (J) is quite smooth. Consequently, one can approximate the J dependence of the partial cross-section by
J. Phys. Chem. A, Vol. 107, No. 37, 2003 7201
Figure 1. Comparison of partial cross-sections from interpolation and explicit calculations. The cross-sections are for the collision between two o-H2 from the initial state j1 ) j2 ) j12 ) 1, Ω ) 1, p ) +1, x ) +1 to all the final rotational states in the vibrationally relaxed channel. The J- and E-dependent curves are obtained by cubic spline interpolation with ∆J ) 6, while (+) represent the explicit results not included in the interpolation.
Figure 2. Final-state-summed partial cross-section for transitions from the initial state j1 ) j2 ) 1, j12 ) 0, Ω ) 0. The even and odd exchange symmetry results are strongly oscillatory.
interpolation in the J space, which avoids the calculation of all the partial wave contributions. Indeed, this scheme works quite well as demonstrated by Figure 1, in which the curves represent results obtained by cubic spline interpolation with ∆J ) 6. As shown in the upper panel, the explicitly computed values (pluses) not included in the interpolation fell nicely on top of the interpolated curves. In the lower panel, the dependence of the cross-section on the collision energy is shown for J ) 4 and 33. Again, the interpolated curves coincide with the explicitly computed cross-sections (pluses) very well. Some special treatments were needed in interpolating the results for Ω ) 0 because the smooth J dependence at Ω * 0 no longer holds. As shown in Figure 2, the final-state-summed partial cross-sections for transitions from j1 ) j2 ) 1, j12 ) 0, Ω ) 0 strongly oscillate. For even exchange symmetry (x ) +1), contributions from even J dominate. Similarly, the odd J components are much larger for odd exchange symmetry (x )
7202 J. Phys. Chem. A, Vol. 107, No. 37, 2003
Lin and Guo
Figure 3. Final-state-summed total cross-section for transitions from the initial state j1 ) j2 ) 1, j12 ) 0, Ω ) 0. The even and odd exchange symmetry results are nearly identical.
-1). Our approach to this special situation is to interpolate separately for J ) even block and for J ) odd block. Interestingly, nearly identical interpolation curves are obtained for the J ) even block in even exchange symmetry and for the J ) odd block in odd exchange symmetry, and vice versa. This can be clearly seen in Figure 2. Unlike the para-para collision where only the even exchange symmetry contributes, both the even and odd exchange symmetries should be included in the ortho-ortho collision, cf. eq 21. However, one needs to compute the results of only one symmetry as the integral cross-section is independent of x even when the partial cross-section is drastically different in the two exchange symmetry species as shown in Figure 2. This is shown in Figure 3, where the final-state-summed integral cross-sections for j1 ) j2 ) 1, j12 ) 0, Ω ) 0 are plotted for both the even and odd exchange symmetries. The phenomenon has previously been noted by Zarur and Rabitz23 in the study of rotationally inelastic H2 + H2 collisions, and also by Huo and Green69 in the study of the N2 + N2 collision. We are not aware of any mathematical proof of this numerical observation. In Figure 4, the calculated thermal rate constants for the vibrational relaxation [(V1 ) 1,V2 ) 0) f (V1′ ) 0, V2′ ) 0)] are displayed for the ortho-ortho (koo), ortho-para (kop), paraortho (kpo), and para-para (kpp) collisions. This figure represents the major results of this work. The calculated rate constants are initial-state-resolved (V1 ) 1, j1,V2 ) 0, j2, with average over j12) and summed over the rotational states in the V1′ ) V2′ ) 0 channel:
kV1′V2′rV1j1V2j2(E) ) 1 (2j1 + 1)(2j2
wxkVΩpx ∑ ′j ′V ′j ′j + 1)Ωpxj j ′,j ′,j ′
1 1 2 2 12′rV1j1V2j2j12
12 12
1
(22)
2
The rate constants are displayed in two panels to better compare with the available experimental data,6-8 which are also plotted in the figure. At the high temperature range, experimental rate constants are available only for collisions between normal H2 (knn), which is some kind of combination of koo, kop, kpo, and kpp. Qualitatively, the calculated rate constant reproduces the overall trend of the experimental measurements spanning over 5 orders of magnitude. Quantitatively, however, the agreement is less promising, particularly at low temperatures. The calculated rate constant overestimates the experimental data at low temperatures by about a factor of 3∼10 and underestimates them at high temperatures by about a factor of 2∼3. The underestimation of the experimental data at high temperatures is reasonable, given the fact that only the lowest initial rotational
Figure 4. Temperature dependence of calculated thermal rate constants (solid lines) for the vibrational relaxation of para-para, ortho-para, para-ortho, and ortho-ortho H2 collisions at low (lower panel) and high (upper panel) temperatures. The open (O)8 and solid circles (b)6 are experimental data for collisions between normal H2, while the squares (0) and triangles (4) are measurements for collisions between p-H2 and between o-H2, respectively.7
states were considered in our calculations. Previous studies17,31,37,39 have indicated significant rate enhancement by rotational excitation of the collision partners. Hence, a better agreement between theory and experiment is expected at high temperatures after the Boltzmann averaging in the theoretical calculations. On the other hand, the Boltzmann averaging will exacerbate the discrepancy between theory and experiment at low temperatures. The above conclusions indicate that the inclusion of o-H2 does not change qualitatively the picture obtained in our earlier work on the para-para collision.39 Our calculations are also consistent with the low-temperature results of Pogrebnya and Clary,37 who used a rather different (CC) numerical method but the same PES. The only dynamical approximation, namely, the CS approximation, is unlikely to be responsible for such large discrepancies with experimental data. One possible source of error is the PES. Indeed, Pogrebnya and Clary37 have attributed the theory-experiment discrepancies to the large anisotropy of the BMKP PES, which stems from the pairwise contributions in the potential energy function. A modified version of the BMKP PES with weaker anisotropy was proposed by these authors and a better agreement was obtained. Although the large anisotropy pointed out by these authors probably represents a genuine artifact of the BMKP PES, their modification was a little superficial. It might be worthwhile to add more ab initio points and perform a better fit in order to confirm the speculation. The theoretical results show that the ortho-ortho rate constant is larger than the para-para one and the difference between them decreases with increasing temperature. Although this feature is in qualitative agreement with experimental observa-
Collisional Vibrational Relaxation of H2 tions, quantitative differences exist. For example, the gap between koo and kpp is much larger than experimental data at most temperatures. The calculated results also show that the ortho-ortho rate constant is always larger than all others and the ortho-para rate constant is significantly larger than both the para-para and para-ortho ones at temperatures below 1200 K. The para-para rate constant is slightly larger or equal to the para-ortho one at this temperature range. At higher temperatures (>1200 K), the para-para rate constant exceeds the other two. Unfortunately, experimental results of rate constants for the ortho-para and para-ortho collisions are not available. In the experimental work of Audibert et al.,7 a relationship between different rate constants is determined within (25%, koo + kpp ) kop + kpo. From an analysis based on the calculated rate constants, we found the relationship koo + kpp ) 1.32(kop + kpo), holds very well over a broad range of temperatures (0∼3500 K). This theoretically determined relationship is in reasonable accord with the experimental observation. IV. Conclusions In this work, we reported accurate full-dimensional quantum dynamical calculations of vibrational relaxation induced by collisions between hydrogen molecules. A time-independent Chebyshev wave packet method and a mixed grid/basis representation were used to compute the S-matrix elements. The approach used in this work has more favorable scaling laws than the traditional CC method; thus it should find more applications in other systems. State-resolved rate constants for both the para-para, ortho-ortho, and ortho-para/para-ortho collisions were obtained up to a high temperature (3500 K) by use of the ab initio based BMKP PES. Calculated rate constants were compared with available experimental measurements and previous theoretical results. The calculated rate constants reproduce qualitatively the experimental trend. However, quantitative comparison, particularly at low temperatures, indicates that the calculated rate constants overestimate the experimental data by as much as an order of magnitude. Such disagreement was also reported recently by Pogrebnya and Clary37 using the same PES. The consistency of the two theoretical calculations rules out the possibility of inaccuracy in dynamical calculations as a possible source of error. Further work on the PES and on the experimental measurements is needed to identify the source of such disagreement. Acknowledgment. This work was supported by National Science Foundation (CHE-0090945). This article is dedicated to Professor Don Kouri to honor his many pioneering contributions to theoretical chemistry. References and Notes (1) Danby, G.; Flower, D. R.; Monteiro, T. S. Mon. Not. R. Astron. Soc. 1987, 226, 739. (2) Sternberg, A.; Dalgarno, A. Astrophys. J. 1989, 338, 197. (3) Le Bourlot, J.; Pineau des Forets, G.; Flower, D. R. Mon. Not. R. Astron. Soc. 1999, 305, 802. (4) Martin, P. G.; Schwartz, D. H.; Mandy, M. E. Astrophys. J. 1996, 461, 265. (5) Kiefer, J. H.; Lutz, R. W. J. Chem. Phys. 1966, 44, 668. (6) Audibert, M. M.; Joffrin, C.; Ducuing, J. Chem. Phys. Lett. 1974, 25, 158. (7) Audibert, M.-M.; Vilaseca, R.; Lukasik, J.; Ducuing, J. Chem. Phys. Lett. 1975, 31, 232. (8) Dove, J. E.; Teitelbaum, H. Chem. Phys. 1974, 6, 431.
J. Phys. Chem. A, Vol. 107, No. 37, 2003 7203 (9) Dove, J. E.; Teitelbaum, H. Chem. Phys. 1979, 40, 87. (10) Kreutz, T. G.; Gelfand, J.; Miles, R. B.; Rabitz, H. Chem. Phys. 1988, 124, 359. (11) Boothroyd, A. I.; Keogh, W. J.; Martin, P. G.; Peterson, M. J. J. Chem. Phys. 1991, 95, 4331. (12) Diep, P.; Johnson, J. K. J. Chem. Phys. 2000, 112, 4465. (13) Boothroyd, A. I.; Martin, P. G.; Keogh, W. J.; Peterson, M. J. J. Chem. Phys. 2002, 116, 666. (14) Aguado, A.; Suarez, C.; Paniagua, M. J. Chem. Phys. 1994, 101, 4004. (15) Mandy, M. E.; Martin, P. G.; Keogh, W. J. J. Chem. Phys. 1998, 108, 492. (16) Cacciatore, M.; Billing, G. D. J. Phys. Chem. 1992, 96, 217. (17) Zenevich, V. A.; Billing, G. D.; Jolicard, G. Chem. Phys. Lett. 1999, 312, 530. (18) Zenevich, V. A.; Billing, G. D.; Jolicard, G. Mol. Phys. 2000, 98, 1691. (19) Zenevich, V. A.; Billing, G. D. J. Chem. Phys. 1999, 111, 2401. (20) Davison, W. D. Discuss. Faraday Soc. 1962, 33, 71. (21) Takayangi, K. AdV. At. Mol. Phys. 1965, 1, 149. (22) Allison, A. C.; Dalgarno, A. Proc. Phys. Soc. 1967, 90, 609. (23) Zarur, G.; Rabitz, H. J. Chem. Phys. 1974, 60, 2057. (24) Green, S. J. Chem. Phys. 1975, 62, 2271. (25) Alexander, M. H.; DePristo, A. E. J. Chem. Phys. 1977, 66, 2166. (26) Heil, T. G.; Green, S.; Kouri, D. J. J. Chem. Phys. 1978, 68, 2562. (27) Roberts, C. S. Phys. ReV. 1963, 131, 209. (28) Schaefer, J.; Meyer, W. J. Chem. Phys. 1979, 70, 344. (29) Flower, D. R.; Roueff, E. J. Phys. B: At. Mol. Opt. Phys. 1998, 31, 2935. (30) Flower, D. R. Mon. Not. R. Astron. Soc. 1998, 297, 334. (31) Flower, D. R.; Roueff, E. J. Phys. B: At. Mol. Opt. Phys. 1999, 32, 3399. (32) Flower, D. R. J. Phys. B: At. Mol. Opt. Phys. 2000, 33, L193. (33) Flower, D. R. J. Phys. B: At. Mol. Opt. Phys. 2000, 33, 5243. (34) Monchick, L.; Schaefer, J. J. Chem. Phys. 1980, 73, 6153. (35) Schwenke, D. W. J. Chem. Phys. 1988, 89, 2076. (36) Billing, G. D.; Kolesnick, R. E. Chem. Phys. Lett. 1993, 215, 571. (37) Pogrebnya, S. K.; Clary, D. C. Chem. Phys. Lett. 2002, 363, 523. (38) Lin, S. Y.; Guo, H. J. Chem. Phys. 2002, 117, 5183. (39) Lin, S. Y.; Guo, H. Chem. Phys. 2003, 289, 191. (40) Zhang, J. Z. H.; Dai, J.; Zhu, W. J. Phys. Chem. 1997, 101, 2746. (41) Gatti, F.; Iung, C.; Menou, M.; Justum, Y.; Nauts, A.; Chapuisat, X. J. Chem. Phys. 1998, 108, 8804. (42) Gatti, F.; Iung, C.; Menou, M.; Chapuisat, X. J. Chem. Phys. 1998, 108, 8821. (43) Goldfield, E. M.; Gray, S. K. J. Chem. Phys. 2002, 117, 1604. (44) Kosloff, D.; Kosloff, R. J. Comput. Phys. 1983, 52, 35. (45) Echave, J.; Clary, D. C. Chem. Phys. Lett. 1992, 190, 225. (46) Wei, H.; Carrington, T. J. Chem. Phys. 1992, 97, 3029. (47) Zare, R. N. Angular Momentum; Wiley: New York, 1988. (48) Condon, E. U.; Shortley, G. H. The Theory of Atomic Spectra; Cambridge: London, 1964. (49) Zhang, D. H.; Zhang, J. Z. H. J. Chem. Phys. 1994, 101, 1146. (50) McGuire, P.; Kouri, D. J. J. Chem. Phys. 1974, 60, 2488. (51) Pack, R. T. J. Chem. Phys. 1974, 60, 633. (52) Light, J. C.; Hamilton, I. P.; Lill, J. V. J. Chem. Phys. 1985, 82, 1400. (53) Corey, G. C.; Lemoine, D. J. Chem. Phys. 1992, 97, 4115. (54) Corey, G. C.; Tromp, J. W. J. Chem. Phys. 1995, 103, 1812. (55) Goldfield, E. M. Comput. Phys. Commun. 2000, 128, 178. (56) Chen, R.; Ma, G.; Guo, H. Chem. Phys. Lett. 2000, 320, 567. (57) Chen, R.; Ma, G.; Guo, H. J. Chem. Phys. 2001, 114, 4763. (58) Tannor, D. J.; Weeks, D. E. J. Chem. Phys. 1993, 98, 3884. (59) Kouri, D. J.; Huang, Y.; Zhu, W.; Hoffman, D. K. J. Chem. Phys. 1994, 100, 3662. (60) Huang, Y.; Zhu, W.; Kouri, D.; Hoffman, D. K. Chem. Phys. Lett. 1993, 214, 451. (61) Huang, Y.; Kouri, D. J.; Hoffman, D. K. Chem. Phys. Lett. 1994, 225, 37. (62) Mandelshtam, V. A.; Taylor, H. S. J. Chem. Phys. 1995, 103, 2903. (63) Mandelshtam, V. A.; Taylor, H. S. J. Chem. Phys. 1995, 102, 7390. (64) Chen, R.; Guo, H. J. Chem. Phys. 1996, 105, 1311. (65) Chen, R.; Guo, H. J. Chem. Phys. 1996, 105, 3569. (66) Zhang, J. Z. H. Theory and Application of Quantum Molecular Dynamics; World Scientific: Singapore, 1999. (67) Kosloff, R. J. Phys. Chem. 1988, 92, 2087. (68) Gray, S. K.; Balint-Kurti, G. G. J. Chem. Phys. 1998, 108, 950. (69) Huo, W. M.; Green, S. J. Chem. Phys. 1996, 104, 7572.
