The isotopic exchange reactions that can occur during the collision between an oxygen atom and an oxygen molecule involve, as the intermediate, the metastable ozone O 3 * in excited rovibrational states above the dissociation threshold:
However, an adequate modeling of this process represented significant difficulties for theory [1][2][3][4][5][6][7][8][9][10][11][12][13] for several decades since the related experimental measurements [14][15][16][17] were published. The probability of this reaction depends on the properties of excited ozone O 3 *, which depend on the potential energy surface (PES) supporting the dynamical process.
The ozone PES has been the subject of many ab initio studies over the years (refs 18-27 and references therein) that revealed a complicated electronic structure 20,23,[28][29][30] and rovibrational patterns [31][32][33][34] of this seemingly simple triatomic molecule. Extensive spectroscopic works (refs 35-40 and references therein) have been carried out both for remote sensing applications in the atmosphere and for the validation of ab initio predictions in a large spectral range from the fundamental bands up to high overtones and combination rovibrational transitions toward the dissociation threshold. 31,41,42 A description of the ozone formation is among the main incentives for dynamical studies, particularly in the upper atmosphere where the so-called "nascent population" 43 of highly excited vibrational ozone states is not yet well-known. This information is mandatory for a correct interpretation 44,45 of ozone measurements by satellite instruments at conditions of nonlocal thermodynamic equilibrium.
In most studies, it is believed that ozone formation under the stratospheric low-pressure conditions proceeds via a threebody recombination process. According to the Lindemann mechanism, the first step in the recombination, leading to the formation of an excited complex O + O 2 → O 3 *, which is the same as in eq 1. The second step is a subsequent stabilization by a collision with another partner O 3 * + M → O 3 + M′, which absorbs an excess of the kinetic energy and makes the metastable complex O 3 * "falling down" into the potential well. Meanwhile, the group of Troe 46,47 has introduced the radical complex ("chaperon") mechanism, 48 which may also play an important role at higher pressure.
Since the experimental discovery of surprising isotope anomalies in the ozone formation, 49,50 many efforts have been devoted to a theoretical modeling of these processes, but a full understanding of related "strange and unconventional" effects 9,[51][52][53][54][55] is still lacking. This involves an unusually high enrichment of heavy isotopomers in the stratosphere and the mass-independent fractionation (MIF) of oxygen isotopes in the ozone formation, which was considered as a milestone in the study of isotope effects. 55 The detailed reviews of the problem can be found in status reports by Schinke et al., 9 Marcus, 55 and Thiemens. 56 A large symmetry selection in the formation of the ozone molecule required the introduction of ad hoc factors 51,55 to fit observed deviations from a simple statistical behavior. Recent experiments and their modeling 57 have shown that purely statistical theories cannot explain the isotopic effects, which therefore require dynamic state-specific studies. A series of theoretical work performed by the group of Babikov [58][59][60][61] explored the role of different types of resonances formed in collisions of O and O 2 , and showed that the isotope effects could be linked to the formation of Feshbach-type resonances.
The isotopic exchange reaction (1) is in a competition with the ozone formation process and exhibits strong isotopic effects 9,11,12,[14][15][16][17]62 as well. Because it is only a two-body process, the corresponding dynamics is easier to investigate rigorously. Since both reactions proceed on the same PES, modeling the reaction of eq 2 in the time domain would help in understanding both the dynamics of the isotopic exchange and the ozone recombination, and also would be a step forward in the interpretation of the MIF effects.
In this work, we focus on qualitative state-specific features in wave propagation dynamics in the collisional process (2) related to the formation of the metastable resonances of the excited ozone O 3 *. The study aims at better understanding of the impact of the accuracy of the PES in the transition state region, in particular, on the role of a hypothetical "reef" structure on the minimum energy path (MEP) in O + O 2 collisions.
Since many previous studies, [8][9][10][11][12][13]62 it has been recognized that an accurate potential energy surface is a prerequisite for quantum dynamical studies of these processes. A perfect PES must be also applicable for an interpretation of spectroscopic experiments that allow reliable validation of ab initio calculations. It should be able to predict the dissociation threshold and the vibrational bands correctly, and it should have a physically meaningful shape at the transition state (TS) on the way between the molecule and the fragments. Earlier ab initio calculations 18 predicted an activation barrier on the MEP at the TS geometries between r = 3 and r = 5 bohr, where r is the O-O bond distance.
The first accurate global three-dimensional (3D) ab initio PES of the ozone molecule has been constructed by Siebert et al. 19,20 in the group of Schinke, which is usually referred to as the SSB PES. 19 This PES was calculated using the multireference configuration interaction internally contracted (ic-MRCI) method 64 with the Dunning atomic basis set corresponding to the quadruple cardinal number (l = 4) augmented by diffuse functions. The SSB PES predicted quite good low energy frequencies (Table 1), but it has a significant drawback: a strong underestimation of the dissociation energy (O 2 ( 3 Σ g -) + O( 3 P)) D e by about -11%. With respect to the most accurate experimental bond dissociation (O 3 → O 2 + O) energy D 0 (obs) = 8485 cm -1 calculated by Ruscic et al. 63,65 (as cited in Holka et al. 22 ), the D 0 (SSB) was lower by about 900 cm -1 . This disagreement was too large for an interpretation of highly sensitive laser absorption experiments 31,37,41 extended up to nearly 93% of D 0 . With the SSB value for D 0 , this range would fall in the vibrational continuum of the SSB surface, which is clearly not the case. As in the previous ab initio studies, the SSB PES exhibits a small barrier for the bond O-O length around 4 bohrs of about 50 cm -1 above the D e asymptotic energy.
One-dimensional (1D) and two-dimensional (2D) calculations have shown 7,66 that the barrier would be converted to a "submerged reef" below D e if larger quintuple (l = 5) or sextuple (l = 6) atomic basis sets were applied. This reef was then followed by a shallow van der Waals (vdW) minimum at larger O-O distances of about 4-5 bohrs. An empirically modified mSSB 3D PES was constructed 21 with the "reef" height E reef = -114 cm -1 with respect to the D e asymptotic energy.
The SSB and mSSB PESs were employed by many groups for the theoretical modeling of the isotope exchange dynamics, [6][7][8][9]57,67 using a quantum statistical model, quasiclassical trajectories, and wave propagation method, but it did not produce a satisfactory agreement with observations. [14][15][16][17]57 On the other hand, it was shown that a complete basis set (CBS) limit (either l = 4, l = 5 → ∞, 16,22,25 or l = 5, l = 6 → ∞, 22,26 or l = 3, l = 4 → ∞ 23 ) was necessary to approach the experimental D 0 value. 63,65 At this point, the situation was quite confusing because at the same time the increasing cardinal number of the atomic basis set beyond l = 5 conducted to stretching vibrational ν 1 and ν 3 frequencies, which are significantly larger than the experimental values (Figure 2 in Tyuterev et al. 26 ). For example, the Ayouz-Babikov (AB) PES 25 computed with the l = 4, l = 5 →∞ CBS limit had a correct D 0 value exhibiting nearly the same "reef" shape at the TS as the mSSB. 20 However, it overshoots the fundamental ν 3 stretching band center by 31 cm -1 with increasing error for the overtone bands, being much less accurate for the spectroscopy than the SSB PES (Table 1).
An important question was whether a reef barrier appeared as an artifact of ab initio approximations presumably caused by an avoided crossing with an excited electronic states. 23 Muller et al. have shown (as cited in ref 26) that the reef would disappear on the 1D MEP if noncontracted MRCI calculations were implemented using the Columbus software. 68 However, a corresponding ansatz was too demanding for a construction of a full 3D PES. Dawes et al. 23 have investigated an impact of the inclusion of several electronic states in the ic-MRCI procedure on the TS range for the ground state. Their first work considered 1D and 2D PES corrections using the l = 3, l = 4 → ∞ CBS extrapolation that changed the shape of the MEP. Dawes et al. 23 have found that due to an optimization of orbitals including several electronic states the reef was transformed to a kind of smooth shoulder. Preliminary estimations using this result together with the approximate quantum statistical model 23 gave possible hints to obtain more consistent temperature dependence of the isotope exchange rates.
Another question concerning the 3D ozone PES was how to combine both a reasonably good D e and reasonably good vibrations. To this end, two full dimensional ab initio surfaces were constructed by Tyuterev et al., 26 hereafter referred to as TKTHS PESs. One of the versions has been constructed via hybrid ic-MRCI calculations by sewed results from different atomic basis sets. Within the main C 2v well up to the TS, the quintuple atomic basis set (with the cardinal number l = 5) was used to ensure reasonably good fundamental vibration frequencies. For the energies above the TS the l = 5, l = 6 → ∞ CBS limit was used to ensure a correct D e asymptotic energy near the experimental value (Table 1). This first PES was called "R_PES" because it maintained a reef structure, but a much less pronounced one: 3 times lower barrier height E reef = -320 cm -1 with respect to the D e than in mSSB 20 or in AB PESs, which had E reef = -114 and -102 cm -1 , respectively.
Another TKTHS PES was constructed using a conceptually more consistent ab initio approach. [The FORTRAN subroutine for the TKTHS PES 26 in the symmetry adapted coordinates is available at the S&MPO information system dedicated to ozone spectroscopy via the link http://smpo.univ-reims.fr/files/codes/.] The CBS l = 5, l = 6 →∞ limit was employed in the entire range of nuclear geometries from the C 2v well to the dissociation accounting for 2D "Dawes corrections" 23 on the stretching manifold due to the effect of excited electronic states. The TKTHS ab initio surface constructed in this way was called "NR_PES" (no-reef) because it did not exhibit the reef structure. A comparison of the MEP cuts for these R_PES and NR_PES potentials in the range from the dissociation asymptotic energy to the transition state is shown in Figure 1. The TKTHS 26 NR_PES provides currently the most accurate ab initio predictions for highresolution spectra assignments: an average error of vibrational bands was about 1 cm -1 for various ozone isotopologues 13,31,[39][40][41][42] measured up to 93% of the dissociation threshold. The validity of the corresponding ab initio ansatz has been recently confirmed by an excellent agreement of predicted band intensities with experiment. 69 Later on, Dawes et al. 27 published a full 3D PES (usually referred to as DLLJG PES) using the ansatz of ref 23 and including the long-range interaction part from Lepers et al. 24 and the spin-orbit correction. The DLLJG PES has been employed for the modeling of reaction 2 by Li et al. 10 and by Sun et al. 11 using the time-dependent wave packet method. This permitted for the first time reproduction of a negative temperature dependence of the rate constants, though their predicted absolute values were still underestimated with respect to the observations. 7,14,15,17 These results on the DLLJG PES have been confirmed by Rajagopala Rao et al. 12 using a time-independent approach.
The reader can find a comparison of the MEP cut for various ab initio ozone PESs in Figure 1 of the recent work by Guillon et al. 13 The TKTHS NR_PES 26 was more attractive in the TS range than all other published PESs including that of Dawes et al. 27 Guillon et al. 13 have shown that this NR_PES 26 enabled one to obtain for the first time an excellent quantitative agreement of the theoretical rate coefficient K 866 with observations. [14][15][16][17] Using full quantum mechanical timeindependent dynamics calculations, Honvault et al. 62 have published similar work for the theoretical rate coefficient K 688 , but a time-dependent study on this PES, which proved to be accurate both for spectroscopy and for rate constants, has not been yet reported.
In this work, we show that the time evolution of the wave packet significantly depends even on small topographic features situated relatively deep in the potential well below the dissociation threshold. To this end, we compare timedependent wave functions obtained using two TKTHS PESs (NR_PES and R_PES described above), which are very similar in the bottom of the C 2v well and at the D e asymptotic energy. A comparison of the high-resolution ozone spectra analyses 31 indicated a clear preference to the NR_PES with respect to the R_PES, but the collisional dynamics on the latter one has not been studied up to now. Minimum energy pathways of the potential energy surface for R_PES 26 (red lines) and NR_PES 26 (black lines) in Jacobi coordinates.
A previous work of Guillon et al. 13 has shown that the TKTHS NR_PES 26 provides an excellent agreement with the experimental isotope exchange rate (2) using a computationally demanding state-averaging procedure within the timeindependent formalism. Here we consider time evolution of the nuclei-exchange process involving a formation of metastable state patterns for the excited ozone O 3 * during the scattering reaction (2) starting from the initial state (ν = 0, j = 1) state of 32 O 2 . This brings a supplementary insight in the corresponding state-specific process.
Figure 1 shows the MEP of the NR_PES and the R_PES from Tyuterev et al. 26 along R, where R is the distance between the 18 O atom and the center of mass of the 32 O 2 molecule. The internuclear distance between two 16 O atoms and the Jacobi angle in the body fixed frame are denoted as r and θ. Along the MEP, θ takes different values in the circled regions in Figure 1, while r is at the equilibrium distance of the 32 O 2 molecule. One can see that relatively small differences between the values of the NR_PES and R_PES in the reef region are only within 50 cm -1 . The barrier reef height for the R_PES lies much deeper inside the well (E reef = -320 cm -1 ) in comparison with other PESs, 20,21,25 which exhibited either a true barrier or a reef structure and had been used for the dynamics calculations in the previous works. [6][7][8][9]57,67 For the mSSB or AB surfaces the reef was only about -100 cm -1 below D e (see Table 1). Another important difference is that the R_PES is more attractive on the ozone MEP than the NR_PES just above the TS range (Figure 1), whereas all other published PESs were significantly less attractive than both of them (see Figure 1 of Guillon et al. 13 ). Despite the presence of a deeply submerged reef, it is unlikely that the R_PES would support the dynamics in the same way as was reported in earlier works. [6][7][8][9]57,67 To study the time evolution of the wave packet, the multiconfiguration time-dependent Hartree [70][71][72][73] (MCTDH) method is used. It is an algorithm to solve the time-dependent Schrodinger equation, which can be considered as a timedependent version of the multiconfigurational self-consistent field (MCSCF) method applied to the nuclei. Within this method the wave function Ψ(Q, t) of the system is written as a sum of products of single-particle f unctions (SPFs), forming a time-dependent orthonormal basis set. SPFs are low-dimensional functions: When they contain more than one degree of freedom (DOF), the combined coordinates Q κ ≡ q 1,κ , ..., q d κ ,κ that comprises d κ physical DOFs are introduced.
The ansatz of the MCTDH wave function reads
where f and p denote the number of degrees of freedom and number of combined modes respectively of the system, A J ≡ A j 1 ,...,j p denotes the MCTDH expansion coefficients, and the configuration or Hartree products Φ J are products of SPFs defined in relation 3. The SPFs are finally represented by linear combinations of time independent primitive basis functions as
( ) ,
where the primitive functions χ l i (κ) (q i,κ ) are usually discrete variable representation (DVR) 74,75 functions with the time-
. The equation of motion for the coefficients and the SPFs are then obtained from the Dirac-Frenkel variational principle. The strength of the method relies on the fact that the number of optimized SPFs used for the propagation is usually much smaller than the number of functions of the primitive basis.
In our particular case, the modes are simply the three nuclear degrees of freedom, R, r, and θ, and the MCTDH wave function reads
The number of grid points (N r /N R /N θ ) and SPFs (n r /n R /n θ ) are given in Table 2: They guarantee a high level of convergence.
Since the form of the original PES is not adapted to MCTDH, we have utilized the POTFIT [76][77][78] algorithm to recast the operator in a sum of products of one-dimensional functions. The original multidimensional function is assumed to be given by its values on a multidimensional product grid. POTFIT operates on the grid points only. We thus assume that the values of a potential V are given on a product grid
where {q} i κ (κ) denotes a grid point of the κth coordinate on a grid with 1 ≤ i κ ≤ N κ . Here N κ denotes the number of grid points of the grid for the κth coordinate and p is the index for coordinates. Similarly as in MCTDH, the variable {q} may be a one-dimensional or multidimensional (i.e., collective) coordinate. Next we define potential density matrices ϱ (κ) : , with components v i κ j κ (κ) as well as their corresponding eigenvalues λ j κ (κ) called natural weights. The natural weights are assumed to be in decreasing order, λ j κ (κ) ≥ λ j κ +1 (κ) , and we introduce the notation v j κ (κ) ({q} i κ (κ) ) ≔ v i κ j κ (κ) . We
then may now approximate the potential as follows:
...
(1) ( )
with m κ ≤ N κ . Here, the expansion coefficients C j 1 ...j p are the overlaps between the potential and the natural potentials
...
(1) ( )
In order to make the representation as compact as possible, a contraction over the νth mode is performed. The following functions are defined:
This allows us to rewrite the potential energy surface as
(1) (
... ... 11) and to reduce the number of expansion terms by the factor m ν . The original potential is exactly reproduced (on the grid points) when m κ = N κ . However, a sufficiently accurate approximation of the potential is usually obtained with much smaller values of m κ . POTFIT along with MCTDH has already been applied for the calculation of eigenvalues and quantum resonances at a spectroscopic accuracy 79,80 on other systems than O 3 .
Here, the contracted mode corresponds to the coordinate R and the values of m r and m θ are given in Table 2 resulting in a root-mean-square error on the grid points in the relevant (physical) region of 0.6917 meV.
For total angular momentum J = 0, the total Hamiltonian in the body-fixed frame is
where
and μ R and μ r are the reduced mass of the 18 O atom with respect to the 32 O 2 molecule and the reduced mass of the 32 O 2 molecule, respectively. In order to compute the flux going into the rearrangement channel and to limit the grid size of the system, two complex absorbing potentials (CAPs) W r and W R are introduced into the Hamiltonian, [81][82][83][84][85] such that
where
is the Heaviside step function, and β Q is the order of the CAP.
Together with the Hamiltonian of the system, once the initial wave function is set up, it will be propagated in time and the flux going into the rearrangement channel will be computed. 86 The initial wave function Ψ 0 is expressed as the product of a Gaussian wave packet of 18 O atom χ 0 (R), the vibrational wave function of the 32 O 2 molecule, and the associated Legendre polynomial, which is an eigenfunction of the operator j 2 with eigenvalue j 0 (j 0 + 1). The initial rovibrational state of the dimer is chosen to be ν 0 = 0 and j 0 = 1. The initial momentum and width of the Gaussian wave packet are chosen with the initial energy distribution desired, which is given by 87
With Ψ(t) ≡ e -iH ̃tΨ 0 , the reaction probability is given by
The parameters used in our calculation are shown in Table 2. The starting position of the CAP is adjusted such that the tail of the initial wave packet does not overlap with the CAP. The strength of the CAP is selected based on the inspection of the reflected wave packet from the CAP. Propagation time is selected based on the convergence test on the reaction probability at the selected initial wave packet.
Using eq 20, we compute the probability of reaction 2 with the total angular momentum J = 0 using the NR_PES and the R_PES. Figure 2 compares the results obtained using the NR_PES and the R_PES, as well as the results from the timeindependent calculation by Guillon et al. 13 The probability obtained with the R_PES has a different fine structure pattern compared to the NR_PES result, mainly due to different positions of resonances present in the collisional spectra for the two surfaces. Also, from the lower panel of Figure 2, one can see that the NR_PES probabilities are, in general, larger than the R_PES result.
The result by Guillon et al. is shifted slightly to adjust the rovibrational energy of the 32 O 2 molecule. We see that from a collision energy larger than 0.034 eV, the result with the NR_PES agrees well with the time-independent calculation. At lower collision energy, positions of resonances in the NR_PES probabilities agree with the time-independent result, but the magnitude is about 50% larger. This is possibly due to the influence of artificial reflections of the low-energy components of the wave packet in the time-dependent calculations. These components are partially reflected from the CAP, which is unable to absorb a wide range of wavenumbers of the outgoing flux. One can overcome this problem by extending the CAP length and the grid, and by lowering the CAP strength. However, this would increase the computational cost significantly. Since we are interested in the dynamics before the fluxes enter the CAP, such improvement is not essential for the present discussion. The agreement between the NR_PES result and the numerically accurate result by Guillon et al. at higher collision energies justifies the validity of the present time-dependent approach at higher energies.
Figure 3 displays the initial probability density of the timedependent wave function as a function of (r, R) and (θ, R) respectively. The density function is computed by integrating the wave function Ψ(r, R, θ, t) over the third coordinate. The PES in the r-R plot has the potential energy path being minimized in θ, while the PES in the θ-R plot has the potential energy path being minimized in r. The contour lines for both PESs begin at -7161 cm -1 and end at 8903 cm -1 with an increment of 2008 cm -1 . Because the difference between values of the NR_PES and R_PES in the reef region is less than 50 cm -1 , the contour lines on this scale are indistinguishable among the two PESs. The initial wave packet in the r-R plot is centered at R = 7 bohrs with a width 0.5 bohrs while being distributed along r = 2 bohrs-2.5 bohrs, as the density function is just the product of the vibrational wave function of the 32 O 2 molecule and the Gaussian wave packet 0 (R) times a constant. The density function in the θ-R plot shows similar features along R, but delocalized along the range of θ since the initial rotational quantum number j 0 of the 32 O 2 molecule is 1.
Figure 4 shows the time-dependent dynamics of the density function in coordinates (r, R) at t = 120, 240, 300, and 440 fs for the NR_PES and R_PES. At t = 120 fs, the incoming flux arrives at the C 2v well without any reflection for both PESs. At t = 240 fs, there are more fluxes in the C 2v well for the NR_PES than for the R_PES. At t = 300 fs, one can see that the incoming fluxes continue to travel to the C 2v well for the NR_PES, while significant parts of fluxes are reflected for the R_PES. At t = 440 fs, the metastable complex O 3 * in the C 2v well begins to decay toward the 16 O 18 O molecule. The reflected fluxes in the R_PES are ahead of those in the NR_PES, suggesting that the fluxes with larger kinetic energy are also reflected in the R_PES, but to a lesser extent in the NR_PES.
Figure 5 displays the dynamics of the density function in coordinates (θ, R) at the same propagation times as in Figure 4. At t = 120 fs, the reaction pathway is mainly near θ ≈ 0.9 rad, which is consistent with the MEP. An additional feature observed at t = 240, 300, and 440 fs is that the reflected flux is localized near θ ≈ 0, and moves preferentially along the MEP to exit the interaction region.
Figure 6 shows the expectation value of the coordinate ⟨R⟩ for the two surfaces, NR_PES and R_PES, as a function of time. As one can see, the minimum value of ⟨R⟩ ∼ 4.5 bohrs is reached at a time around 150 fs. Before that time the (initially identical) wave packets move along the two surfaces in the same way. During the period of time between 150 and 200 fs the packets are being reflected, such that ⟨R⟩ starts to increase, and a time delay between the two wave packets is building up: The wave packet traveling along NR_RES is delayed compared to the R_PES wave packet; i.e., it spends more time in the region of small values of ⟨R⟩. The relative time delay is about 20-30 fs, which is about 10% of the total time that the wave packets spend in the region of small ⟨R⟩, where reaction 2 of (upper panel) Reaction probability for reaction 2 with the total angular momentum J = 0, using the NR_PES (blue solid line), R_PES (green dashed line). and results from time-independent calculation from Guillon et al. 13 (red circles). (lower panel) Ratio of the reaction probability computed with the NR_PES to that one obtained with the R_PES. nuclei exchange takes place. The longer time spent in the interaction region means a larger reaction probability, which is consistent with the results shown in Figure 2.
From the time dynamics of the density functions, one can see that, although the topological difference between the NR_PES and the R_PES is rather subtle, features lying relatively deep in the well could have a significant consequence on the reaction dynamics, thus explaining the result in Figure 2. This gives a hint for understanding why the true or submerged barriers near the D e asymptote can change the temperature dependence for the thermally averaged rate coefficients. 13 Besides, although there is no barrier along the MEP for the reaction on the two surfaces, reflected fluxes are still observed due to the multidimensional nature of the reaction. Therefore, one can see that an one-dimensional study along the MEP for the reaction does not provide an appropriate picture for the reaction.
Our results also link to the collision time of the reaction, which is essential in some of the statistical approaches for reactive scattering. One 89,90 of such approaches implies that the lifetime of the collision complex is long enough such that the collision complex lives sufficiently long to redistribute the energy of the system randomly between different reaction channels. However, from our results, the time it takes for the system to arrive at the C 2v well is about 100 fs, while at about 400 fs the complex begins to dissociate. Therefore, the ratio between the average lifetime of the complex to the collision time is ≈ 300/100 = 3, which is not large enough for the assumption of randomized exit channels. For example, the wave packet dynamics shown in Figures 4 and5 clearly demonstrates that there is a preferential exit channel, along θ = 0, i.e., along the collinear geometry of 34 O 2 -16 O, which means that the scattering is very anisotropic and cannot be accounted for by the statistical approaches based on randomized exit channels. An anisotropy in the nuclei-exchange reaction was also observed in the experiment. 8,57,91 V. CONCLUSION In this work, we studied the wave packet dynamics for reaction 2 with 32 O 2 (ν = 0, j = 1) and total angular momentum J = 0 on the TKTHS PES 26 with and without the reef structure using the MCTDH method. For the NR_PES, a good agreement in the reaction probability at collision energy above 0.034 eV was found between our time-dependent result and the numerically well converged time-independent result obtained by Guillon et al., 13 who have also used this potential energy surface. But the time-dependent dynamics permits a supplementary insight into the problem.
One of the conclusions of the wave packet propagation concerns the time scale of the process, which in our study is 0.3-0.5 ps for both PESs. This is by an order of magnitude larger than a typical period of O-O bond vibration, confirming the results of the previous works that the exchange reaction (2) is not a direct one-dimensional kick-off process. The resonance scattering states of the intermediate metastable ozone O 3 * corresponding to times 0.24-0.3 ps in Figures 4 and5 thus play a crucial role for the process. Van Wyngarden et al. 57,67 have observed an anisotropy of angular distribution ("forward bias") of the reaction product in their experiment of crossing O and O 2 beams. This feature has also been found in theoretical investigations of the same reaction. 88 According to the interpretation of these experimental and theoretical results, the time scale of the exchange reaction was not sufficient for a full randomization of vibrational energy distribution within the O 3 * complex, which was a necessary condition of a simplified statistical approach. 51,89,90 On the theoretical side, Sun et al. 91 concluded that their state-specific quasi-classical calculations for the trajectories, which have lifetime t < 2 ps, would explain this "forward-bias" observation 67 (though they used the symmetry-forbidden j = 0 state of O 2 for simplicity). Our results provide a quantum mechanical counterpart consistent with these classical trajectory conclusions, as the wave packet reaction time of 0.4 ps in Figure 4 is indeed shorter than the t < 2 ps criterion of Sun et al. 91 and could thus contribute to the interpretation of experiments. 57,67 Our wave packet scattering lifetimes are also consistent with the lifetime distribution of a significant number of vibrational states of metastable ozone 48 O 3 * (main isotopologue) computed by Lapierre et al. 33 using the NR_PES. However, a full map of the metastable vibrational states of the 50 O 3 * complex enriched by 18 O isotope and their lifetimes were not yet calculated.
Another conclusion concerns the impact of the ozone PES shape on the reaction and the scattering resonances of the O 3 * intermediate complex. Both TKTHS PESs 26 are very similar at the bottom of the main C 2v well and possess the same dissociation threshold D 0 . The R_PES has a small submerged barrier much deeper in the well (Table 1) and is more attractive than other PESs possessing the reef structure. 5,21,25,66 The corresponding effect on the dynamics has not been yet used for modeling the exchange reaction (2) in previous studies.
In their quasi-classical study, Janssen et al. 17 have considered the exchange reaction cross section σ as a product of two factors, σ = σ cap P reac , where σ cap is the capture cross section, i.e., the cross section for trajectories entering the ozone well. Figure 4 shows that, for our quantum wave packet propagation study, σ cap is bigger for NR_PES than for R_PES. This is a nontrivial result as the R_PES is more attractive than the NR_PES (Figure 1) at large range distances near the D e asymptote. A quite small topographic reef feature with the amplitude of about 50 cm -1 sitting relatively deep in the well could thus have a significant impact on the reaction probabilities. In general, that the reaction probability is larger for the NR_PES (lower panel of Figure 2) explains the good agreement between experiments and the total reaction rates obtained by Guillon et al. 13 with this PES. At this point we see a correlation between dynamical results and the spectroscopy, as the NR_PES gave clearly more accurate prediction for the observed bands in the energy range approaching the TS. 31 We found that the topological differences between the NR_PES and the R_PES manifested mainly in the reflection of the wave packet near the reef region. Surprisingly, despite the contrast in the reaction probability between the two PESs, the visible differences in Figures 5 and6 are subtle as there is only slightly more flux reaching the C 2v well for the NR_PES than for the R_PES. Therefore, we concluded that resonance structures and the reaction probability of reaction 2 depend sensitively on the topological structure of the PESs, which is consistent with the conclusion by Guillon et al. 13 Consequently, since the formation of ozone at low pressure proceeds mainly through the Lindemann mechanism, which involves the metastable O 3 *, we expect that the three-body recombination rate coefficient at low pressure may also critically depend on the shape of the ozone PES. It is therefore important to use an accurate PES to calculate resonance energies and widths of metastable O 3 * with different isotopic substitutions. To this end, new spectroscopic measurements and assignments of 16 O 16 O 18 O and 16 O 18 O 16 O ozone isotopomers in the energy range near the dissociation threshold would be extremely helpful to validate ab initio PESs.
Ultimately, using reliable calculations of scattering resonance energies and widths of the metastable 48 O 3 * complex enriched by 18 O isotope (not yet available on spectroscopically accurate PESs), one should be able to explain the MIF effect, 49,50,62 which is the major challenge for the ozone formation dynamics. 9,[46][47][48]51,54,55 In this context an investigation of long-lived resonances will be important, for example, those linked to the "roaming" nuclear motion 92,93 considered in the classical orbit studies by Mauguiere et al. 94 Although the present results suggest that the assumption about the random distribution over exit channels, used in previous statistical approaches 51,55 applied to O 3 , is not appropriate for the 32 O 2 + 18 O → 34 O 2 + 16 O reaction, the other types of statistical approaches, taking into account the branching ratios over the final channels, may be appropriate and can simplify significantly the dynamics of O + O 2 collisions and the ozone formation in three-body collisions. Theoretical studies, taking advantages of quantum scattering and some reasonable statistical assumptions, describing the ozone dynamics and the dynamics of three-body collisions, such as O 2 + O + N 2 → O 3 + N 2 , are highly desirable.
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.9b06139.
Movies of the density function of the system in (r, R) and (θ, R), propagated along the NR_PES of the TKTHS PES (ZIP)
Corresponding Author *E-mail: vladimir.tyuterev@univ-reims.fr.
Chi Hong Yuen: 0000-0002-0544-4976
The authors declare no competing financial interest.
2 + → * → +(1)If a heavy oxygen atom18 O hits the symmetric diatomic molecule32 O 2 composed of two16 O atoms, one16 O could be kicked out and replaced by isotopic substitution: