A time-dependent wave packet approach for reaction and dissociation in H2+H2

6 July 2001

Chemical Physics Letters 342 (2001) 177±184

www.elsevier.com/locate/cplett

A time-dependent wave packet approach for reaction and dissociation in H2 ‡ H2 Daniela di Domenico, Marta I. Hern andez *, Jose Campos-Martõnez Instituto de Matem aticas y Fõsica Fundamental (CSIC), C/ Serrano 123, E-28006 Madrid, Spain Received 27 March 2001; in ®nal form 2 May 2001

Abstract A time-dependent wave packet method to study four center (4C) reactions in competition with collision-induced dissociation is presented and applied to the H2 ‡ H2 system, using the potential of Aguado et al. [J. Chem. Phys. 101 (1994) 4004] and a reduced dimensionality model. The calculated probabilities are better converged than previous timeindependent quantum calculations and, in addition, compare quite satisfactorily with quasiclassical trajectory calculations. Tunneling e€ects near reaction thresholds are, however, important. For the 4C reaction in H2 …v1 ˆ 10† ‡ H2 …v2 ˆ 0†, it is estimated that such e€ects are dominant for temperatures below 1200 K, with decreasing values for more excited initial states or for the dissociation process. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction Some years ago, a quantum-mechanical calculation of the four center (4C) reaction and collision-induced dissociation (CID) in H2 ‡ H2 was reported by Hern andez and Clary [1]. This was the ®rst quantum-mechanical treatment of a 4C reaction, where two bonds break and form in the course of the reaction. A reduced dimensionality approach was proposed where the atoms are restricted to be within a trapezoidal geometry, so that CID and 4C probabilities were obtained for reactants with selected vibrational states, H2 …v1 † ‡ H2 …v2 † using a time-independent hyperspherical coordinate treatment and a realistic potential energy surface [2]. Although several

*

Corresponding author. E-mail address: marta@ima€.cfmac.csic.es (M.I. Hernandez).

computational limitations prevented a satisfactory convergence of the initial-state selected probabilities, several clear trends were found: (a) only reactants with v1 ˆ high and v2 ˆ 0 yielded large probabilities; (b) as v1 increases, lower dynamical thresholds and higher probabilities are found (vibrational enhancement); and (c) dominance of CID over 4C at low energies and a reverse trend at higher energies. In contrast to the well-studied AB ‡ CD ! ABC ‡ D reactions [3,4], 4C reactions [5] usually have quite high barriers, involve more complex mechanisms, and present competing dissociation processes. This makes quantum-mechanical calculations a big theoretical and computational challenge, and nowadays full-dimensional calculations are still far from possible. It is also worth noting that the quantum-mechanical treatment of CID is already a tremendous challenge for triatomic systems [6,7]. The use of reduced dimensionality approaches (as that of [1]) has at least a

0009-2614/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 1 ) 0 0 5 5 8 - 9

178

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

double interest. On the one hand, the most relevant quantum features can be understood and, at the same time, computational methods can be progressively developed to reach a quantum description each time more accurate and richer. On the other hand, full-dimensional quasiclassical trajectory (QCT) calculations [8±11] are indeed possible and provide detailed information of magnitudes of interest such as cross-sections and rate coecients, but the validity of the classical approach is sometimes unknown. The extent of quantum e€ects can be inferred, for instance, from comparisons of restricted dimensionality QCT calculations with equivalent quantum computations. The latter strategy has been recently followed by Ceballos et al. [8]. They carried out both full and reduced dimensionality QCT calculations for H2 ‡ H2 , also using the potential of [2]. They ®rst compared QCT and quantum [1] probabilities ± within the same model± and found that, although the classical probabilities follow the same trends outlined in the ®rst paragraph, a quantitative agreement was not reached. Ceballos et al. then compared reduced and full-dimensional QCT calculations and noticed several similarities (such as the energy dependence of the cross-sections and the vibrational enhancement e€ect) and some differences (a major relevance of CID and appearance of reactive±dissociative processes). They also found that nearly zero impact parameter collisions are most e€ective in promoting the 4C reaction, in agreement with the assumptions of the reduced dimensionality model. Large impact parameters are, on the contrary, quite important for dissociative processes. The somewhat poor agreement between model classical and quantum calculations could be due to the onset of quantum e€ects but also to some computational limitations of the quantum calculations. As discussed in detail in [1], only one of the four symmetry blocks of the Hamiltonian system was considered and, in addition, the asymptotic states converged too slowly using hyperspherical coordinates. To answer the question raised above, we solve here the dynamical problem of [1] with a di€erent quantum±mechanical method, the time-dependent

