Reaction cross sections and thermal rate constant for Cl + CH3Br ClCH3 + Br from J-dependent quantum scattering calculations

Employing dimensionality-reduced time-independent quantum scattering theory and summation over all possible total angular momentum states, initial-state selected reaction cross sections for the exothermic gas-phase bimolecular nucleophilic substitution (SN2) reaction Cl + CH3Br ClCH3 + Br have been calculated. The carbon–halogen bonds and the rotation of the methyl halides are taken into account. In agreement with previous calculations for J = 0, initial rotational motion of CH3Br decreases the reaction probability and consequently the cross sections. The experimentally obtained thermal rate constant for 300 K is reproduced within the experimental error. For lower temperatures, it is calculated to be below the experimental values but shows the same strong increase for T 0.


I. Introduction
Nucleophilic bimolecular substitution (S N 2) reactions [1][2][3][4] which are characterized by the simultaneous making and breaking of two single bonds are not only of high practical interest but also can contribute to elucidate the very fundamentals of chemical reactions. In particular, this holds for gas-phase reactions which are not complicated by the surrounding solvent molecules. Extensive material, both experimental and theoretical, is available on the simple halogen exchange S N 2 reactions (with halogen atoms X and Y as nucleophiles and leaving groups).  As is well-known, [1][2][3][4]6,9,14 these reactions are governed by two relatively deep wells (ca. 4000 cm À1 ) before and after the central reaction barrier in the respective potential energy profiles, originating from the strong electrostatic attraction between ion and molecular dipole, yielding the two complexes X À Á Á ÁCH 3 Y and XCH 3 Á Á ÁY À . It should be noted that also hydrogen-bonded complexes, e.g. the Cl À Á Á ÁHCH 2 Br prereaction complex and the ClCH 2 HÁ Á ÁBr À post-reaction complex play an important role. 65,66 While identity S N 2 reactions (X = Y) are characterized by low reaction rates due to sharp barriers in the reaction profiles, exothermic halogen exchange reactions usually exhibit a barrier located below the reactant energy. From a naive point of view, ion-dipole capture theory should suffice to describe such a reaction, but even in statistical theories and classical mechanical treatments it is shown that a central barrier lower than the energy of the reactants influences the reaction due to large collision angular momenta and the corresponding large rotational energy (centrifugal) barriers which are added to the central barrier height. 67 Comparisons between experiment, chemical dynamics, and ion-molecule capture rate constants are made for OH À + CH 3 I in ref. 68 and for OH À (H 2 O) + CH 3 I in ref. 69. Good agreement is found between experiment and chemical dynamics, and for OH À + CH 3 I it is quantitative. The ionmolecule capture rate constants are too large. For these reactions the central barrier is about 20 kcal mol À1 lower in energy than that for the reactants.
Quantum-mechanically the potential wells reflect incoming wave packets, thus setting high demands on the theoretical description. The wells support considerably long-lived metastable states bound by intermolecular ion-dipole forces which in the reaction probability vs. energy profiles show up as extremely narrow lines and strongly influence not only differential and total reaction cross sections but also the overall rate constant. These long-lived states can be of dynamical nature (energy stored in higher-frequency modes: Feshbach resonances) or, in exothermic systems, e.g. Cl À + CH 3 Br, additionally can be quasi-bound by the potential well (shape resonances). A quantum-mechanical approach is indispensable to calculate state-to-state reaction probabilities which then can be summed up to differential or integral cross sections and finally turned into rate constants.
In this paper, we present a quantum study of the exothermic reaction Cl À + CH 3 Br (a-channel) -ClCH 3 + Br À (b-channel). ( Experimentally, Viggiano et al. 12 found that in the Cl À + CH 3 Br reaction a higher CH 3 Br temperature does not increase the thermal rate constant, and Graul and Bowers 11,13 reported that the products are vibrationally hot. Angel and Ervin 21 found that the S N 2 mechanism is ineffective compared with predictions from phase-space theory and ion-dipole capture theory. However, experiments on the detailed, state-selected dynamics of this reaction are scarce. The thermal rate constant down to 23 K was measured by LeGarrec et al. 43 who also report a calculational result based on a dimensionality-reduced 3D quantum model, employing an approximate Hamiltonian for the umbrella bending mode of the methyl group and the rotating bond approximation (RBA), yielding a good, yet not perfect agreement with the experimental data. Reaction cross sections in a 2D model were calculated and analyzed by Schmatz and Clary. 46 Even in these dimensionality-reduced models, the low intermolecular frequencies of the complexes, the heavy mass of the bromine nucleus and the far-reaching ion-dipole potentials render the Cl À + CH 3 Br S N 2 system computationally very demanding because both memory and CPU time requirements set limits. The small rotational constants strongly increase the densities of both reactant and product states. Thus, one has to resort to dimensionality-reduced models 42,48 where only those degrees of freedom are included which are absolutely indispensable for a qualitatively correct description of the reaction coordinate. In a second step, the influence of additional degrees of freedom can be investigated in detail. This analytical approach is only viable by theory whereas experiments deliver a more or less complete information, often in the form of highly averaged quantities like thermal rate constants. It has to be pointed out that the molecular beam ion-imaging experiments of the Wester research group have provided quite detailed non-thermal dynamics for S N 2 reactions. [22][23][24][25][26][27][28][29] While in our previous 4D calculations carried out under the constraint of C 3v symmetry, 55,57,[59][60][61] i.e. with inclusion of the umbrella bending and C-H stretching vibrations, a lot of information on mode-selectivity and the influence of the excitation of different internal degrees of freedom on the reactivity was obtained, the overall rate constant could not be reproduced exactly. The calculated value turned out to be too small 57 because the computed reaction probabilities were determined to be too low, and important resonances, in particular broad ones were missing. As the employed model was exact within the dimensionality-reduced framework and all computations were well converged, the reason for the discrepancy can only be found just in this reduced dimensionality. Thus, in ref. 62 rotation of both reactant and product methyl halide molecules was included in state-to-state reaction probability calculations to properly take into account translation-to-rotation energy transfer. While the rotation around the axes perpendicular to the axes of molecular symmetry was included, the two internal totallysymmetric vibrational modes (C-H stretch and umbrella bend) were not longer taken into account. Particular numerical challenges in this type of calculation are (i) the number of energetically accessible states which is enhanced by one order of magnitude to 10 4 -10 5 , (ii) the necessity to propagate the wavefunctions out to much larger distances where the interaction between the reactant and product channels is zero, and finally (iii) the very fine energy resolution imposed by the narrow resonances.
In contrast to the reaction probability P(E), the cross section s(E) is a quantity that is relatively easily accessible by appropriate experiments. State-selected reaction cross sections provide valuable information on both the influence of initial vibrational and/or rotational excitation. In the present work, we report on total angular-momentum dependent quantum scattering calculations for initial state-selected total cross sections s(E) and the thermal rate constant k(T) for reaction (2). The cross sections are obtained as weighted sums over J-dependent reaction probabilities P J (E) for given energy E and all possible total angular momentum quantum numbers J.
The paper is organized as follows: In Section II, we briefly report on the theoretical background, in particular the employed Hamiltonian, the reactive scattering formalism and the reduceddimensionality theory to recover the full-dimensional thermal rate constant while in Section III numerical details of our computations are given. Section IV presents the results and their discussion and Section V contains our conclusions. Throughout, energies are quoted in wavenumber units.

II. Theory
We use only a single Jacobian coordinate system {R, r, g} connected to the product channel with the angle g describing the rotation of the pseudo-diatom CH 3 Cl or, in other words, the position of Br À with respect to this moiety. With the reduced masses m 1 = m Br m CH 3 Cl /(m Br + m CH 3 Cl ) and m 2 = m Cl m CH 3 /m CH 3 Cl , the full four-dimensional Hamiltonian for this pseudo-triatomic system in the space-fixed frame 70,71 can be simplified by transforming the wavefunction C into c/(Rr), yieldinĝ In eqn (3), the squared orbital angular momentum operator l 2 (as a function of its two defining angles in the space-fixed coordinate system, y and j) is given bŷ while a corresponding equation for the angular momentum ĵ of the pseudo-diatom holds with g taking the role from y in eqn (4). Employing total angular momentum conservation, Ĵ = l + ĵ, and changing to the body-fixed frame with reference axis z 0 (along the Jacobi vector R, thus l z = 0) we use with raising and lowering operators Ĵ AE and ĵ AE . Arranging terms and integrating over spherical harmonic functions the Hamiltonian reads: 70,71 where 72 Here, | jki denotes a spherical harmonic function and hÁ Á Ái g indicates integration over g, while J, j = |k|, |k| + 1,. . . and k = ÀJ,. . ., J are the quantum numbers of the conserved total angular momentum, of the rotating diatom and of the projection of the total angular momentum Ĵ onto the body-fixed axis defined by the third atom and the center of mass of the diatom.
In this work, we treat the rotation of the reactive system for J Z 0. The infinite rotational energy of the third atom with respect to the center of mass of the diatom is responsible for a singularity. Only in the J 4 0 setting, this is a coordinate singularity due to the non-Cartesian coordinate system (Eckart singularity).
In Jacobian coordinates each of two variables R and r can become infinite in different asymptotic configurations. Thus, a transformation to a single unbound reaction coordinate only, the coordinate-system independent hyperradius r (ref. [73][74][75][76][77][78] is applied, where d = arctan(r/R) is the corresponding hyperangle and R and r are scaled according toR ¼ ffiffiffiffiffiffiffiffiffiffi The Hamiltonian operator then becomeŝ with the parametrically r-dependent surface Hamiltonian Including nonvanishing total angular momenta, we have now and We compute the wavefunctions based on only one of the (biased) hyperspherical coordinate systems in the interaction region, a method already employed directly in Jacobi coordinates for the standard reaction D + H 2 -DH + H. 79 All arrangement configurations still can be described accurately by a proper basis set adaptation. For an exact Hamiltonian, at fixed hyperradius r, all arrangement channel coordinates yield the same eigenstates (with different parametrization), such that matching onto the asymptotic coordinate systems is feasible once the interaction region is left. An additional intermediate region is used to smoothen the transition between these two regimes.
To describe the eigenfunctions of Ĥ surf , discrete variable representations (DVRs) [80][81][82] are employed for d and g; while for the hyperangle d, the ''sinc-DVR'' 83,84 is chosen, for the physical angle g the Legendre polynomial basis {P j (g)} is employed. As within each sector, the Ĥ surf matrices of size up to 30 000 render a direct diagonalization very difficult, we employ an iterative method by taking advantage of the relative sparsity of the matrices in matrix-vector multiplications and calculate the eigenvalues and eigenvectors of the Ĥ surf matrices using Lanczos iteration. 85,86 The basis for the matrix representations of Ĥ surf was constructed making use of potential-optimized DVRs (PODVRs) 87 in both d and g.
In For given J and total energy E, the time-independent Schrödinger equation Ĥ J C J = EC J , is solved in two steps to obtain the partial waves C J . After expanding the partial waves C J n 0 for initial vibrational state n 0 in close-coupling form with N channels, 55,59 the eigenvalue equation for Ĥ surf is solved. The propagation of the wave functions in all three regions is carried out according to the R-matrix propagation algorithm. 88 The eigenfunctions in hyperspherical coordinates and the corresponding R-matrices are projected from the representation in hyperspherical coordinates onto that in Jacobian coordinates 89 and the propagation is continued in the specific Jacobian coordinates for the a and b configuration. In the asymptotic region, larger sector sizes can be chosen, and we used a step size algorithm similar to the description in ref. 88 to determine appropriate sector sizes. Finally, proper boundary conditions were applied to extract the complex S-matrix elements S J i, f (E), employing the single-sector projection method.
When only reaction probabilities for J = 0 are available, the centrifugal term Ĥ J = Ĥ + h 2 J( J + 1)/2mr 2 can be added to Ĥ J=0 in the propagation of the scattering wave function (rotating line approximation (RLA) 90 ). To be more accurate, here we make use of probabilities P J i, f (E) = |S J i, f (E)| 2 calculated employing Ĥ J from eqn (10) for all possible values of J. Reactive scattering cross sections are then computed according to with k i given by where e i is the asymptotic energy of channel i, and E À e i = E trans . While from the full-dimensional cumulative reaction probability the exact thermal rate constant could be obtained, the dimensionality-reduced theory based calculation of k(T) requires energy shifting procedures. 89,[91][92][93][94][95] However, since the saddle point in the Cl À + CH 3 Br reaction is located below the reactant energy, usual transition state theory cannot be applied due to a formally negative barrier height. However, characteristic data of the transition state as harmonic vibrational frequencies and rotational constants can be used in the calculation of partition functions. The azimuthal rotation in the transition state is included in the sum over the total angular momentum J. In a modification of eqn (13) from ref. 61, we calculate k(T) according to Here, the quantum number K gives the body-fixed projection of the total angular momentum on the molecular axis of symmetry, A † and B † are the rotational constants of the prolate symmetric top saddle point geometry [ClÁ Á ÁCH 3 Á Á ÁBr] À , and Q y red 0 denotes the partition function of the two non-degenerate totally symmetric modes (C-H stretch and umbrella bend) and the three doubly degenerate internal modes (antisymmetric C-H stretch, antisymmetric C-H bend and the ''rocking motion'') in this transition state structure. Note that the degenerate bending mode is explicitly included in our treatment of the rotation. The difference of the zero-point energies between the saddle point structure and the reactants with respect to the two a 1 modes and the three e modes is omitted due to the absence of a positive net barrier height. Q trans denotes the three-dimensional translational partition function and Q int the partition function for all internal degrees of freedom of the reactants.

III. Numerical details
The potential energy surface of Hase et al. 34 was used and optimized with respect to the position of the CH 3 group. 62 A local symmetry approximation (LSA) was employed keeping the four nuclei of the methyl group always in the geometry of a right circular cone. 62 The final sector-dependent potential-optimized (PO) basis sets in d and y were constructed as described in ref. 62. The number of basis functions in both coordinates and the respective energies E PO for the interaction, intermediate and asymptotic regions are collected in Table 1 of ref. 62. The cutoff energy was set to V cut = 21 000 cm À1 . Depending on the region, we either employed the Hamiltonian (10) in hyperspherical coordinates with kinetic energy operators eqn (12) and (13) or Ĥ according to eqn (6) in scaled Jacobi coordinates (R,r) with kinetic energy operators according to eqn (7) and (8). All channels up to a total energy of E diag max = 2500 cm À1 have been computed and included in the R-matrix-propagation with a maximum number of N ch,max = 1466 channels. All channels below 4500 cm À1 were additionally included in each sector in computations with an energy below this limit. The calculations started at r = 6.5 a 0 in the interaction region and switched to the intermediate region at r = 13.028 a 0 (where the actual separation of the wave functions according to eqn (15) has already occurred at lower values of r and a minimum separation value of r = 13 a 0 was imposed). At r = 40.889 a 0 , the transition to the asymptotic region took place. The wave functions were propagated to R asym = 3316.575 a 0 . The final S-matrix was computed at energy intervals of DE = 5 cm À1 from E min = 300 cm À1 to E max = 2000 cm À1 . The convergence parameter for the localization of the wavefunctions, e sep , was set to 2 Â 10 À6 . In the computation of k(T), the summation over K runs from ÀN to +N. This has only a very small influence because the sum over K is fast converging (at least from J = 15) and the probabilities considerably decrease from J = 100. The summation over J is carried out until ð2J þ 1Þ P f P J i;f ðEÞ o 10 À6 . Unless stated otherwise, all energies are counted from the asymptotic potential energy value of CH 3 Br and are quoted in wavenumber units throughout. The convergence of the calculated cross sections was subject to careful examination, particularly with respect to the final value of the hyperradius and the number of partial waves considered in the expansion of the wavefunction. The geometrical parameters employed to calculate the symmetric top rotational constants as well as the vibrational frequencies used in the k(T) computation are given in Tables 1 and 2.

IV. Results and discussion
The integral cross sections for the reaction Cl À + CH 3 Br -ClCH 3 + Br À for total energies up to 1900 cm À1 , calculated within the dimensionality-reduced model, are graphically This journal is © the Owner Societies 2016 displayed in Fig. 1 and 2. Initial state-selected total cross sections, summed over all product states, are shown. Due to the summation over all partial waves that contribute to s(E) for a given total energy, a first automatic averaging occurs; an additional smoothing of the s(E) curves was carried out according to our earlier work on reaction probabilities P(E) for J = 0. 62 Even though many partial waves are summed up, a coarse-grained oscillating pattern is still seen in the s(E) curves. However, in some cases, the smoothing results in plateau-like structures which have to be considered as artifacts.
The low-energy cross sections contribute most to the thermal rate constant. Consequently, a high resolution with respect to energy is necessary in the calculations.
Most of the cross sections show a typical well-known behaviour: 49 An immediate rise when the channel opens and then a decay that behaves inversely proportional to k m 2 or E trans , respectively. The sharp maximum of s(E) is located close to that energy where the channel opens. The cross sections can thus very well be compared with respect to their magnitude, much better than the corresponding reaction probabilities. 62 It is important to note the logarithmic scale in Fig. 1 and 2 when comparing the amplitudes of the cross sections. As can be seen from Fig. 1, the integral reaction cross sections strongly decrease for rotational excitation in line with the observations made in ref. 62. It is remarkable that throughout the shapes of the different cross section curves in Fig. 1 are very similar, in line with the rotating line approximation that just shifts the energies by including a centrifugal term in the Hamiltonian which is proportional to the rotational energy of the whole system. This effect is clarified by considering the broad maximum in s (0,0) (E) at about 1400 cm À1 that gradually shifts to higher total energies up to 1750 cm À1 for j = 40. The shifting towards higher energy, correlated with J( J + 1), defines a trend which, however, is somewhat too oversimplified. E.g., the respective maximum for j = 4 is found below that for j = 5.
This indicates that the quantum character still prevails and that resonance effects between partial waves play a very important role, giving rise to deviations from the simple rotating line approximation picture.
Due to the data processing and averaging, the well-resolved resonance structure is lost. The plateau-like features seen in the plots originate from the data sampling. Interesting effects in the cross sections are the rectangular features in the cross sections for j = 0, j = 1 and j = 2. The plateau in the range around 1000 cm À1 does not show up in the comparable energy regimes for the other initial-state selected cross sections. The same holds for the broad peak located between 600 and 800 cm À1 . On the other hand, the broad feature at about 1700 cm À1 is also found in the cross section for j = 2 at somewhat higher energies. In the curve for j = 1, however, this pattern is not seen in between these two energy values, but instead shows up at an even higher total energy.  1 Integral reaction cross sections summed over all product states for the reaction Cl À + CH 3 Br (v,j) -ClCH 3 + Br À , for different initial rovibrational states (v = 0,j) with initial rotational excitation of CH 3 Br up to j = 40. Fig. 2 Integral reaction cross sections for the reaction Cl À + CH 3 Br (v,j) -ClCH 3 + Br À , for different initial rovibrational states (v,j) with initial rotational excitation of CH 3 Br up to j = 40, summed over all product states. Results are shown for v = 1 and 2 quanta in the C-Br stretching vibration and various additional initial rotational excitations of CH 3 Br.
A further interesting feature is that the overall shape of the reaction cross sections changes dramatically for very high angular momenta, e.g. for j = 40. Note that the cross sections for high rotational excitation do not show a pronounced rise directly at the opening of the channel. As the sharp rise is lost, the s(E) curve resembles much more typical cross section vs. energy curves in neutral-neutral reactions.
In Fig. 2, reaction cross sections for additional excitation of the C-Br stretching vibration v, localized in the coordinate of the bond being broken, is shown. For comparison, the cross section out of the rovibrational ground state (0,0) is also displayed. For excitation in the C-Br stretching vibration with one or two quanta, the cross sections rise very strongly as known from our 2D calculation presented in ref. 49. The effect of rotational motion in the reactant molecule is the same as discussed before, i.e. in general a decrease is observed for higher initial j values. It is worth noting that for higher initial j, the maxima of the cross sections become smaller, but as a function of energy no further decay of the functions is observed. There is a sharp peak when the channel opens followed by a moderate decay, but then the cross section curves remain on average almost constant -apart from coarse oscillations -for a relatively broad energy interval.
The thermal rate constant for the reaction Cl À + CH 3 Br -ClCH 3 + Br À as obtained from this work is graphically displayed in Fig. 3. The plot shows an atypical non-Arrhenius behavior with a negative temperature coefficient, which, however, strongly resembles that for ion-dipole reactions without a barrier. The older experimental thermal data, measured by Tanaka et al., 5 Olmstead and Brauman, 6 Caldwell et al., 7 Bohme and Raksit 8 as well as Viggiano et al. 12 are very well met by the results of the previous calculation.
For lower temperatures, however, quantitative agreement with the experimental data could not be reached. While experiment reveals a strong negative temperature dependence, a dramatic increase of the rate constant by more than two orders of magnitude going down from 300 K to 23 K, the calculations give only a one-order-of-magnitude increase. This is a similar behaviour as obtained by Clary in his RBA treatment 43 which, however, yielded results somewhat closer to experiment. While Clary made use of the same potential energy surface 34 as employed in our work, he scaled the geometries of the C-Br and C-Cl distances close to the saddle point structure.
The reasons for the discrepancies of both theoretical treatments are the relatively low quality of the potential energy surface which however is full-dimensional. Another potential surface, based on higher-level ab initio data may provide better results. As long as full-dimensional quantum scattering calculations are not possible, dimensionality-reduced settings have to be employed. These however rely on reduced-dimensional potential surfaces, optimized with respect to the remaining coordinates. The computation of optimized potential energies on the fly is still very costly.
A further point are the already mentioned approximations in theory. While Clary included the umbrella bend, we explicitly considered the degenerate transition-state bending vibration which has a much lower frequency leading to considerably more adiabatic channel curves in the scattering calculation. The respective other modes are taken into account by transition state theory, i.e. via the partition functions and zero-pointenergy differences, as are the remaining modes which are often considered to be spectators. These have much higher frequencies and are relevant for higher temperatures.
Clary pointed out that the quantum treatment of the umbrella bend is responsible for increase of k(T) for low temperatures because the adiabatic curves can have activation energies which may become relevant at higher T. In this work, we showed that inclusion of the degenerate transition state bending motion has the same effect. This means that the general behaviour of the k(T) curve is dominated by the two-dimensional scattering, including only C-Br and C-Cl coordinates. 49 The total cross sections reveal sharp peaks for low translational energies which are responsible for the behaviour of the rate constant that increases for lower temperatures. This can be explained by the scattering cross sections which in general become smaller for higher translational energies.

V. Conclusions
We have carried out time-independent quantum scattering calculations for the gas-phase reaction Cl À + CH 3 Br -ClCH 3 + Br À for all possible total angular momentum values up to energies of about 2000 cm À1 above the asymptotic potential energy in the reactant channel. The converged S-matrices were used to compute initial-state selected total reaction cross section and therefrom, by thermal averaging, also rate constants. The latter are the only quantities for which reliable experimental data are available. In the past, several degrees of freedom have been included in addition to the distances of the bonds being formed and broken (which are indispensable for a quantum-mechanical description of the reaction): the umbrella bending mode Fig. 3 Thermal rate constant for the reaction Cl À + CH 3 Br -ClCH 3 + Br À as calculated in this work (full line). Also shown are calculational results carried out by Clary 43 and several experimental data. [5][6][7][8]12,43 Note the logarithmic scale of the ordinate.
(within an approximate treatment), 42 the azimuthal rotation of the methyl halides, 45 the symmetric C-H stretching mode and the umbrella bending mode (exact within in the reduceddimensional treatment) 55 and in the present paper the doublydegenerate low-frequency bending mode. 62 The latter, which corresponds to the rotational motion in the reactant or product molecules, should be expected to be most sensitive in a rate constant calculation due to its small vibrational frequency in the transition state which leads to a considerable population of excited internal states. However, it can be seen that the thermal rate constant is mostly influenced by the simple 2D picture which already yields the characteristic shapes of the reaction cross section curves. 49 The additional degrees of freedom are responsible for a stepwise more and more accurate description of the rate constant. It is expected that the full cumulative reaction probability will give a rate constant that perfectly agrees with the experimental results. However, this cannot be the main target of the research within quantum scattering theory applied to molecular encounters. Much more valuable are the detailed features that originate from the individually investigated degrees of freedom. From a quantitative point of view, the most important issue will be the availability of a better potential energy hypersurface, at CCSD(T) quality or better. Further research could also explicitly include the umbrella bending mode.