I. Introduction

2

2

2

2

The above reaction belongs to the family of hydrogen exchange reactions. In contrast to the neutral reaction H + H₂ H₂ + H it has a shallow potential well in the entrance channel, but the reaction barrier in the interaction region is of comparable magnitude. Compared to the positively charged system H⁺ + H₂ H₂ + H⁺ with a deep minimum of 4.61 eV for the equilibrium structure of H₃⁺, the well depth of 0.05 eV is relatively small. What distinguishes ionic systems from neutral systems is that because of the long-range inductive interaction potential (V ~ -1/R⁴) the important range of the PES is much more extended (at least for our systems: R (atom-diatom)> 16 a₀) than in case of neutral systems.

In recent years, H₃^- has been investigated rather intensively experimentally and theoretically. There exist two types of calculations for the potential energy surface: (a) the ab initio calculations of the ground-state PES by Stärck and Meyer (SM)⁶ and (b) the diatomics-in-molecules (DIM) investigations of Belyaev et al.^7-9 in which several electronic states were included. The first quantum mechanical calculation for the collinear H₃^- system⁷ was performed with the DIM potential using the S-matrix version of the Hulthén-Kohn variational method of Miller and co-workers.^10-12 The two different PESs have been used by Gianturco and Kumar^13,14 investigating vibrational and rotational inelasticity and by Mahapatra,¹⁵ Mahapatra and Sathyamurty¹⁶ concentrating on the reactive dynamics using time-dependent approaches. The most recent experimental investigations are the guided ion beam experiments of Haufler et al.¹⁷ and the crossed beam experiments of Zimmer and Linder.¹⁸ Both groups report integral and differential cross sections for different isotope variants. At total energies below the dissociation limit of the hydrogen molecule, the results of the H^- + H₂ collison can be inelastic excitation, rearrangement, and electron detachment without or including rearrangement. From the analysis of the potential energy surfaces, one knows that the electron detachment channel opens up at an energy of 1.2 eV.

This paper is organized as follows: in section II we present the theoretical approach and in section III we will concentrate on our applications.

II. Quantum Dynamics

J

Figure 1 Sketch of reactant, product and interaction region for the H- + H2 reaction. To represent the different arrangement channels, the symbols for the three different hydrogen isotopes have been used. Two different approaches are shown: (a) the reaction T- + HD starts in reactant Jacobi coordinates (RC), (b) the reaction T- + HD starts in the product Jacobi coordinates (PC) of the TH + D- arrangement.

B. Formalism. The general theory and various numerical aspects of time-dependent reactive scattering have been discussed elsewhere.^19,33-37 Nonetheless, we want to summarize the main features of our own wave packet implementation. We start with the initial wave packet (see Figure 2a)

where R is the distance between the atom A and the center of mass of the BC molecule, {q_i} are the internal degrees of freedom for the BC molecule, and (q) is the initial state of the BC molecule for given vibrational (v) and rotational (j) quantum number. R₀, k₀, and  define the initial location of the center of the wave packet in coordinate and momentum space and the initial width.

Figure 2 Sketch of the interaction region for the H- + H2 reaction: (a) reactant and (b) product Jacobi coordinates with a given initial wave packet. (To represent the different arrangement channels, the symbols for the three different hydrogen isotopes have been used.)

The dynamics of the system is followed by solving the time-dependent Schrödinger equation for nuclear motion numerically

where  is used as an index for the different arrangements ( = A + BC,  = AB + C,  = AC + B), and the formal solution is given by

with Û being the evolution operator for the Schrödinger equation and where (t = 0) and (t) are the wave function of the system at time 0 and t, respectively.

The Hamiltonian operator in Jacobi coordinates (R, r, ) for the body-fixed frame ( = 1) is given as

with

Here, _R is the reduced mass of the A - BC system, _r is the reduced mass of the BC (diatom) molecule, J is the total angular momentum of the system, j is the rotational angular momentum for BC, V(R, r, ) is the potential energy surface and J_±, j_± are usual ladder operators.

The wave function for a particular total angular momentum quantum number J and its space fixed z components can be expressed in terms of body-fixed coordinates in the form