wave packet approach. Initial-state selected CID and 4C probabilities are obtained by evaluating ¯uxes [12] through surfaces separating reactants and the di€erent product channels. In the hyperspherical method of [1], state-to-state probabilities from all initial reactants states are obtained at once in the same run of the code, but achieving high accuracy for all transitions is quite hard. Once this one is interested in a selected set of initial states of reagents, the wave packet approach (using reagent Jacobi coordinates) can provide accurate results at a lower computational cost. Indeed, we have obtained much better converged probabilities over a larger energy range. With further improvements, this method could allow us to include more degrees of freedom in the calculations. We end this section with a word on the H2 ‡ H2 system. It is fundamental in theoretical chemistry and molecular dynamics, since it is the lightest tetraatomic system. The H2 ‡ H2 dynamics and kinetics is of interest in astrophysics [10], combustion [13] and spacecraft modeling [14]. A summary of experimental results is given in [5] and [15]. Previous dynamical studies are reviewed in [1]. More recently, Mandy et al. [11] have compared the role of H2 , He and H as colliders of target H2 , using the QCT method, and, in other work [16], have analyzed the classical trajectories within several restricted geometries. In addition, Ceballos et al. [9] have extended their earlier study [8] to the computation of rate coecients and product distributions.

2. Time-dependent wave packet method for H2 ‡ H2 In the reduced dimensionality model of [1], four AABB atoms are constrained to be in a plane forming an isosceles trapezium. The allowed con®gurations are described by restricted Jacobi coordinates r1 and r2 , the interatomic distances of A2 and B2 , respectively, and R, the distance between the centers of mass of A2 and B2 . It should be noted that coordinate R ranges from 1 to 1 since con®gurations where R  0 are energetically allowed. The vibrational states of reactants, A2 …v1 † ‡ B2 …v2 †, are well de®ned within the model.

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

179

Table 1 Reaction channels of A2 ‡ B2 within the reduced dimensionality model of [1] Channel initial/inelastic A2 dissociation B2 dissociation Four center

A2 ‡ B2 A ‡ A ‡ B2 A2 ‡ B ‡ B AB ‡ AB

R

r1

r2

Large All All Small

Small Large Small Large

Small Small Large Large

Values of the coordinates where the reactants/products valleys are located are also indicated.

The product channels described by the model are listed in Table 1. The Hamiltonian is Hˆ

h2 o2 2lA2 or12

2 o2 h 2lB2 or22

2 h 2lA2

B2

o2 oR2

‡ V …r1 ; r2 ; R†;

where lA2 , lB2 and lA2 B2 are the reduced masses of A2 , B2 and A2 B2 , respectively, and V is the interaction potential. The potential energy surface of Aguado et al. [2] was used, which is an accurate ®t of a large set of (also accurate) ab initio energies [17]. The coordinates not explicitly accounted for in the model (3 angles) are ®xed to those values preserving the trapezoidal geometry. The initial wave packet is given as …2†

