## 1 Introduction

The fast-developing velocity map imaging (VMI) technique is increasingly used in molecular beam experiments on gas-phase polyatomic reactions. VMI is particularly interesting for reactions leading to two molecules [1–6], for it makes the accurate measurement of the angle-velocity distribution for a given quantum state of one of the two products possible. This pair-correlated angle-velocity distribution allows us to probe the dynamics of polyatomic processes at an amazing level of details, generally possible only for triatomic reactions before the invention of VMI. Consequently, this distribution is emerging as one of the most important experimental data on gas-phase chemical reactivity. In the present work, we mainly focus our developments on the theoretical description of the pair-correlated velocity distribution (PCSD), deduced from the previous density by integration over the angle. It should be noted that at the moment, this distribution appears to have been much more often measured than theoretically predicted [7–11]. There is thus an important gap to be filled. We also briefly approach the more complex issue of the calculation of the pair-correlated angular distribution.

We shall perform our study within the framework of the reactions OH+CH_{4}→H_{2}O+CH_{3} and OH+CD_{4}→HOD+CD_{3}, studied ten years ago by Zhang et al. [12–14] using the previous techniques. These seven-atom reactions leading to two polyatomic products are among the largest bimolecular processes considered to date in this type of experiment. By varying both the collision energy and the probed quantum state of the CH_{3}/CD_{3} radical, they could obtain accurate informations on the subtle way the vibrational state distribution of the co-fragment H_{2}O/HOD and the angular distribution may evolve.

If a qualitative interpretation of such results is possible from simple physical and chemical arguments, their quantitative analysis requires the development of theoretical models capable of accurately reproducing them. These models can then be used for detailed understanding of the dynamics. They can also represent (inexpensive) alternative approaches to predict data needed by atmospherical chemists and astrophysicists in order to model planetary atmospheres [15] and interstellar clouds [16], respectively.

The previous reactions take place in the electronic ground state. The first step of their accurate theoretical description requires performing ab-initio electronic structure calculations and using fitting procedures to construct the corresponding potential energy surface (PES) [17,18]. This work has recently been performed by two of us [19]. The surface obtained was called PES-2014. Note that the latter takes into account the full dimensionality of the system, i.e. 15 (reducing the dimensionality generally makes far easier the construction of PESs, but this may induce artificial dynamical behaviors ; treating a molecular collision in its full-dimensionality is often preferable). The second step of the simulation consists in moving nuclei in the previous PES. Ideally, one would like to do it quantum mechanically [20,21], and exactly. For the title reactions, however, the basis sizes necessary to converge quantum scattering calculations are huge, thus leading to exceedingly large computation times. The alternative is to classically move nuclei, a method traditionally called the quasi-classical trajectory method (QCTM) [22,23]. However, it is known that in the case where only a few vibrational states are available to the products, as is often the case when VMI is used, velocity distributions should preferably be estimated only from those classical trajectories ending in the products with integer vibrational actions, i.e. satisfying Bohr's quantization principle in the separated products [24]. During the last decade, the QCTM in this old quantum spirit proved to lead to very satisfying predictions as compared to exact quantum scattering calculations or experimental measurements [24]. This approach has been recently used by two of us [19,25] in order to reproduce some of the measurements of Liu and co-workers. In the present work, we focus on those obtained for a collision energy of ∼10 kcal/mol published in Refs. [12] and [14].

The corresponding QCTM results semi-quantitatively reproduced experimental evidence [19,25], which is already quite satisfying given the size of the processes under scrutiny. In particular, the theoretical vibrational state distributions of the co-fragment H_{2}O/HOD were found to be in good agreement with the experimental ones [25]. On the other hand, the predicted speed distributions associated with given vibrational states of the co-fragment tend to be broader than their measured counterparts and not always well located [25]. As a consequence, the PCSD obtained by summing the previous densities appears to partially or totally loose the structures seen in the experimental PCSD [12,14] (see Section 2.3). Without any further calculations, one might attribute these discrepancies to a possible lack of accuracy of the PES-2014, but we shall see that this is not the case. The goal of the present note is indeed to show that the previous disagreement can be significantly reduced by adding a few more constraints to the analysis of QCTM results so as to better comply with the restrictions imposed by VMI. In Section 2, these changes are detailed within the framework of the OH+CD_{4} reaction, and comparison between the predictions of the new analysis and measurements is made for both title processes. We end up this section discussing the high numerical cost of the pair-correlated angular distribution (PCAD), as compared to the PCSD. We conclude in Section 3.

## 2 Calculations

### 2.1 Trajectories

We used the VENUS code [26] to run 10^{6} trajectories on the PES-2014 for each process. Full computational details on these calculations can be found in refs. [19,25] (note that in these previous studies, only 10^{5} trajectories were run for each reaction). As far as OH+CD_{4}→HOD+CD_{3} is concerned, 90,646 trajectories turned out to be reactive. The product quantities involved in the next developments are the vibrational energy ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{v}\mathrm{i}\mathrm{b}}$ of CD_{3}, its rotational energy ${E}_{C{D}_{3}}^{rot}$, the vibrational energy ${E}_{\mathrm{H}\mathrm{O}\mathrm{D}}^{\mathrm{v}\mathrm{i}\mathrm{b}}$ of HOD, its rotational energy ${E}_{\mathrm{H}\mathrm{O}\mathrm{D}}^{\mathrm{r}\mathrm{o}\mathrm{t}}$, the vibrational actions a_{1}, a_{2} and a_{3} associated with the symmetric stretching, antisymmetric stretching and bending vibrational modes of HOD, respectively, the relative translational energy E_{trans} between HOD and CD_{3}, the speeds (with algebraic sign) v of CD_{3} and v′ of HOD within the center-of-mass system, and the scattering angle θ. The VENUS code provides directly all these quantities except the actions and the velocities. The actions are obtained by means of the normal-mode analysis (NMA) method [27], which includes anharmonicity and Coriolis-coupling terms. The velocities are deduced from E_{trans} as follows; the total linear momentum within the center-of-mass system being zero, we have

$${m}_{\mathrm{C}{\mathrm{D}}_{3}}v+{m}_{\mathrm{H}\mathrm{O}\mathrm{D}}{v}^{\prime}=0\text{,}$$ | (1) |

_{3}and m

_{HOD}is the analogous mass for HOD. Moreover,

$${E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}=\frac{1}{2}{m}_{\mathrm{C}{\mathrm{D}}_{3}}{v}^{2}+\frac{1}{2}{m}_{\mathrm{H}\mathrm{O}\mathrm{D}}{v}^{{\prime}^{2}}\text{.}$$ | (2) |

From the two previous equations, one obtains

$$v={\left[\frac{2{E}_{\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{s}}}{{m}_{\mathrm{C}{\mathrm{D}}_{3}}\left(1+\frac{{m}_{\mathrm{C}{\mathrm{D}}_{3}}}{{m}_{\mathrm{H}\mathrm{O}\mathrm{D}}}\right)}\right]}^{1/2}\text{.}$$ | (3) |

This velocity is the one measured by Zhang et al. using the VMI techniques [12–14].

### 2.2 Vibrational state populations

The vibrational state of HOD is denoted (n_{1}, n_{2}, n_{3}). These quantum numbers are associated with the same vibrational modes as the previous actions (a_{1}, a_{2}, a_{3}). The corresponding frequencies are ω_{1}, ω_{2}, and ω_{3} (their values are given in Ref. [25]). The population of state (n_{1}, n_{2}, n_{3}) is denoted as ${p}_{{n}_{1}{n}_{2}{n}_{3}}$. In the experiments of Zhang et al. [12–14], frozen reagents meet with a collision energy of ∼10 kcal/mol (among others) and P(v) is measured for CD_{3} in its vibrational ground state and its lowest rotational levels (N ≤ 5) [3]. The rotational energy ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{r}\mathrm{o}\mathrm{t}}$ of the probed CD_{3} is thus lower than ∼1 kcal/mol. Three vibrational states are found to be available to HOD (at a collision energy of 10 kcal/mol, to be exact [12]). These states are (1,0,0), (1,0,1) and (2,0,0), by increasing the order of energy.

We shall now estimate their populations from the reactive trajectory results. Since these results are in much larger amount than in Refs. [19,25], the present method will better cope with the constraints imposed by the VMI techniques than before. These constraints will be taken into account as follows:

