State-to-state reaction probabilities within the quantum transition state framework

Rigorous quantum dynamics calculations of reaction rates and initial state-selected reaction probabilities of polyatomic reactions can be efficiently performed within the quantum transition state concept employing flux correlation functions and wave packet propagation utilizing the multi-configurational time-dependent Hartree approach. Here, analytical formulas and a numerical scheme extending this approach to the calculation of state-to-state reaction probabilities are presented. The formulas derived facilitate the use of three different dividing surfaces: two dividing surfaces located in the product and reactant asymptotic region facilitate full state resolution while a third dividing surface placed in the transition state region can be used to define an additional flux operator. The eigenstates of the corresponding thermal flux operator then correspond to vibrational states of the activated complex. Transforming these states to reactant and product coordinates and propagating them into the respective asymptotic region, the full scattering matrix can be obtained. To illustrate the new approach, test calculations study the D + H(2)(ν, j) → HD(ν', j') + H reaction for J = 0.


I. INTRODUCTION
A detailed understanding of chemical reaction processes is one of the central challenges of chemical physics and theoretical chemistry.The most detailed results, fully quantum state resolved reaction probabilities and cross sections, are available for reactions involving only few atoms.Within the last decade, significant progress towards the investigation of polyatomic reactions has been made and experiments studying state-to-state reactive scattering of six atom reactions as F + CH 4 → HF + CH 3 (Refs.1-3) or Cl + CH 4 → HCl + CH 3 (Refs.4-7) have been published.][10][11][12] Here, almost all work focused on the H 2 + OH → H + H 2 O reaction [8][9][10][11] and its isotopic analogs 13 and only very recently Zhang and co-workers 12 presented fully state-resolved calculations for the HO + CO → H + CO 2 reaction.
The aim of the present work is to develop an approach in order to also calculate state-to-state reaction probabilities employing the quantum transition state concept and MCTDH wave packet propagation.This approach could facilitate future calculations of completely reactant and product state resolved observables for six atom reactions as X + CH 4 → XH + CH 3 .(Developing a transition state based description of state-to-state reaction dynamics, the state-tostate-to-state model of Gustafsson and Skodje 44 should be mentioned.While their model is, in contrast to the present work, approximate in nature, it still draws on the same underlying physical pictures and utilizes the vibrational states of the activated complex to analyze the state-selectivity of the reaction process.) Extending the quantum transition state concept towards the description of state-to-state reaction dynamics, two issues not present in calculations of cumulative or initial stateselected reaction probabilities have to be addressed.
First, standard flux correlation functions [45][46][47] use only up to two different dividing surfaces.However, calculating stateto-state reaction probabilities by propagation of thermal flux eigenstates requires three different surfaces: one dividing surface located in the transition state region is used to define the initial wave packets and two other ones located in the reactant and product asymptotic area are employed for the analysis of the quantum states of the reactants and products, respectively.Thus, the existing flux correlation functions are not sufficient for the present purpose and new analytic formulas must be derived.The necessary developments will be presented in Sec.II.
Second, the simultaneous analysis with respect to quantum states of reactants and products requires the use of at least two different coordinate systems.Wave packet dynamics calculations studying state-to-state observables (e.g., Refs.10-13) typically transform between reactant and product coordinate systems while propagating using reactantproduct decoupling schemes. 48,49 n the present approach, initial wave packets representing the vibrational state of the activated complex are obtained by diagonalization of the thermal flux operator.These initial wave packets are then transformed into reactant and product coordinate systems.This approach has two advantages: the coordinate transformations have to be done only at one time and the wave packets to be transformed are compact and well localized.Since the subsequent separate propagations in reactant and product coordinates are connected via the thermal flux eigenstates providing the initial wave packets, all state-to-state scattering information can be obtained by these separate propagations and no further coordinate transformation is required.This feature of the new scheme follows straightforwardly from the equations derived in Sec.II and the technical details of the coordinate transformation in the context of MCTDH wave functions are discussed in Sec.III.
Following the description of the general theory and the MCTDH-related aspects of the present approach in Secs.II and III, respectively, the D + H 2 → HD + H reaction (for J = 0) will serve as an example and a very first application.Numerical details of the specific calculations are given in Sec.IV and the results are presented and discussed in Sec.V. Concluding remarks (Sec.VI) and an appendix discussing selected specific aspects will complete the article.