where M is the quantum number of the projection of total angular momentum J on the space-fixed z-axis and  is the quantum number for the projection of j on the body-fixed z-axis.  denotes the Wigner rotation matrix elements. The resulting Hamiltonian in the body-fixed frame is given in a tridiagonal matrix representation, where, according to eq 5, the diagonal part is of the form

with the coupling term

The action of terms of type H and H_,± on the wave function can be computed independently and is optimal for a parallel implementation.²⁰

The propagation of the wave packet can be performed in two different ways. In the first case we use the Chebyshev propagation method as proposed by Tal-Ezer and Kosloff.²¹ The evolution operator is developed in complex Chebyshev polynomials _n

The Chebyshev polynomials _n are calculated with a recursion formula

where x = -iH_norm, ₀ = (0), ₁ = x(0).

Alternatively, one can perform one long propagation step as described by Mandelshtam and Taylor (see eq 3.1 in ref 23):

The absorbing potential is simulated by (R,r):

In eqs 16 and 17, H_s is a scaled and shifted Hamiltonian operator

where a_s = 2/E, b_s = 1 - a_sE_max as defined in the equivalent recursion formula by Gray and Balint-Kurti (see eq 12 in ref 30), which can be used for the complex wave packet (as we do it in the present paper) or its real part only.³⁰ The recursion relation can be used to create the complex wave function for a given time (information from eq 16 has to be inserted into the e^(-it/) expansion) or is used as a dynamical iteration alone as in the Gray-Balint-Kurti-approach. In our implementation we use the various schemes as described above. For the final analysis we store, in the case of one long propagation, information from each iteration step; within a test-run one can clarify if for each iteration the information has to be stored. Details are given in the work of Gray, Balint-Kurti, and coworkers.^30,31

The second method of propagation of the wave packet is the split operator technique where we approximate the propagator by dividing the total time into N segments t using the Trotter formalism of second order. Each short time propagator e^-iHt/ is again splitted symmetrically in the following way (potential referenced expression)

To get an overview about other ways of performing the action of the propagator the reader is referred to the literature.^40,41

In both schemes (CH, SO) the action of the radial kinetic energy on the wave function is used in conjunction with the fast Fourier method (FFT); in case of angular kinetic energies a discrete variable representation method is used (DVR). The kinetic and potential energy is given as follows. Radial kinetic energy for the R coordinate (and for r, respectively) is (a) for Chebyshev:

where k is the wave vector of a given plane wave. For (b) for SO:

The angular kinetic energy (angle  coordinate) for (a) Chebyshev is

with

and for (b) SO:

with (1/I) = ((1/_RR²) + (1/_rr²)).

Although the last two formulations (eqs 23 and 25) seem to be clear with respect to numerical implementation, there have been some complications if not enough memory is available within the program code. In the case of Chebyshev implementation, the moment of inertia I_il is indexed with respect to R_i and r_l values. If one defines a maximum energy value for the angular kinetic energy contribution, it might be possible to define a j_max = f(E_max) value for the case that the prefactor (1/2I) does not become too large. In that case, one can save computation time by storing a new matrix U_kk'

in which (1/2I) has been separated out. In case of large (1/2I_il) values (i.e., for R 0) several of these matrices have to be stored. We have in our program defined an E_max value which leads to individual j_max values for given i and l; we think that this is more accurate in the numerical implementation than a global definition of j_max. In the case of SO calculations, this trick of separation of (1/2I) is not possible. If enough core memory is available, one can store

otherwise the transformation given in eq 25 becomes the most time-consuming step during the propagation, although H is not as often performed as in the case of Chebyshev propagation.

The potential energy is diagonal in the coordinate grid representation, and the action in the CH approach is just a multiplication

whereas in the SO approach, the wave function has to be multiplied twice with the e^{-iV(R_ir_l,_k)t/2} factor. In the moment, we use the SO approach only in case of diagonal potential energy representation.

As mentioned before, the wave function is discretized on a 3D grid in Jacobi coordinates R, r, and  with typically N_R = 64-128, N_r = 64-128, and N = 32-80 grid points. Ionic systems require a finer discretization, because deep potentials (typically V_min -0.5  -4.5 eV) lead to strong resonances. The propagation time is set (a) by using N segments of t in the split operator formalism, (b) by setting a maximum propagation time and by that fixing the number of N Chebychev iterations for the propagation as proposed by Kosloff, or (c) by doing one long propagation step as proposed by Mandelshtam and Taylor and fixing for that the number of Chebychev iterations. Of course, after the propagation is performed for a given number of steps, the run can be restarted if the reaction probabilities do not seem to be converged.

C. Analysis. In our calculations we can perform energy-resolved state-to-state and flux calculations. The wave packet is propagated until it has "completely" left the interaction region. To perform the analysis we need the wave packet in the correct Jacobi coordinates, depending on the different possible arrangement channels . A basic difficulty in the theory of reactive collisions is that the coordinates appropriate for reactant and product arrangements differ from each other. There have been different ways presented in the literature to solve the problem:^42-45 one natural way would be to transform the wave function to the appropriate coordinates. This approach is time consuming and leads to numerical errors.

Because of that reason, we perform for each reactive product arrangement channel an individual scattering calculation performed in the appropriate product Jacobi coordinates as illustrated in Figures 1 and 2. If only global reaction probabilities are needed, one can calculate within the reactant Jacobi coordinates the state-to-state inelastic transition probabilities and calculate the reaction probabilities as the difference from unity:

The disadvantage of propagating in product coordinates is that the representation of the starting wave packet needs more angular grid points for an appropriate description with a given quality of the norm of the wave function. The analysis for the reaction probability or the S-matrix can be performed in two ways: (a) we can split off the wave packet in the asymptotic region (as proposed by Heather and Metiu²⁴) or (b) calculate the contribution at an asymptotic analysis line as proposed by Balint-Kurti et al.²⁷ In case of splitting off the wave function,  is divided in two parts, the product part _P and the interaction part _I as given in eqs 30 and 31:

R_p and C define the shape of the split-off function. After each time step, the product part is split off and the interaction part is analyzed. Its contribution to the S-matrix is given by

After each splitting, the propagation is continued, with (1 - f) replacing . This procedure has the quality of an absorbing potential so that the wave packet does not reach the end of the grid point area; otherwise, this leads to numerical problems in the FFT approach. In case of split-off functions, no further absorbing potential is needed in the calculation.

In the Balint-Kurti-approach,²⁷ the S-matrix is calculated by the following procedure: in the asymptotic region, i.e., R is large, first the wave packet along a cut (R_ana) is projected onto the final product state _F(r) to produce a set of time-dependent coefficients C_F,I(t) (e.g., using product Jacobi coordinates):

Then these coefficients are Fourier transformed over the time to give energy dependent coefficients A_F,I(E)

and the S-matrix results from the simple relation (derivation given in refs 38 and 39):

If just energy resolved reaction probabilities (E) are desired, (E) is computed from (e.g., using reactant Jacobi coordinates)

where in case of reactant Jacobi coordinates the quantity in the brackets is the energy-resolved flux of the wave packet at the asymptotic dividing surfaces defined at the position r_ana (for reactive analysis) or at R_ana (for inelastic analysis); the angular brackets denote integration over the other two coordinates. The energy-dependent wave function (R, r, , E) is obtained by Fourier transforming (R, r, , t). The described formulation of the Balint-Kurti analysis approach is straightforward if using the wave function as given in eq 11 using the recursion relation in eq 15. In case of using the recursion formula given in eq 16, one has to take into account that the energy scale was shifted and scaled (see ref 30; for flux analysis: see ref 31).

To avoid unphysical reflections of the wave function arising from finite size boundary of the coordinate grid (this is needed at least for the asymptotic analysis procedures at asymptotic cut-lines of flux calculations through a dividing surface), the outer part of the radial part of the grid is surrounded by an optical (absorbing) potential of the form proposed by Vibok and Balint-Kurti;³² we mostly use the type given in eq 18.

To summarize our approach: there are two ways for a state selective analysis: (1) Figure 2a: For inelastic 3D-investigations we use reactant Jacobi coordinates, so that a state-to-state inelastic analysis is possible and energy resolved total reaction probabilities can be calculated. (2) Figure 2b: In case of reactive 3D-investigations we use product Jacobi coordinates, so that state-to-state reaction probabilities can be calculated.