- (i) the zero point energy (ZPE) of CD
_{3}is ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{0}=13.8$ kcal/mol. We assume that all the trajectories violating the ZPE, i.e. such that ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{v}\mathrm{i}\mathrm{b}}<13.8$, do not contribute to the reactivity. We thus discard them from the statistics. - (ii) the energy of the first excited vibrational state of CD
_{3}is 14.9 kcal/mol [25], corresponding to one quantum in the lowest umbrella mode of CD_{3}(396 cm^{−1}or 1.1 kcal/mol) above the ZPE. Calling δ the previous quantum, the intermediate energy between the two previous energies is ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{0}+\delta /2=14.35\phantom{\rule{0.25em}{0ex}}\text{kcal/mol}$. We now assume, in the spirit of standard binning (SB) [24], that all the trajectories leading to ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{v}\mathrm{i}\mathrm{b}}>14.35\phantom{\rule{0.25em}{0ex}}\text{kcal/mol}$ do not contribute to the vibrational ground state of CD_{3}probed by Zhang et al. [12–14], and should thus be ignored. - (iii) finally, we also reject all the trajectories leading to ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{r}\mathrm{o}\mathrm{t}}>1\phantom{\rule{0.25em}{0ex}}\text{kcal/mol}$.

In Fig. 1, the cloud of points of coordinates $\left({E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{v}\mathrm{i}\mathrm{b}},{E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{r}\mathrm{o}\mathrm{t}}\right)$ is represented for the 90,646 reactive trajectories. Gray/red points correspond to rejected/accepted trajectories. Only 722 paths turn out to be useful, thus representing only ∼0.07% of the whole set of calculated trajectories. This clearly illustrates the fact that simulating polyatomic reactions studied by VMI is numerically demanding, even with the QCTM.

A purely classical treatment of vibrational motions is realistic when many vibrational states (ideally, several tens) are available to the products. In the present case, however, only three vibrational states of HOD are found to be populated [12] and it is preferable, as stated in the introduction, to take into account Bohr's quantization principle in the analysis of the final trajectory results [24]. For polyatomic molecules, this is usually done by means of the energy-based Gaussian binning procedure, often called the 1GB procedure [24,25,28–31]. The latter puts strong emphasis on those trajectories leading to HOD with its classical vibrational energy close to one of the quantum energies of the available vibrational states. In practice, the 1GB estimation of ${p}_{{n}_{1}{n}_{2}{n}_{3}}$ is performed as follows: we consider the set of trajectories for which ${\overline{a}}_{1}={n}_{1},\phantom{\rule{0.25em}{0ex}}{\overline{a}}_{2}={n}_{2}\phantom{\rule{0.25em}{0ex}}\text{and}\phantom{\rule{0.25em}{0ex}}{\overline{a}}_{3}={n}_{3}$, where ${\overline{a}}_{i}$ is the nearest integer of a_{i}, i = 1, 2 or 3. In other words, trajectories for which (a_{1}, a_{2}, a_{3}) does not belong to the unit cube centered at (n_{1}, n_{2}, n_{3}) are rejected. This selection is common to both the 1GB and SB procedures. However, contrary to SB which assigns the democratic unit weight to each trajectory, the 1GB weight, within the harmonic treatment of HOD vibrational motions, is given by

$${p}_{{n}_{1}{n}_{2}{n}_{3}}^{1\mathrm{G}\mathrm{B}}\left({a}_{1},{a}_{2},{a}_{3}\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\left[\frac{\sum _{i=1}^{3}{\omega}_{i}\left({a}_{i}-{n}_{i}\right)}{\mathit{\varepsilon}\sum _{i=1}^{3}{\omega}_{i}}\right]}^{2}\right)$$ | (4) |

^{−1}above the ZPE of HOD, respectively [32]. Note that the set of trajectories used in the anharmonic case is still the one for which (a

_{1}, a

_{2}, a

_{3}) belongs to the unit cube centered at (n

_{1}, n

_{2}, n

_{3}). ${p}_{{n}_{1}{n}_{2}{n}_{3}}$ is obtained by summing the 1GB weights over the previous paths and normalizing to unity the set of available populations. The full width at half maximum of the Gaussian in Eq. (4) $\left(\text{in}\phantom{\rule{0.25em}{0ex}}\text{terms}\phantom{\rule{0.25em}{0ex}}\text{of}\phantom{\rule{0.25em}{0ex}}\text{the}\phantom{\rule{0.25em}{0ex}}\text{energy}\phantom{\rule{0.25em}{0ex}}\text{difference}\phantom{\rule{0.25em}{0ex}}\sum _{i=1}^{3}{\omega}_{i}\left({a}_{i}-{n}_{i}\right)\right)$ is proportional to ɛ. In principle, strictly applying Bohr's quantization would thus require taking ɛ as close to 0 as possible [24,29], so the trajectories contributing to ${p}_{{n}_{1}{n}_{2}{n}_{3}}$ would be those leading to classical vibrational energies nearly equal to the quantum energy of state (n

_{1}, n