A. Quantum transition state concept
,40 Starting from flux correlation functions [45][46][47] the CRP can be written as (Note that atomic units and ¯= 1 are used throughout the present work.)Here, Ĥ is the Hamiltonian and F = −i[ Ĥ , h] measures the flux through an arbitrary dividing surface defined by the cutoff function h which vanishes at the reactant side of the dividing surface and equals unity on its product side.Locating the dividing surface close to the transition state geometry, an efficient description within the quantum transition state concept can be utilized to compute this quantity.To this end, one can employ the thermal flux operator 50 where T is a reference temperature and k denotes the Boltzmann constant.Hence, the CRP is evaluated as To calculate the CRP, one first has to find the eigenvalues of the thermal flux operator.This can be implemented using imaginary time propagation.As the thermal flux operator is a purely imaginary and hermitian operator, the eigenvalues f T come in pairs with opposite sign.2][53] For vanishing total angular momentum the number of non-zero eigenvalues of the thermal flux operator is small. 34To finally obtain the CRP, these eigenstates are then propagated in real time.Important advantages of this method are that, first, only few wave packets have to be propagated and, second, only short time propagation is required since the cumulative reaction probability only depends on the dynamics near the barrier.
For more details about the quantum transition state concept and its application to the calculation of cumulative reaction probabilities and initial state-selected reaction probabilities see, e.g., Refs.51-53.

B. Flux correlation functions and state-to-state observables
In the present work, an approach for the rigorous calculation of state-to-state reaction probabilities based on the quantum transition state concept will be developed.While working equations for the calculation of cumulative reaction probabilities and thermal rate constants 27,28,[30][31][32][33][34][35][36][37][38][39]41 or even initial state-selected reaction probabilities 29,32,40 can directly be derived from flux correlation functions, [45][46][47] the calculation of state-to-state reaction probabilities and cross sections within a transition state based approach is not as straightforward.
This can immediately be realized by counting the number of different dividing surfaces required in either calculation.To calculate the thermal rate constant or the cumulative reaction probability, one only has to be able to distinguish between reactants and products.Thus, only a single dividing surface is required in these calculations and the CRP can be written as in Eq. (1).
The calculation of initial state-selected reaction probabilities or cross sections requires the resolution of the specific quantum states of the reactants and therefore one dividing surface must be located in the asymptotic reactant region in these calculations.Thus, the calculation of initial state-selected observables within a transition state based approach requires two different dividing surfaces: one located close to the transition state and another one located in the reactant asymptotic area.Working equations can be derived straightforwardly from flux correlation functions 29,32,40 since two different flux operators, which can be related to different dividing surfaces, appear in flux-flux correlation functions.The initial state-selected reaction probabilities p n (E) then are calculated via where P n projects onto a specific initial quantum state, and F 1 and F 2 measure the flux through the surface located in the barrier region and the asymptotic surface, respectively.Proceeding towards the calculation of state-to-state reaction probabilities and cross sections, one finds now that two different dividing surfaces are needed to resolve the specific quantum states of the reactants and products: one located in the reactant asymptotic area and another one located in the product asymptotic area.Together with the dividing surface located close to the transition state, this makes a total of three different dividing surfaces required in the transition state based calculation of state-to-state observables.Thus, the number of required dividing surfaces exceeds the number of flux operators or dividing surfaces present in the existing flux correlation functions.Consequently, new functions involving flux operators are required to facilitate the calculation of state-to-state observables within the quantum transition state concept.
The identity provides a starting point for the present derivation.
Here, the notation [ Â(t)] Ĥ −e)t , Eq. ( 5) can be derived via where the substitution of the integration variables t = t 1 , t = t 2 + t 1 in the two-dimensional integration is employed.
In the next step, arbitrary incoming or outgoing wave packets initially located in the reactant or the product channel are considered.The incoming reactant wave packets, outgoing reactant wave packets, incoming product wave packets, and outgoing product wave packets are denoted + r , − r , + p , and − p , respectively, and satisfy The matrix elements of the propagator between these wave functions, , can be rewritten employing Eqs. ( 7)- (10) and the commutator Transforming the propagator e −i Ĥ t into the δ( Ĥ − E) distribution via Fourier transformation and using Eq. ( 5), the above matrix elements can be transformed from the time-dependent representation into an energy dependent one and read Since all above matrix elements either vanish or include a flux operator, they can be calculated by propagating eigenstates of the (thermal) flux operator therefore naturally fitting in a transition state based approach.Employing the thermal flux operator in its eigenstate representation, the right hand side of J. Chem.Phys.136, 064117 (2012) Eq. ( 16) can be cast into the form The matrix elements present in the equation can be calculated by propagating the thermal flux eigenstates into the reactant and product asymptotic region.Since − p |e −i Ĥ t |f T and f T |e i Ĥ t | + r vanish at negative times t (see Eqs. ( 7)-( 10)), only forward propagation of the thermal flux eigenstates is required to obtain the matrix elements of Eq. ( 16) and working equations are obtained as

C. Scattering matrix and state-to-state reaction probabilities
Scattering matrix elements can be calculated from the above matrix elements using the wave packet correlation function of Tannor and Weeks. 54Here, one considers wave packets j a = φ j a χ n a , a = r or p, (22)   which are products of a one-dimensional incoming (j = +) or outgoing (j = −) wave functions φ j a located in the asymptotic region depending only on the respective scattering coordinate and a wave function χ n a depending on all other coordinates.Note that j a and derived symbols depend on the internal state given by the quantum number n a , but for the sake of simplicity this dependence is suppressed throughout this work.Asymptotically, the Hamiltonian Ĥ can be separated in a one-dimensional Hamiltonian Ĥ 1D a describing motion in the respective scattering coordinate and a Hamiltonian Ĥ ⊥ a describing the remaining vibrational and rotational motion of reactants or products.χ n a is taken as a normalized eigenfunction of Ĥ ⊥ a with eigenenergy E ⊥ n a and n a are the corresponding quantum numbers.
The eigenstates of the complete Hamiltonian can then be constructed by a linear combination of these wave packets as 54 An energy normalized representation of j a,E can then be derived and reads 54 The η factors provide an appropriate energy normalization of the incoming and outgoing wave packets, Their phases are defined via η j a (E) = k|φ j a , where |k is an energy normalized plane wave state with the appropriate momentum k corresponding to the total energy E.
The correlation function expression that allows one to calculate the scattering matrix elements S n b n a (E) follows as 54 The scattering matrix elements depending on the total energy E are then given by Inserting Eq. ( 21) into Eq.( 27), one obtains an equation for the calculation of S-matrix elements which includes a flux operator measuring the flux through an arbitrarily positioned dividing surface.State-to-state reaction probabilities p n r →n p (E) are given as the squared modulus of the S-matrix elements Introducing the operators and employing Eqs. ( 27) and ( 25), the following equation for the p n r →n p (E) results: (For a more detailed discussion on the above formula and further formal developments see the Appendix.)

D. Asymptotic state analysis
The above equation combined with Eq. ( 21) could straightforwardly be applied to compute p n r →n p (E)'s within the quantum transition state concept.It employs arbitrary incoming and outgoing wave packets located in the asymptotic region to facilitate the initial and final state analysis.However, choosing specific analysis wave functions and revising the operators Pj r (E) and Pj p (E) employed in Eq. ( 30), one can derive particularly convenient equations as will be shown in the following.
The operator Pj a (E) consists of two different parts: P ⊥ n a is the projection operator which projects onto the specific internal and rotational quantum states considered in p n r →n p (E).This operator depends only on the asymptotic vibrational and rotational states under investigation and does not depend on the total or collisional energy.
The second operator present, Pj;1D a (E), acts only on the scattering coordinate R a describing the relative translational motion of the reactants or products.Due to the presence of the energy-normalizing factors η j a (E), it explicitly depends on the collisional energy This explicit energy dependence, which in particular complicates the subsequent computation of state-averaged quantities, can be avoided if an analysis wave function with an energy-independent normalization factor |η , where μ a is the reduced mass of the two colliding reactants or products and R a is their distance, appropriate analysis wave functions centered around R a = R 0 read where the function f(k) is defined as Here, f (k) can be any function with f (p max ) = 1 which decays sufficiently fast for k → ∞.These wave functions cover the complete momentum spectrum up to a maximum total momentum p max .For collisional energies E coll < (p 2 max /2μ a ) the absolute values of their energy-normalization factor equals unity: where