where /v1 …r1 † and /v2 …r2 † are the vibrational eigenstates of the diatomic reactants and v…R† is a general Gaussian wave packet [18] for the translational degree of freedom. The wave packet is propagated in time subject to the time-dependent Schr odinger equation, using the split operator method [19]. Wave packets are absorbed in the asymptotic regions by means of a wave packet splitting algorithm [20,21], using absorbing func2 tions of the form f …q† ˆ exp… b…q qabs † †, for q > qabs and one otherwise, where q  r1 ; r2 or R. Initial-state selected probabilities are obtained by evaluating the time-independent ¯ux through surfaces separating reactants and the di€erent product channels [12,22]. In this way, the dissociation probability of A2 is written as Z 1 Z r2s 2† …E† ˆ dR dr2 J …r1f ; r2 ; R; E†; …3† PvCID…A 1 ;v2 1

0

r2s

1

Z …1†

W…r1 ; r2 ; R; t ˆ 0† ˆ /v1 …r1 †/v2 …r2 †v…R†;

(analogously for the dissociation of B2 ), and the 4C probability as Z 1  Z r2f 4C Pv1 ;v2 …E† ˆ dR dr2 J …r1f ; r2 ; R; E† ‡

r1f

r1s

 dr1 J …r1 ; r2f ; R; E† ;

…4†

where J …r1f ; r2 ; R; E† is the probability current at r1f :   h dWE  J …r1f ; r2 ; R; E† ˆ Im WE …r1f ; r2 ; R† jr ˆr1f ; lA 2 dr1 1 …5† and analogously for J …r1 ; r2f ; R; E†. WE is the timeindependent wavefunction, obtained by Fourier transformation of the wave packet, Z 1 dt eiEt=h W…r1 ; r2 ; R; t†; WE …r1 ; r2 ; R† ˆ …6† a…E† and a…E† is the energy-component of the initial wave packet, a…E† ˆ Ch/v1 /v2 e

ikR

j W…t ˆ 0†i;

…7†

C being a normalization constant such that hWE jW0E i ˆ 2phd…E E0 †. The choice of the ¯ux surfaces becomes clear with the help of Table 1: r1f and r2f are suciently large values of the interatomic coordinates, such that diatomic bound states present a negligible probability density. To obtain the probability of dissociation of A2 , for instance, it is assumed that all the ¯ux passing through the dividing surface r1f at a suciently small r2 value (from 0 to r2s , such that the B2 bond is not broken) will end up in the A ‡ A ‡ B2 asymptotic valley. Analogously, it is assumed that the ¯ux through a surface where both r1 and r2 are suciently large will end up in

180

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

the AB ‡ AB asymptotic region (as far as the full dissociation channel, A ‡ A ‡ B ‡ B, is closed, as in the present case; otherwise, a more detailed sectioning of the ¯ux surfaces would be needed). A proper separation of CID and 4C ¯uxes is assured by choosing the values of r1s and r2s such that the potential energy is suciently high at the set of points …r1s ; r2f ; R† or …r1f ; r2s ; R†. Parameters used in the calculations are listed in Table 2. Note that the number of grid points is rather large, making computations rather expensive. For instance, a large number of r1 grid points are needed to properly describe initially excited vibrational states (with several nodes in a small region of the r1 interval). In addition, all the three coordinates necessarily range up to quite high values in order to account for the di€erent asymptotic channels (see Table 1). The only exception is that, in order to reduce the computational cost, the negative limit of the R grid was chosen to  which is not an asymptotic value. be R ˆ 3:5 A, This restriction does not a€ect the results, since we have checked that the norm of the wave packet, restricted to the corresponding absorbing region, is rather small (below 0.3%) in all cases, except for H2 …v1 ˆ 14† ‡ H2 …v2 ˆ 0†, where it is about 3%, within our estimated global convergence error. After several tests, where the parameters of Table 1 were varied, the computed probabilities are considered to be well converged (within about 3%).