_{2}, n

_{3}). In the following, we identify these paths as those for which ${p}_{{n}_{1}{n}_{2}{n}_{3}}^{1\mathrm{G}\mathrm{B}}\left({a}_{1},{a}_{2},{a}_{3}\right)$ is larger than 0.5 (see Eq. (4)) and call them 1GB paths. Their number, among the 722 paths complying with constraints (i–iii), is represented in Fig. 2 as a function of ɛ (blue curve). At ɛ = 0.041, the number of 1GB trajectories is equal to 100, which can be roughly considered as the minimum statistical ensemble for a reasonable estimation of the three 1GB populations. The latter are represented in terms of ɛ in Fig. 3. The vertical red dashed line is defined by ɛ = 0.041. On the right side of this line, we note that the populations vary almost linearly with ɛ. On the left side, the smaller ɛ, the smaller the statistics and the more irregular the variations of the populations. Since the experimental populations are assumed to be the best described by the 1GB populations corresponding to ɛ = 0, and the only satisfying 1GB populations are for ɛ ≥ 0.041, we deduce the former by linearly extrapolating the latter from ɛ = 0.041 down to 0. Approximating the 1GB populations by straight lines passing through the points of the curves corresponding to ɛ = 0.041 and 0.1 (see Fig. 3), we arrive at p

_{100}= 23, p

_{101}= 12 and p

_{200}= 65% (rounded to their nearest integer), to be compared with the experimental ones, i.e. p

_{100}= 19, p

_{101}= 11 and p

_{200}= 70% [12,33]. Note that at ɛ = 0.1, the number of 1GB trajectories is equal to 218.

Constraints (i) and (ii) define an energy-based windowing, which is related somehow to SB and 1GB [24,28–31]. The main difference is that the former procedure takes into account paths for which the nearest integers of the vibrational actions of CD_{3} are non 0, contrary to the latter. This is probably a bit less physical, but significantly increases the number of trajectories contributing to the populations, in addition to avoiding the determination of the six vibrational actions of CD_{3}.

In order to improve the statistics, we also made calculations with the maximum value of ${E}_{\mathrm{C}{\mathrm{D}}_{3}}^{\mathrm{r}\mathrm{o}\mathrm{t}}$ kept at 2 kcal/mol (instead of 1 kcal/mol), but this affected the accuracy of the vibrational state populations. Moreover, this artificial shift, whenever used, should be as small as possible, since it allows for more energy channeled into the rotational motion of CD_{3}, at the expense of its translation motion. This is expected to shift the speed distribution towards small velocities.

### 2.3 Pair-correlated speed distribution

The contribution of state (n_{1}, n_{2}, n_{3}) to the PCSD P(v) is denoted as ${P}_{{n}_{1}{n}_{2}{n}_{3}}\left(v\right)$. This quantity will be normalized to unity. From the previous developments, we have

$$P\left(v\right)={p}_{100}{P}_{100}\left(v\right)+{p}_{101}{P}_{101}\left(v\right)+{p}_{200}{P}_{200}\left(v\right)\text{.}$$ | (5) |

We now concentrate on the calculation of P_{100}(v), P_{101}(v) and P_{200}(v). We use the same constraints on the rovibrational energy of CD_{3} as previously used. Like ${p}_{{n}_{1}{n}_{2}{n}_{3}}$, ${P}_{{n}_{1}{n}_{2}{n}_{3}}\left(v\right)$ should better be calculated in a quantum spirit. In order to illustrate this, we first focus on the most probable state (2,0,0) and calculate P_{200}(v) in the standard following way ignoring any quantization: we apply the SB procedure, i.e. we focus on trajectories for which (a_{1}, a_{2}, a_{3}) belongs to the unit cube centered at (2,0,0). We then divide in 80 equal intervals the range [0.5, 2.5] comprising the velocities measured in km/s [12] and count the number of trajectories per interval (the resulting curve is made of points separated by 0.025 km/s, thus leading to a satisfying resolution). The resulting SB distribution, normalized to unity, is represented by the blue curve in Fig. 4, to be compared with the magenta experimental curve [12]. The SB curve appears to be much broader than the magenta curve. We thus calculate P_{200}(v) by means of the 1GB procedure, following the same route as above, with the difference that trajectories are now assigned 1GB weights with ɛ = 0.041 (see Eq. (4)) instead of unit weights. The resulting distribution is given by the green curve in Fig. 4. The agreement with the experiment appears to be much more satisfying than for the SB procedure. Moreover, we repeat this calculation with ɛ = 0.022 (not shown) and observe that despite more irregularities in P_{200}(v) due to worse statistics, its width is nearly unchanged. The distribution at ɛ = 0.041 is thus expected to be close to its extrapolation at ɛ = 0. The analogous distributions for the states (1,0,1) and (1,0,0) are displayed in Figs. 5 and 6. As a matter of fact, the 1GB distributions are slightly shifted towards the large velocities, revealing some inaccuracies in the theoretical treatment. The three distributions, weighted by their corresponding 1GB vibrational state populations, are represented in Fig. 7 by the blue, green and orange curves for the (2,0,0), (1,0,1) and (1,0,0) states, respectively. Summing them up leads to P(v) (see Eq. (5)).