III. QUANTUM DYNAMICS
Combining Eqs. ( 21), (30), and (36), the working equation is obtained.The initial wave functions are provided by the thermal flux eigenstates localized in the transition state region.These are obtained using imaginary time propagation.Then the thermal flux eigenstates have to be propagated into the product and reactant asymptotic regions.The MCTDH approach 42,43 provides an efficient scheme for the imaginary and real time propagations.Since no single coordinate system allows one to efficiently simulate the dynamics in the reactant and product regions simultaneously, different coordinate systems are required.Consequently, one must transform the wave functions between two different coordinate systems.This can be achieved by initially transforming the thermal flux eigenstates to reactant and product type coordinate systems.This approach is particularly efficient since the thermal flux eigenstates are quite compact and the transformation has to be done only once at the initial time.
In this section the MCTDH approach will be introduced and the problem of coordinate transformation within this approach will be discussed.

A. MCTDH approach
The ansatz for a set of MCTDH wave function reads 42,43 where 1;k j k (x k , t) are time-dependent basis functions, called single-particle functions (SPFs).The A 1 j 1 ...j m ,w (t) are the expansion coefficients for the wave function in this basis.In this work the state-averaged MCTDH approach, 41 which simultaneously describes several MCTDH wave packets by using a common SPF basis, is employed.The different wave packets are denoted by the index w in Eq. ( 38).The singleparticle functions are represented by a time-independent basis J. Chem.Phys.136, 064117 (2012) One can use FFT and discrete variable representation (DVR) schemes here.The equations of motion are derived by employing the Dirac-Frenkel variational principle. 55Matrix elements of a general potential energy surface are obtained using the correlation DVR (Ref.56) scheme throughout this work.An efficient propagation scheme for the MCTDH approach, the constant mean-field (CMF), has been developed. 57In this work a revised version, the CMF2 scheme of Ref. 58, is used.

B. Coordinate transformation
Since the MCTDH approach employs a layered structure for the wave function representation, performing a coordinate transformation of the wave packets is not straightforward.The basic idea will be illustrated using the D + H 2 → DH + H reaction (more details about this example can be found in Sec.IV) and the transformation from hyper-spherical coordinates to Jacobi coordinates as an example.
The MCTDH wave function in hyper-spherical coordinates (ρ, α, ϑ) reads The SPFs 1;κ j κ (x κ ) can be represented in a basis of timeindependent functions as in Eq. (39).The basis functions χ κ i κ (x κ ) are known analytically and are related to the DVR or FFT grid representations by a unitary transformation (for details of the use of corresponding grid and basis representations see, e.g., Ref. 59).
The equivalent wave function in the Jacobi coordinate system is denoted by ¯ w (R , r , ϑ ) and reads Thus, a direct product grid of size N R × N r × N ϑ in the new coordinate system can be generated and the value of the wave function at each grid point (R k , r l , ϑ m ), k = 1, . . ., N R , l = 1, . . ., N r , m = 1, . . ., N ϑ , is given by Introducing DVR or FFT basis functions χ R k (R ), χ r l (r ), χ ϑ m (ϑ ) in the new coordinate system which are located at the corresponding grid points (R k , r l , ϑ m ), the wave function can be rewritten in the new basis as Note that for most DVRs weight functions have to be considered, which have been suppressed in the above equation for reasons of simplicity.Then this wave function is converted to a MCTDH representation.As a first step one adds as many SPFs as there are time-independent grid points in the respective coordinate.The wave function is thus written as with 1;3 Finally, the number of SPFs is reduced by analyzing the population of the natural orbitals.The natural orbitals are obtained by diagonalizing the single-particle density matrices 43 ρ(1;κ) .. ,j κ−1 ,m,j κ+1 ... j d .(48) Their eigenvalues give the populations of the uniquely defined natural orbitals.All the natural orbitals which are not significantly populated can be removed and an efficient MCTDH representation is obtained.
Before closing this section some technical details have to been mentioned which previously have been disregarded for the sake of clarity.The MCTDH implementation used in this work incorporates the volume element in the radial coordinates, v = rR or v = r R , into the wave function, i.e., ˜ = √ v .Hence, the change of the volume element has to be taken into account.It should also be noted that the above procedure results in a wave function with undefined norm.However, this does not constitute a problem since the norm of the initial wave functions can be stored and recovered after the transformation.
For high-dimensional systems no direct product grid in all coordinates can be stored.However, in a polyatomic reaction typically only a subset of the coordinates need to be transformed.Considering, for instance, the H + CH 4 → H 2 + CH 3 reaction with its 12 internal degrees of freedom, for example, the six coordinates describing the methyl group do not change and only the remaining six coordinates have to be transformed.Here one only would need to construct a direct product grid in these six coordinates.