3. Results In Fig. 1, initial-state selected CID probabilities are presented and compared with previous quantum [1] and quasiclassical [8] calculations, for the initial states H2 …v1 ˆ 10; 12; 14† ‡ H2 …v2 ˆ 0†. Analogously, 4C reaction probabilities are shown in Fig. 2. The present results are referred in the Figures and in what follows as TD (time-dependent), those of [1] as TI (time-independent), and those of [8] as QCT. All calculations have been performed with the same model Hamiltonian and potential energy surface. The origin of the energy scale corresponds to H2 ‡ H2 , with the two molecules in®nitely separated from each other and in their equilibrium position (bottom of the reactants valley). When TI results were recovered, it was noticed that the origin actually used in [1] is 0.026 eV higher than the one used in the present as well as in the quasiclassical [23] studies. In the ®gures, this di€erence has been corrected by displacing the TI probabilities accordingly. Before discussing the comparison among TD, TI and QCT results, it should be said that the three methods give the same qualitative behavior for the energy dependence of the probabilities. First, dynamical thresholds for dissociation are larger than the static one (at 5.018 eV, which is equal to the dissociation energy of H2 , 4.748 eV, plus the H2 …v2 ˆ 0† vibrational energy, 0.270 eV), the e€ect being quite pronounced for v1 ˆ 10 but almost

Table 2 Parameters used in the time-dependent wave packet calculations for initial vibrational states (v1 ˆ 12, v2 ˆ 0) r1

r2

R

Number of grid points  Grid boundaries (A)

144 [0.05, 8.50]

144 [0.05, 8.50]

256 [)3.50, 8.50]

Absoption parameters:  qabs (A)  2) b (A

6.00 0.19

6.00 0.19

)2., 6. 0.19

Flux parameters:  rf (A)  rs (A)

5.00 1.50

5.00 1.50

Propagation time step Absorption time step

0.06 fs 1.50 fs

Parameters for other initial states can change slightly and are not shown here for simplicity. Note the use of two absorbing functions, for the negative and positive grid edges in the R coordinate, respectively.

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

Fig. 1. Collision-induced dissociation (CID) probabilities vs. total energy for state selected H2 …v1 † ‡ H2 …v2 ˆ 0† reactants. Present time-dependent wave packet calculations (TD, solid line) are compared with previous time-independent quantum calculations from [1] (TI, dashed line) as well as with quasiclassical trajectory calculations from [8] (QCT, crosses). (a) v1 ˆ 10; (b) v1 ˆ 12; (c) v1 ˆ 14.

absent for v1 ˆ 14. After threshold, dissociation probabilities rise more steeply and reach larger values as the initial vibrational state is higher, so that vibrational energy is more e€ective than translational energy for dissociation (vibrational enhancement). This behavior has been already observed in several studies of CID in non-reactive

181

Fig. 2. Same as Fig. 1 but for the four center (4C) reaction probabilities.