The advantage is that no coordinate transformation of the wave packet is necessary. The transformation of the wave packet, from one coordinate system (e.g., reactant coordinates) to the one that is optimal for asymptotic analysis (e.g., product coordinates) is time-consuming and leads to numerical inaccuracies.

For most of the calculations in this present work, the width of the wave packet is chosen such that within one WP calculation accurate results for a collision energy range of ± 1 to 2 eV around the initial starting collision energy can be achieved. The quality of the propagation depends, for a given collision energy, on the propagation time length and all internal parameters, such as grid size for the coordinates, width of the wave packet, etc. We accept a propagation as being reasonable if the contributions coming from the r_ana and R_ana analysis line add up to 100 ± 1-3%. This test has been performed by adding up (a) fluxes through the r_ana and R_ana intersection line, (b) adding up flux through the r_ana and state-to-state contributions calculated at the R_ana intersection line (or vice versa), or (c) adding up from two individual runs the state-to-state contributions calculated at the R_ana intersection line for reactant and product Jacobi coordinates.

III. Results and Discussion

Reaction a will be discussed in the present paper, reactions b and c in a forthcoming paper.⁵¹ The H₃^- system is of great importance in hydrogen plasmas.¹⁸ We are using the potential energy surface of Stärck and Meyer⁶ and the diatomics-in-molecules (DIM) potential of Belyaev.⁷ A comparison of scattering calculations with the experiment will be a test for the quality of DIM surfaces, which nowadays are often used for reactive dynamics, especially on coupled surfaces. In case of inelastic investigations for H₃^- at high energies (up to keV), the DIM results had been reasonable compared to experiment.¹⁴

The 3D DIM PES for the ground state and some low lying excited states was reported by Belyaev et al.^8,9 Ground and first excited state form a conical intersection at D_3h symmetry. The ab initio surface of Stärck and Meyer⁶ has been determined by MR-CI and CEPA(2) calculations for 403 nuclear configurations. It has been cast in an analytical form by a fit with 23 parameters and an rms error of 0.23 mE_h. The PES fit presented in ref 6 includes some misprints, so we were happy that W. Meyer could provide us with his code for the PES. This code had been tested by using quasiclassical trajectory calculations,⁵² and no numerical problems from the fit were encountered. Within our WP calculations we found that in a small region in the strong repulsive part of the potential, i.e., at short H-H-distances, the energy tends to go to minus infinity (this results from eq 1 in ref 6). We corrected this behavior, so that the only minimum region on the global surface is the linear van der Waals minimum in the entrance channel (r = 1.416a₀, R = 6.183a₀,  = -0.0476 eV). The position and height of the barrier have been found at r_bar = 1.997a₀ and  = 0.454 eV. The electron detachment seam has been determined and its lowest point was found for R = 2.86a₀, r = 1.42a₀, and  = 90 with an energy of 1.2 eV, which agrees experimental findings.¹⁸ The barrier height of the DIM potential (at the position r_HH = 1.74a₀) is by 0.17 eV ( = 0.624 eV) larger than the one of the SM potential. The DIM potential also exhibits a shallow van der Waals minimum in the entrance channel ( -0.05 eV) and has a energy threshold value of 1.46 eV for electron detachment.

Because of the existence of several electronic excited states at energies E > 1.2 eV, a nuclear dynamics calculation involving several electronic states is necessary for E > 1.2 eV. In this paper we will present calculations on a single surface, which should be reasonable at least up to E = 1.2 eV. The influence of higher excited states on the dynamics has to be investigated in the future. Collinear reactions have been investigated in recent years using both the DIM and the SM potential. Inelastic scattering calculations for higher collision energies has been reported by Gianturco and Kumar.¹³ The first energy resolved reaction probabilities in three dimensions on the DIM PES using the time-dependent WP approach has been presented by Mahapatra;¹⁵ some preliminary calculations using the SM PES have been reported by Mahapatra and Sathymurthy.¹⁶

In this paper we will concentrate mostly on time-dependent WP calculations using the SM potential, because, as will be seen in the next subsections, the DIM PES leads to very different results concerning the reaction probabilities. To check our own code and to compare with time-independent approaches, we performed in addition calculations using the hyperspherical coordinate method of Manolopoulos et al.⁵⁴