IV. SYSTEM AND NUMERICAL DETAILS
For a first numerical test of the described method, the state-to-state reaction probabilities of the D + H 2 → DH + H(J = 0) reaction are calculated.The present calculations employ the BKMP2 (Ref.60) potential energy surface.

A. Coordinate systems
Three different coordinate systems are employed in the calculations.For the real time propagation reactant and product Jacobi coordinates are used.Here, r and r denote the bond distance of the diatom (H 2 in case of reactant Jacobi coordinates and DH in case of product Jacobi coordinates), R and R denote the distance of the third atom to the center of mass of the diatom, and ϑ and ϑ are the angles between the two distances, respectively.The calculation of the flux eigenstates employs hyper-spherical coordinates based on the reactant Jacobi construction: where μ r and μ R are the corresponding reduced masses and the unchanged angle ϑ as third coordinate.
The kinetic energy operators employed read and

B. Thermal flux eigenstates
The flux operator is employed.denotes the Heaviside function.Here, a Heaviside function is employed to eliminate flux at ϑ angles larger than π/2.Thus, only the reaction of one of the two equivalent hydrogen atoms is explicitly considered.The dividing surface is located at α 0 = 0.45 and up to 20 flux eigenstates are considered in the state-averaged MCTDH calculations.A reference temperature for the thermal flux operator of 4000 K is chosen.All other details are summarized in Table I.The thermal flux eigenvalues obtained in the largest calculation are listed in Table II.

C. Real time propagation
The thermal flux eigenstates are transformed to the reactant and product coordinate systems and propagated 250 fs into the reactant and product asymptotic regions, respectively.For details of the wave function representation see Table III.The number of grid points and SPFs guarantees converged results.
Quartic absorbing potentials of the form are used.x a denotes the start and x e denotes the end of the absorber.The absorbing potentials act in r and R in the reactant coordinate system and in r and R in the product one.In all calculations presented below, the end of the absorbers is placed at the end of the grid.The strength parameter c is chosen to be 5 eV throughout this work.x a = 3.1 bohr is used as a starting location for the absorbers in r and r in the reactant and product coordinate system, respectively.For the absorbers in R and R , a starting location of x a = 30 bohr in both coordinate systems is used.For higher energies also the insertion channel begins to open and small parts of the wave packet start to move into the region of small R and R .Since the kinetic energy operators in reactant and product coordinates show singularities at R = 0 and R = 0 , respectively, numerical difficulties could result.
To avoid these problems a barrier of 10 eV is introduced at  Legendre-DVR 20 56 0.0 -π R = 0 and R = 0. Consequently, the insertion channel may not be properly described in the present calculations.
The dividing surfaces for the operators P+ r and P− p are located at R 0 = 7.5 bohr and R 0 = 7.5 bohr, respectively, and 60 ro-vibrational states for reactants and products are considered.Convergence tests have been performed to assure that the analysis wave functions are located sufficiently far in the asymptotic regions.The function f (k) occurring in Eq. ( 34) is chosen as with p values of 2.47 a.u. and 4.95 a.u.for the reactant and product channels and employing a p max value of 12.2 a.u..

D. ABC calculations
The reference data is obtained performing timeindependent scattering calculations with the program ABC. 61 The parameters used for the J = 0 calculations can be found in Table IV.For more details of these parameters see Ref. 61.