systems [24] (where `non-reactive' means, in a broad sense, that no reactive channel is e€ectively open in the corresponding energy range). After reaching a maximum value, dissociation probabilities decrease at roughly the same rate than that of increasing of the 4C probability after its threshold. This suggests an intimate competition between both processes. In addition, the 4C reaction probabilities also exhibit a vibrational enhancement e€ect as well as large dynamical thresholds.

182

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

It can be seen from Figs. 1 and 2 that the TD probabilities are more similar to the QCT than to the TI ones. TD and TI calculations di€er more notoriously at the CID threshold for v1 ˆ 10 (Fig. 1a), and in the di€erent rates of decrease (increase) of the CID (4C) probabilities for v1 ˆ 14 (Figs. 1c and 2c). It is somewhat surprising that two quantum approaches agree more poorly than one quantum and one classical treatment. The origin must be in some computational limitations of the TI calculations [1]. First, they were performed within only one of the (four) symmetry blocks of the system. To check if this is the origin of the discrepancies, we have performed time-dependent calculations within the same block of the TI study (using a properly symmetry-adapted initial wave packet), but obtained probabilities are almost identical to those of the non-symmetrized initial wave packet. Therefore the reason must be exclusively the convergence problems due to the use of approximate boundary conditions with hyperspherical coordinates [25], discussed in detail in [1]. Indeed, the uncertainty of the TI probabilities was estimated to be quite high, in the range 20±70%. The computations should have been extended to larger values of the hyperspherical coordinate or projections to Jacobi coordinates should have been performed, both alternatives being computationally quite expensive. The QCT probabilities compare much better with the new quantum-mechanical calculations (TD) than with the older ones: both TD and QCT probabilities agree quite well in the whole energy range studied. Although the energy grid used in the classical calculations is too sparse for a detailed analysis, we compare the dependence on energy of classical and quantum probabilities near thresholds, where major quantum e€ects are expected. For the dissociation process (Fig. 1), TD and QCT thresholds compare quite well for v1 ˆ 14 and well for v1 ˆ 12, whereas for v1 ˆ 10, the classical threshold seems to be a bit higher and sharper than the corresponding quantum one. A similar behavior has been noted in the CID of nonreactive systems [24,26], where it was found that, near thresholds, quantum probabilities are smoother than classical ones (`quantum tails'), the e€ect being more pronounced for less excited ini-

tial vibrational states. This type of e€ect is also found, and in a more pronounced manner, in the 4C reaction probabilities: quantum thresholds are smoother and compare worse with the classical thresholds as the initial state is less excited. To gain further insight into the extent of quantum e€ects, we have computed the product of the 4C or CID quantum probabilities and the translational Maxwell±Boltzmann distribution (such a product would take part in the computation of initial state selected rate constants), Fv1 ;v2 …E; T † ˆ Pv1 ;v2 …E† exp‰ …E

ev1 ;v2 †=kB T Š;

…8†

where ev1 ;v2 is the vibrational energy of reactants, and studied the dependence of F on energy and temperature, T. This function is shown in Fig. 3 for the case of the 4C reaction with …v1 ˆ 10; v2 ˆ 0† and a temperature of 1200 K. The dependence on the energy of TD quantum and QCT probabilities is also displayed in the ®gure. The distribution peaks near 5.6 eV. At this energy, the calculated classical probability is zero or almost zero but the quantum probability is, however, non-negligible. At a rough guess, the area under a classical distribution, P QCT …E† exp‰ …E e†=kB T Š (or integral over E), would be less than half of that of the quantum distribution. In this rough way, we estimate that quantum e€ects are important below 1200 K. In the same manner, we have estimated

Fig. 3. Dependence on the total energy of Fv1 ;v2 …E; T † (see Eq. 8), for the case of the present (TD) quantum four center reaction probability with …v1 ˆ 10; v2 ˆ 0† and a temperature of 1200 K. Reaction probabilities themselves (quantum and classical [8]) are also displayed, but referred to the right-hand side vertical axis. See text for more details and discussion.

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

that below 800 and 500 K, quantum e€ects are important for the 4C reaction with v1 ˆ 12 and 14 initial vibrational states, respectively. For the CID process, quantum behavior seems to dominate below 500 K for v1 ˆ 10, whereas for v1 ˆ 12, 14, the temperatures for which quantum e€ects would play a role in the rate coecients would be rather low, below 100 K. These predictions are more grounded for the 4C than for the CID process, since full-dimensional quasiclassical calculations [8] suggest that the present model is more realistic for the reactive process. Before ending, we would like to note the double-peaked structure (just after threshold) of the v1 ˆ 14 TD CID probability (Fig. 1c). It has been found that such a structure is due to the competition between two mechanisms: in one, the pairs of diatoms recoil after their closest approach, while in the other, the bounded diatom inserts into the nearly dissociated one [27]. It would be interesting to know if a similar e€ect is also found classically.

183

dissociation) specially for the 4C reaction probabilities and the lowest initial vibrational states. Temperatures where these e€ects would be manifest have been estimated. However, to draw a de®nitive conclusion, it would be convenient to perform further classical calculations on a ®ner energy grid and perhaps increasing the number of trajectories. A study of isotopic substitution effects, both classical and quantum-mechanically, would also contribute to the understanding of the mechanisms and the quantum e€ects of the H2 ‡ H2 reaction and dissociation. Acknowledgements We thank Prof. E. Garcõa for sending us ®les with the quasiclassical trajectory results, and A. Ceballos and M.E. Mandy for discussions and for sending their works prior to publication. This work has been supported by the DGICYT (Spain) grant number PB95-0071 and EU contract HPRNCT-1999-0005.