The numerical grid parameters and properties of the initial wave function used in the calculations of total and state-to-state reaction probabilities are summarized in Table 1.

B. H^- + H₂. In the following three figures (Figures 3, 4, and 5), we present results of energy-resolved total reaction probabilities of H^- + H₂(v,j) for different intial vibrational (v) and rotational (j) states. In Figures 3 and 5 we show in addition the results for the neutral reaction H+H₂ using the LSTH potential⁴⁷ and compare the two potentials SM and DIM. For the SM potential, the reaction starts at a total energy of E  0.55 eV (E_trans = 0.46 eV), similar to the neutral reaction H+H₂. Because the barrier height is larger in case of the DIM potential, there the reaction starts 0.17 eV later. The forms of the onset of the reaction probability are completely different if one compares SM and DIM. In case of the SM potential, the reaction probability for initial H^- +H₂ (v = 0, j = 0) (Figure 3) increases steeply up to 70-80% and falls down slowly around 1.8 eV, whereas for DIM the reaction probability reaches its maximum for E > 2 eV. For the SM potential, the reaction probability is a relatively smooth curve comparable with the neutral H+H₂ reaction, except that resonances occur in some smaller energy ranges for E > 1.2 eV. If one compares Figures 3 (with initial H₂ (v = 0, j = 0)) and 4 (with initial H₂ (v = 1, j = 0)), one will see that for H^- + H₂ (v = 1, j = 0) the resonance structure is much more pronounced between E = 1.2 to 1.3 and E = 1.6 to 1.8 eV. Figure 5 includes results for the different starting conditions for H₂ in different vibrational (v = 0, 1,...) and rotational states (j = 0, 1,...). The most impressive feature is that for H₂ (v = 0, j = 1) the reaction probability increases steeply up to 95% with some small resonance features around E 1.3 eV. The influence of initial rotational excitation has been experimentally investigated only in the case of inelastic scattering.

	Figure 3 H- + H2 (v = 0, j = 0): total reaction probabilities calculated for the two different potential energy surfaces (SM and DIM) and with two Jacobi coordinate systems (RC and PC). The analysis is made either from fluxes with the SO propagator (SO) or from summed state-to-state reaction probabilities (S-F: state specific in the product region and flux in the reactant region) using the Chebychev propagator. For comparison, the result for the neutral H + H2 reaction (LSTH potential47) is given. Etrans = 0.7 eV; Rmax = 15.5(12.0)a0 for H3- (H3).
	Figure 4 H- + H2 (v = 1, j = 0): total reaction probabilities calculated for the two different potential energy surfaces (SM and DIM) and with two Jacobi coordinate systems (RC and PC). The analysis is made either from fluxes with the SO propagator (SO) or from summed state-to-state reaction probabilities (S-F: state specific in the product region and flux in the reactant region) using the Chebychev propagator. (Rmax = 15.5a0).
	Figure 5 H- + H2 (v = 0,1, j = 0,1,2)(SM potential): total reaction probabilities calculated for the SM potential energy surface and for different starting conditions of the diatomic H2 using RCs and PCs. The analysis is made from summed state-to-state reaction probabilities (S-F). For comparison, the result for the neutral H + H2 reaction (LSTH potential) is given. (Rmax = 15.5a0).

Within Figures 3-5 we compare different wave packet approaches (different coordinates, propagators, etc.). In principle the results are qualitatively the same (in order not to overload the figures, not every result is shown), but if one is interested in very accurate, detailed information, the Chebychev method is clearly the one to be preferred. To get a rough overview of the reaction probability, the split-operator approach performs better, compared to the Chebychev method (this results from fewer actions of the Hamiltonian within the propagation). For the H₃^- system we find only few resonances, which are smoothed a bit when using the split-operator approach (presumably the time step t has to be shorter than 10 [au]). This smoothing behavior is especially a result that we have found within our investigations of the Ne + H₂⁺ reaction,⁵⁵ where we have a much deeper potential minimum leading to many sharp resonances. The comparison of reactant (RC) and product coordinates (PC) reveals that if one is interested just in total reaction probabilities for systems such as H^- + H₂, calculations with RC are good enough. The strong difference in reaction probabilities for H₂ in different rotational states (v = 0, j = 0, 1, 2,...) reveals the strong influence of the anisotropy of the PES on the dynamics. For initial vibrational states v = 2, 3, 4, further total reaction probabilities are given in Figures 6 and 7.

	Figure 6 H- + H2 (v = 2, j = 0)(SM potential): total reaction probabilities using two Jacobi coordinates (RC and PC).
	Figure 7 H- + H2 (v = 3,4, j = 0)(SM potential): total reaction probabilities using two Jacobi coordinates (RC and PC).