A. Cumulative reaction probability
As a starting point, the cumulative reaction probability N(E) is calculated by summing initial state-selected reaction probabilities calculated using a scheme based on Eq. ( 4) and detailed in Ref. 25.After transforming the thermal flux eigenstates computed in hyper-spherical coordinates to a reactant or product Jacobi coordinate system, respectively, the CRP is calculated in the respective Jacobi coordinates and the second dividing surface is located at R d = 7.5 bohr or R d = 7.5 bohr.The results are shown in Fig. 1. Results obtained with the ABC program are given for comparison.Up to an energy of 1 eV, the results are in very good agreement with the timeindependent scattering calculations indicating the reliability of the coordinate transformation.The agreement also proves that a sufficient number of ro-vibrational states in the reactant and product channels are included in the present calculations.Furthermore, the cumulative reaction probability is calculated by summing all state-to-state reaction probabilities computed according to Eq. (37).Here the results are in excellent agreement for energies up to 0.8 eV.For higher energies good agreement is found, however, tiny deviations are visible at about 0.85 eV.

B. Convergence with respect to flux eigenstates
The convergence with respect to the number of eigenstates of the thermal flux operator used throughout the calculation is investigated in Fig. 2, displaying the D + H 2 (ν = 0, j = 0) → H + DH(ν = 0, j = 0) reaction probability as a function of the total energy.Considering the whole range of energies up to 1 eV, 6 pairs of thermal flux eigenstates are required to obtain convergence.Including only three pairs of eigenstates convergence up to an energy of about 0.85 eV is obtained.This result can be rationalized using a harmonic model of the transition state.The harmonic frequencies at the transition state are ω b = 868.7 cm −1 (bend) and ω s = 1768.6cm −1 (symmetric stretch).Since only even quantum numbers in the bending mode are allowed for J = 0, three pairs of thermal flux eigenstates include the vibrational states (0, 0), (2, 0), (0, 1), where (n b , n s ) denote the quantum numbers in the bending and symmetric stretching mode, respectively.Thus, in harmonic approximation vibrational states up to an excitation energy of 1ω s = 1768.6cm −1 = 0.22 eV are included.The state (4, 0) with an excitation energy of 3474.8 cm −1 = 0.43 eV is the lowest one neglected.Thus, the results obtained with three pairs of thermal flux eigenstates can be expected to yield converged results for energies up to about 0.4 eV above the threshold for reaction.This is consistent with the findings of Fig. 2. Using the same arguments, one can expect that six pairs of thermal flux eigenstates are sufficient to obtain converged results for energies up to about 0.6 eV above the threshold.Energy / eV ν=0, j=0, ν'=0, j'=0 ν=0, j=0, ν'=0, j'=1 ν=0, j=0, ν'=0, j'=2 ν=0, j=0, ν'=0, j'=5 FIG. 5. State-to-state reaction probabilities for D + H 2 (ν = 0, j = 0) → H + DH (ν = 0, j ).The ABC results are given as symbols.

C. State-to-state reaction probabilities
Figure 2 also shows that the D + H 2 (ν = 0, j = 0) → H + DH(ν = 0, j = 0) reaction probabilities computed with the present approach (Eq.( 37)) agree well with the reference results obtained with the ABC scattering for energies up to about 0.8 eV.For energies above 0.8 eV, the present approach starts to yield reaction probabilities somewhat smaller than the ABC reference results.These differences increase up to about 15% at the upper end of the energy scale, 1.0 eV.However, one should note that the magnitude of this probability difference, 0.004, is small on an absolute scale and thus consistent with the tiny differences seen Fig. 1.We assume that this difference is mainly caused by the neglect of the insertion channel which starts to open in this energy range.Furthermore, it should be noted that the accurate calculation of very small probabilities is always a difficult task for a mean-field based approach like the MCTDH scheme which uses non-linear equations of motion.Inspite of careful convergence testing, as done for the present calculations, tiny remaining errors in the propagated wave function might remain undetected.In the present calculations, numerical inaccuracies generally tend to become more relevant for higher energies due to the thermal weighting used in the present state-averaged MCTDH propagation of the thermal flux eigenstates.
It is thus interesting to study how different state-to-state reaction probabilities are affected by numerical inaccuracies.To this end, different partial summed and state-to-state reaction probabilities are shown in Figs.3-7.In all figures, the results are in good agreement with the benchmark results up to an energy of about 0.8 eV.For higher energies the results tend to differ more strongly from the ABC results.These findings are very similar to the ones already seen in Fig. 2 for the D + H 2 (ν = 0, j = 0) → H + DH(ν = 0, j = 0) reaction probability.It is noteworthy that results summing over final vibrational (Fig. 3) or rotational channels (Fig. 4) show similar accuracy as the corresponding completely state resolved probabilities of Figs. 5 and 6, respectively.Moreover, reaction probabilities starting and/or ending in vibrational excited states (Figs. 4 and 6) are obtained with about the same accuracy.This also holds true for fully state-resolved reaction probabilities starting in a rotational excited state and ending in different rotational excited states (Fig. 7).Thus, for a given total energy all state-to-state reaction probabilities are obtained with comparable accuracy.