At this point, a comment on the experimental measurement of P(v) is in order. In a hypothetical experiment where both reagents would be in the rovibrational ground state and the collision energy would be perfectly controlled, the total energy available to the products would take the same value for the whole set of reactions detected. Since E_{trans} is equal to the previous energy minus the quantized internal energy of the products, E_{trans} would also be quantized. Consequently, P(v) would be given by a set of Dirac peaks with given weights. In our calculations, however, we pseudo-quantize the vibration motion through the 1GB procedure, but not the rotation motion. P(v) is thus a rough approximation of what would be observed if the measurement of the final velocities were introducing a blurring of the peaks slightly exceeding the average rotational spacing between them. In real experiments, however, the blurring is much larger, for the collision energy and the reagent rotational states are not fully controlled and VMI itself significantly contributes to the blurring. We thus take this into account in the calculation of the PCSD through the Gaussian convolution

$${P}_{c}\left(v\right)=\int \mathrm{d}{v}^{\prime}\frac{1}{\sqrt{\pi}\eta}exp\left[-{\left(\frac{{v}^{\prime}-v}{\eta}\right)}^{2}\right]P\left({v}^{\prime}\right)\text{.}$$ | (6) |

The value of η is chosen so as to reproduce as satisfyingly as possible the experimental speed distribution at the threshold and cut-off. With η kept at 0.08, one obtains the red curve in Fig. 7, to be compared with the experimental one (magenta curve), extracted from Fig. 2a of Ref. [12]. The convolution is welcome, as a million trajectories is not enough to get a smooth P(v). Without blurring, good agreement between experiment and theory would likely be beyond reach. Compared with the previous predictions, displayed in the lower left panel in Fig. 2 of Ref. [25], the PCSD is clearly improved. In particular, when summing up the three bells in the previous figure, only two bumps are found, instead of the three bumps observed both in the experimental and present theoretical results.

As far as the OH+CH_{4} reaction is concerned, five states of H_{2} O are found to be available (at a collision energy of 10.2 kcal/mol, to be more precise [14]). These states are (0,0,1), (0,1,0), (1,0,0), (0,1,1) and (1,0,1), by increasing the order of energy. We checked that the harmonic normal mode description of water vibrations is very satisfying for the energies of the previous states. For the PES-2014, the vibrational frequencies of the OH stretching, OH^{∗} stretching and HOH bending are ω_{1} = 3806, ω_{2} = 3762 and ω_{3} = 1651 cm^{−1}, respectively. Moreover, the ZPE of CH_{3} is 18.8 kcal/mol, and the lowest umbrella frequency is 510 cm^{−1} or 1.5 kcal/mol. The same type of calculation as before has been performed. Over 10^{6} trajectories run on the PES-2014 [19,25], 123567 are found to be reactive, and 1013 comply with the three constraints imposed by VMI. The dependence of the 1GB populations on ɛ is displayed in Fig. 8. The region of statistical confidence roughly lies on the right side of the red dashed line, defined by ɛ = 0.053. By linearly extrapolating the 1GB populations from ɛ = 0.053 down to 0, we arrive at p_{001} = 6, p_{010} = 53, p_{100} = 19, p_{011} = 10 and p_{101} = 12%. The speed distributions P_{001}(v), P_{010}(v), P_{100}(v), P_{011}(v) and P_{101}(v) weighted by the previous populations are represented in Fig. 9. With η kept at 0.09, Gaussian convolution of the sum of these densities leads to the PCSD represented by the red curve in Fig. 9, to be compared with the experimental magenta curve extracted from Fig. 2a of Ref. [14]. Note that states (1,0,1) and (0,1,1) have similar energies, so the speed distributions associated with these states overlap (see the blue curves in Fig. 9). This comment holds for states (1,0,0) and (0,1,0) (see the green curves in Fig. 9). This is the reason why the vibrational populations deduced from the speed distributions were previously denoted ${p}_{{01}^{\text{*}}1}$ for states (1,0,1) and (0,1,1), ${p}_{{01}^{\text{*}}0}$ for states (1,0,0) and (0,1,0) and p_{001} [25]. From the present calculations, ${p}_{{01}^{\text{*}}1}=22$, ${p}_{{01}^{\text{*}}0}=72$ and p_{001} = 6 (see above), to be compared with the experimental populations, i.e. ${p}_{{01}^{\text{*}}1}=19$, ${p}_{{01}^{\text{*}}0}=75$ and p_{001} = 6 [14]. Like previously observed, the PCSD is clearly improved as compared with its previous estimation, displayed in the upper left panel in Fig. 2 of Ref. [25]. In particular, when summing up the three bells in the previous figure, a single bell is found, instead of the three bumps observed both in the experimental and present theoretical results.