The work of Mahapatra and Sathyamurthy¹⁶ performed with the SM potential shows much more pronounced resonances, which we attach especially to a too early starting position for the wave packet, the use of the SO method, and a smaller number of grid points. In the work of Mahapatra¹⁵ using the DIM potential, most of these ingredients have been improved, but the DIM potential is less demanding with respect to resonances, except for H^- + H₂ (v = 2), where some resonance features are seen in the onset of the reaction probability. On the other hand, we see strong resonance features for all initial H₂ (v > 0), but one has to keep in mind that most of these features do happen at energies above the opening of the detachment channel.

In Figure 8 (a smaller energy section is seen in Figure 9), we test the influence of increasing the range of the interaction potential. The maximum value for R has been chosen between 15.5 and 19.5a₀. Results presented in Figures 3-5 have been calculated for R_max = 15.5a₀, and the new results presented in Figure 8 reveal that R_max was too short, because the reaction probability around E = 0.9 eV is too low. In addition, we compare (a) our code for complex wave packets, (b) the real wave packet code of Stephen Gray,⁵³ and (c) the time-independent "ABC"-code of Manolopoulos et al.;⁵⁴ for the three different codes all calculations have been performed by ourselves. The results for the complex and real wave packets are nearly identical if the same numerical parameters have been used (Figure 9). There is a slight difference in the results when we compare WPs with different initial starting collision energies (Figure 9: E_trans = 1 eV, 2 eV). If the hyperspherical radius  (rhomax) in the time-independent calculations (abc) is chosen large enough (i.e., compare rhomax = 12 and 15 a₀), then the comparison with "S-F: PC, R_max = 19.5, E_trans = 1 eV" is nearly perfect. The resonance structure in the time-independent calculation (abc) is a little bit different from the WP results, but this comes from the less dense energy grid in case of the abc calculation. In addition, one can see that with initial E_trans = 2 eV we can get reasonable results in one run up to the dissociation energy of H₂. In case of initial E_trans = 1 eV, values only up to E = 3.5 eV are acceptable. At the moment, it is not clear why we have a "constant" shift between the results calculated with E_trans = 1 and 2 eV for R_max = 17.5a₀.

Figure 8 H- + H2 (v = 0, j = 0)(SM potential) energy-range: 0.5 5 eV, total reaction probabilities calculated in PCs from summed state-to-state reaction probabilities (S-F). Different grid sizes and initial collision energies have been tested: Rmax = 15.5, 17.5 and 19.5 a0, Etrans = 1.0 and 2.0 eV. Comparison with a real wave packet (real WP) calculation and with a time-independent hyperspherical approach (abc: time-ind.) is presented. In the time-independent calculation, two different maximum values for the hyperspherical radius max have been tested.

Figure 9 H- + H2 (v = 0, j = 0)(SM potential): as Figure 8, but energy range: 0.7  2 eV.