VI. CONCLUDING REMARKS
A new method for the calculation of state-to-state reaction probabilities is proposed in this work.Analytical formulas required for the calculation of S-matrix elements and state-to-state reaction probabilities within the quantum transition state concept are derived.The numerical approach proposed employs the eigenstates of the thermal flux operator, which are located near the transition state.These flux eigenstates are then transformed to reactant and product coordinate systems and subsequently propagated into the asymptotic regions, where they are analyzed with respect to the asymptotic states.An analysis scheme specifically adapted to this approach has been derived.State-to-state reaction probabilities are finally calculated by combining the results of both asymptotic analysis.The new method was tested for the D + H 2 → DH + H(J = 0) reaction and validated by a comparison to exact results calculated with the ABC program.
While existing wave packet dynamics approaches for state-to-state reactive scattering employing initial wave packets corresponding to specific quantum states of the reactants and localized in the reactant asymptotic region, the new scheme presented here employs initial wave packets located in the transition state region.This has distinct advantages: First, one obtains the full S-matrix from propagating only a small set of wave packets defined by the thermal flux operators.Employing the state-averaged MCTDH approach, only two calculations are required.Second, the coordinate transformation into reactant or product coordinates has to be done only once.Third, one can decouple different reaction channels at the transition state and use coordinate systems adapted only to reaction via a specific transition state geometry.
However, also an important restriction resulting from this strategy must be mentioned.Since one has to account for all vibrational channels through the transition state which are open at the total energy considered, the new scheme proposed is practically limited to the calculation of reaction probabilities at total energies which exceed the threshold energy of reaction only moderately.At higher energies, the number of accessible channels through the transition state becomes prohibitively large for polyatomic reaction.Additionally, due to the thermal weighting employed in the definition of the initial wave packets via the thermal flux operator, the numerical accuracy required in the MCTDH calculations also tends to become a critical factor at high energies.
The new scheme can be applied to direct polyatomic reactions proceeding via a potential energy barrier and facilitates further developments of the quantum transition state concept.The previous work 24,25 demonstrated that the quantum transition state concept permits the calculation of initial state-selected (J = 0) reaction probabilities for the H + CH 4 → H 2 + CH 3 reaction in full dimensionality.Using the scheme presented, these calculations could now be extended to also compute state-to-state reaction probabilities for this system.
The present work focused on the calculation of reaction probabilities for J = 0.The approach could straightforwardly be used to compute S-matrix elements required for the calculation of differential cross sections.However, proceeding to the calculation of integral and differential cross sections, which are of particular interest in comparison with experiment, S-matrix elements for all values of the total rotational quantum number J are required.Work along these lines is in progress.

TABLE I .
Number of grid points (N) and SPFs (n) used in the flux eigenstate calculation.

TABLE II .
Eigenvalues of the thermal flux operator.

TABLE III .
Number of grid points (N) and SPFs (n) used in the propagation to the reactant and product asymptotic region.

TABLE IV .
Parameters used in the ABC calculation.FIG. 1. Cumulative reaction probability of D + H 2 → DH + H(J = 0) calculated in the two different Jacobi coordinate systems with the second dividing surface located at R d = 7.5 bohr or R d = 7.5 bohr and by summing all state-to-state reaction probabilities.The ABC results are given as symbols.