4. Conclusion We have presented a time-dependent wave packet approach to compute 4C reaction and collision-induced dissociation probabilities for H2 ‡ H2 in selected initial vibrational states and within a reduced dimensionality model. It has been found that the method provides more accurate results than a time-independent hyperspherical treatment previously applied to the same problem [1]. The present results compare quite well with quasiclassical trajectory calculations [8] over a broad range of collision energies, which is somewhat surprising in view of the light atoms involved in the reactive process. This result, within a reduced dimensionality model, provides a ground for full-dimensional classical calculations (some already reported in the literature [8,11]), since fulldimensional exact quantum calculations are very dicult to carry out at present. A comparison of quantum and classical behavior near reaction thresholds suggests, however, the presence of tunneling e€ects (or wall penetration, in the case of

References [1] M.I. Hernandez, D.C. Clary, J. Chem. Phys. 104 (1996) 8413. [2] A. Aguado, C. Suarez, M. Paniagua, J. Chem. Phys. 101 (1994) 4004. [3] J.Z.H. Zhang, J. Dai, W. Zhu, J. Phys. Chem. A 97 (1997) 2746. [4] D.C. Clary, J. Phys. Chem. 98 (1994) 10678. [5] S.H. Bauer, Ann. Rev. Phys. Chem. 30 (1979) 271. [6] K. Sakimoto, J. Chem. Phys. 112 (2000) 5044, and references therein. [7] R.T Pack, R.B. Walker, B.K. Kendrick, J. Chem. Phys. 109 (1998) 6701, and references therein. [8] A. Ceballos, E. Garcõa, A. Rodrõguez, A. Lagana, Chem. Phys. Lett. 305 (1999) 276. [9] A. Ceballos, E. Garcõa, A. Rodrõguez, A. Lagana, J. Phys. Chem. A 105 (2001) 1797. [10] P.G. Martin, W.J. Keogh, M.E. Mandy, Astrophys. J. 499 (1998) 793. [11] M.E. Mandy, P.G. Martin, W.J. Keogh, J. Chem. Phys. 108 (1998) 492. [12] D. Zhang, J.Z.H. Zhang, J. Chem. Phys. 101 (1994) 1146. [13] D.W. Schwenke, J. Chem. Phys. 89 (1988) 2076. [14] C. Park, J. Thermophys. Heat Transfer 7 (1993) 385. [15] A. Lifshitz, M. Bidani, H.F. Carroll, J. Chem. Phys. 79 (1983) 2742.

184

D. di Domenico et al. / Chemical Physics Letters 342 (2001) 177±184

[16] M.E. Mandy, T.A. Rothwell, P.G. Martin, J. Chem. Phys., in press. [17] A.I. Boothroyd, J.E. Dove, W.J. Keogh, P.G. Martin, M.R. Peterson, J. Chem. Phys. 95 (1991) 4331. [18] E.J. Heller, J. Chem. Phys. 62 (1975) 1544. [19] M.D. Feit, J.A. Fleck Jr., A. Steiger, J. Comput. Phys. 47 (1982) 412. [20] R. Heather, H. Metiu, J. Chem. Phys. 86 (1987) 5009. [21] P. Pernot, W.A. Lester Jr., Intern. J. Quantum Chem. 40 (1991) 577.

[22] W.H. Miller, J. Chem. Phys. 61 (1974) 1823. [23] A. Ceballos, private communication. [24] See, for instance, C. Leforestier, in: D.C. Clary (Ed.), The Theory of Chemical Reaction Dynamics, Reidel, Dordrecht, 1986, p. 235 and references therein. [25] J. R omelt, Chem. Phys. Lett. 74 (1980) 263. [26] K. Nobusada, K. Sakimoto, Chem. Phys. Lett. 288 (1998) 311. [27] D. di Domenico, M.I. Hernandez, J. Campos Martõnez, in preparation.