In the next figures (Figures 10, 11, 12, 13), we present results for state-to-state reaction probabilities. If one starts initially with v = 0, j = 0, then for the product molecule mainly vibrational state v = 0 is preferred (Figure 10), whereas in case of the final product H₂ (v = 0), especially the rotational states j = 1 and j = 0 (Figure 11) dominate; similar results are seen for for final vibrational v = 1 (Figure 12) and v = 2 states (Figure 13).

	Figure 10 H- + H2 (v = 0, j = 0)  H2 (v', j') + H- (SM potential): state-to-state reaction probabilities (in PCs) for different rotationally summed vibrational product states (Rmax = 17.5 a0, Etrans = 1.0 eV).
	Figure 11 H- + H2 (v = 0, j = 0)  H2 (v' = 0, j') + H- (SM potential): state-to-state reaction probabilities (in PCs) for different rotational product states (with v' = 0) (Rmax = 17.5 a0, Etrans = 1.0 eV).
	Figure 12 H- + H2 (v = 0, j = 0)  H2 (v' = 1, j') + H- (SM potential): state-to-state reaction probabilities (in PCs) for different rotational product states (with v' =1) (Rmax = 17.5 a0, Etrans = 1.0 eV).
	Figure 13 H- + H2 (v = 0, j = 0)  H2 (v' = 2, j') + H- (SM potential): state-to-state reaction probabilities (in PCs) for different rotational product states (with v' = 2) (Rmax = 17.5 a0, Etrans = 1.0 eV).

In Figure 14, we present results for the total reactive cross section using the J-shifting method of Bowman et al.^48-50 Time-dependent and time-independent calculations lead to similar results. The onset of the cross-section (up to E = 1.2 eV) is the same, regardless of how many J-states are included. J-shifting does not take into account that the reaction probability decreases considerably at higher total angular momentum. Our newest unpublished results show that, up to J = 8, the magnitude of the reaction probabilty does not change significantly and that from then on the reaction probability consistently decreases, reaching for J = 20 (including Coriolis coupling) a value that is already a factor of 2-3 smaller than the one for J = 0 with a maximum around E = 2-2.5 eV. This is the reason why in Figure 14 mostly results for the J = 0  20 summation are presented. In a forthcoming paper we will present results where integral and differential cross sections are derived explicitly from calculations for total angular momentum J 0.

Figure 14 H- + H2 (v = 0, j = 0)(SM potential): total cross sections calculated using the J-shifting method (JS). A comparison of wave packet and time-independent results (abc) is performed. The contribution of total angular momentum J was fixed to 20 and 50.

IV. Summary

Part of our project is to clarify, which numerical technique is the most efficient for solving wave packet propagation, especially when using the code on multiprocessor machines. For the present applications, we can say that the Chebychev method is clearly the one to be preferred, because the accuracy can be improved by propagating in time as long it is needed. The recursion formula in eq 16 is independent of the magnitude of the time step. If high accuracy is preferred, the split-operator method needs time-steps of t < 10 [au] and is not reliable for long propagations, as are needed in the case of deep potentials. In addition there are problems in handling eq 25 optimally for R 0. The efficiency of the wave packet approach compared to the time-independent hyperspherical approach of Manolopoulos et al.⁵⁴ has not been analyzed in detail, but in case a fine energy grid for the reaction probability was needed, the time-independent approach was more time-consuming (but on the other hand the full S-matrix was calculated).

The advantages of wave packet (WP) calculations are that WPs are relatively easy to apply for systems with many reaction channels, that WPs are suitable for the use on several PESs, and that with WPs one gets, within one calculation, information about a large collision energy range (E_coll > 1 eV). The present disadvantage is that, in the case of deep potential energy minima, WP calculations can become very CPU intensive, because the wave packet is trapped. Applications of wave packet calculations are performed now on an IBM-SP2 using parallel architecture. Work of optimizing the "parallel" code is in progress and will be published soon.

In the case of the H₃^- system, we presented energy-resolved and state-to-state reaction probabilities for different starting conditions. It is for the first time that the nuclear dynamics within the H₃^- system has been investigated in such a detail using mostly the ab initio potential of Stärck and Meyer (SM)⁶ and that the results have been checked by us with two other time-dependent and time-independent approaches. We calculated reaction probabilities for the DIM and ab initio PES: the results are qualitatively different.

We could show that our results are quantitatively comparable with the real wave packet code of S. Gray⁵³ and with the "ABC" code (time-independent hyperspherical coordinate method) of Manolopoulos et al.⁵⁴ A direct comparison between theory and experiment can be made when our integral and differential cross sections are available. We are in the process of finishing this work. Most of the experimental work for reactive analysis is performed for the isotopic variants HD₂^- and DH₂^-; our results for these systems will be soon available. We are planning now to perform the WP calculations on several surfaces, initially using the DIM approach.

Acknowledgment