Last but not least, we wish to emphasize that standard QCT-SB calculations do not even allow us to reveal the structures seen in the experimental PCSDs, as can be seen in Fig. 10. The QCT-SB distributions are indeed simple bell-shape curves broader than the experimental distributions. In comparison, QCT-1GB distributions are much more realistic (see Figs. 7 and 9).

### 2.4 High numerical cost of the pair-correlated angular distribution

Following Ref. [19], we estimate the PCAD by limiting the restrictions imposed by VMI to the second constraint in Section 2.2 and using the 1GB procedure with ɛ = 0.041. This distribution is called P_{a}(θ). For the OH+CD_{4} reaction, the number of 1GB trajectories is then equal to 4417. P_{a}(θ) is found to nearly reproduce the previous prediction represented by the solid black curve in Fig. 7 of Ref. [19]. This distribution turns out to be in qualitative agreement with the experimental PCAD of Zhang et al. represented by the dashed black curve in the same figure (see also Fig. 3a in Ref. [12]). The main difference is that the former involves a forward contribution not seen in the latter. One may thus wonder whether calculating the PCAD in a way consistent with the estimation of the PCSD might make this forward contribution disappear. Moreover, it is logical to use the same approach for estimating different observations measured under the same conditions. This calculation would require adding constraints (i) and (iii) (see Section 2.2) to the method leading to P_{a}(θ) and deduce the PCAD at ɛ = 0, called P_{b}(θ), from its values at ɛ = 0.1 and 0.041 by linear extrapolation. However, the PCAD is determined by means of 218 and 100 trajectories at ɛ = 0.1 and 0.041, respectively. Now, these numbers are too small for the PCAD to be reliable, especially in the neighborhood of the poles, as is clearly shown in the study by Truhlar and Blais [34]. In particular, it is shown in Figs. 7 and 8 of the previous work that at least one or two thousand trajectories are necessary to get significant results around the poles. The accurate calculation of P_{b}(θ) for the title processes would thus imply running at least 10 million trajectories.

## 3 Conclusion

The development of the velocity map imaging (VMI) technique [1] over the last twenty years [2] has made possible the study of polyatomic reactions at an unprecedent level of details through the measurement of pair-correlated speed distributions. For the title processes, for instance, one may extract from the previous observations accurate data on the subtle way H_{2}O/HOD vibrates and rotates for a given quantum state (or set of closely spaced quantum states) of CH_{3}/CD_{3} [12–14]. If one is able to accurately reproduce these data by an adequate theoretical treatment of nuclear motions on a high-quality ab-initio potential energy surface, one may then use this treatment to get a deep understanding of the dynamics of the processes under scrutiny and improve the synergy between experiment and theory. We have confirmed in this work that the QCT method can be such an approach provided that Bohr's quantization of product vibration motions is taken into account in the final analysis of trajectory results [24,25,28–30]. This is an encouraging result considering (a) the scarcity of theoretical predictions on pair-correlated speed distributions, and (b) the computational challenges of quantum scattering calculations for processes involving large systems. Moreover, we have confirmed that the PES-2014 correctly describes the dynamics of the title processes. We plan to apply this QCT method in a quantum spirit to several polyatomic bimolecular reactions and photodissociations for which pair-correlated speed distributions have never been theoretically reproduced. It should be noted that pair-correlated angular distributions can also be predicted by this method, though the number of trajectories required should be much larger than for the previous distributions.

## Acknowledgments

LB is grateful to Prof. Kopin Liu for clarifying explanations on the cross-beam experiments performed in his laboratory.