Sugata
Goswami
a,
Juan Carlos
San Vicente Veliz
a,
Meenu
Upadhyay
a,
Raymond J.
Bemish
b and
Markus
Meuwly
*ac
aDepartment of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel, Switzerland. E-mail: m.meuwly@unibas.ch
bAir Force Research Laboratory, Space Vehicles Directorate, Kirtland AFB, New Mexico 87117, USA
cDepartment of Chemistry, Brown University, RI, USA
First published on 14th September 2022
The dynamics of the C(3P) + O2(3Σ−g) → CO(1Σ+) + O(1D) reaction on its electronic ground state is investigated by using time-dependent wave packet propagation (TDWP) and quasi-classical trajectory (QCT) simulations. For the moderate collision energies considered (Ec = 0.001 to 0.4 eV, corresponding to a range from 10 K to 4600 K) the total reaction probabilities from the two different treatments of the nuclear dynamics agree very favourably. The undulations present in P(E) from the quantum mechanical treatment can be related to stabilization of the intermediate CO2 complex with lifetimes on the 0.05 ps time scale. This is also confirmed from direct analysis of the TDWP simulations and QCT trajectories. Product diatom vibrational and rotational level resolved state-to-state reaction probabilities from TDWP and QCT simulations agree well except for the highest product vibrational states (v′ ≥ 15) and for the lowest product rotational states (j′ ≤ 10). Opening of the product vibrational level CO(v′ = 17) requires ∼0.2 eV from QCT and TDWP simulations with O2(j = 0) and decreases to 0.04 eV if all initial rotational states are included in the QCT analysis, compared with Ec > 0.04 eV obtained from experiments. It is thus concluded that QCT simulations are suitable for investigating and realistically describe the C(3P) + O2(3Σ−g) → CO(1Σ+) + O(1D) reaction down to low collision energies when compared with results from a quantum mechanical treatment using TDWPs.
One of the latest investigations for the entire C(3P) + O2(3Σ−g) ↔ CO2 ↔ CO(1Σ+) + O(1D)/O(3P) reaction was carried out using quasi-classical trajectory (QCT) simulations with kernel-represented PESs based on extensive multi reference configuration interaction (MRCI) calculations.21 The focus was more on the high-temperature regime of the reaction also because experimental and computational data from shock tube experiments for the C + O2 → CO + O reaction between 1500–4200 K was available.9,13 The latest QCT simulations reported the forward and backward thermal rates for the lowest five electronic states together with vibrational relaxation times.21 However, the state-to-state dynamics was not considered and the question remains whether a classical framework is suitable when compared with results from a quantum mechanical treatment of the nuclear dynamics.
One recurrent theme for atom + diatom reactions concerns the range of applicability of quasi-classical-based dynamics approaches for computing thermal rates and final state distributions. Cross sections and thermal rates are averaged over initial and/or final ro-vibrational resolved state-to-state information. Thus, although more or less heavily averaged quantities may favourably agree between different approaches, it is possible that such agreement arises from the averaging process. Hence, it is also relevant to compare properties at the state-to-state level. All such quantities are essential as input to more coarse-grained investigations of reaction networks as they appear in hypersonics, combustion, atmospheric and astrophysical chemistry.30 Because QCT simulations are computationally more efficient, they are often used instead of and also along with more time-consuming quantum nuclear dynamics simulations.31–35
Quantum mechanical and QCT approaches were used to evaluate reaction probabilities and rates for the C + OH reaction and were found to be in fairly good agreement.36–38 For the S + OH reaction quantum mechanical and QCT simulations produced excellent agreement for cross sections but only fair agreement for the total reaction probabilities and thermal rates.32,39 Comprehensive QCT and QM dynamics investigations for the Br + H2 and O + HCl reactions reported favourable agreement between QM and QCT results for properties such as reaction probabilities, integral and differential cross sections.40,41 Similar to the S + OH reaction, for N + H2 results from QM and QCT simulations agree favourably for cross sections but not so well for reaction probabilities.42
For the S + OH reaction it has been specifically reported that for reactive collisions a more general relationship between the mechanistic details of the dynamics and the ensuing rates can be difficult to obtain.32 Nonetheless, it is generally believed that the dynamical features of an exoergic and barrierless reaction primarily depend on the masses of the participating atoms, the exoergicity and the topographical details of the underlying PES.32,43,44 To better understand the origins of the observed dynamics systematic investigations have been carried out for different systems.38,43–47 As an example, for the C + OH reactive collisions38,43,44 it was found that reagent vibrational excitation decreased reactivity on the first electronically excited PES and enhances reactivity on the second excited state. The excess vibrational energy is transferred into product translational energy for the first excited PES whereas for the second excited state it is transferred into product vibration and rotation. These dynamical effects on the final states are caused by the topology of the underlying PESs which have no barrier to products for the first electronically excited state but involve a late barrier for the second excited state.
The reaction C(3P) + O2(3Σ−g) → CO(1Σ+) + O(1D)/O(3P) using its five lowest lying states (1A′, 3A′, 3A′′, 1A′′ and (2)1A′) is a particularly suitable system for such an investigation due to the availability of high-quality, full-dimensional PESs which were recently validated vis-a-vis experiment.21 The exoergicity of the three singlet and two triplet states are ∼3.8 eV and ∼5.9 eV, respectively.21 Therefore, the difference in the dynamical attributes on the singlet and triplet states is expected to arise solely from the topographical details of the underlying surfaces.32
Quantum mechanics-based dynamics approaches can become computationally expensive for barrierless exoergic reactions with deep wells.32,38,43,44,48,49 Thus, QCT approaches are an attractive alternative, at least at the qualitative level. In the present work the initial state-selected and state-to-state reaction probabilities for the C(3P) + O2(3Σ−g) → CO(1Σ+) + O(1D) reaction on the electronic ground (1A′) state are determined for low to moderate collision energies (Ec ≤ 0.4 eV) by TDWP and QCT approaches. Total reaction probabilities and final state distributions are determined for a range of collision energies and allow to compare the two dynamics methods in terms of their findings and with experiment. The work is organized as follows. First the methods used are described. Then, the convergence of the TDWP simulations is assessed, followed by results for the total reaction probabilities and the product vibrational resolved reaction probabilities from TDWP and QCT simulations. Next, product vibrational and rotational state distributions are determined and compared. Finally, the computational results are discussed in the context of experiments and in a broader sense.
![]() | (1) |
![]() | ||
Fig. 1 Reactant (left) and product (right) collision geometries in body-fixed Jacobi coordinates for the present TDWP simulations. |
The initial WP (eqn (1)) is subsequently transformed to product state Jacobi coordinates (RC, rC, γC) with C as atom and AB as diatom, see right-bottom part of Fig. 1, according to52
![]() | (2) |
The action of the nuclear Hamiltonian Ĥnuc, formulated in product Jacobi coordinates, Ĥnuc, on the WP is52,55
![]() | (3) |
The propagated WP at the first time step, τ1, is evaluated according to
![]() | (4) |
![]() | (5) |
As is customary for TDWP simulations, spurious reflections of the WP components from the finite sized grid edges at longer propagation time need to be damped. This is done by multiplying the wavepacket with a double exponential damping function52a(y) along coordinates y = Rc, rc
a(y) = exp[−cabsexp(−2(ymax − yabs)/(y − yabs))] | (6) |
The propagated WP is projected on ro-vibrational levels of the product diatom at the product asymptote (R = R∞) to obtain the time-dependent coefficients52
![]() | (7) |
![]() | (8) |
Hamilton's equations of motion were solved using a fourth order Runge–Kutta method. The initial conditions for initiating the trajectories were sampled by using standard Monte Carlo methods63 and the ro-vibrational levels of reactant and product diatoms are calculated using semiclassical quantization. The simulations are run at several collision energies with a time step of Δt = 0.05 fs which guarantees the conservation of total energy and angular momentum. At each collision energy, 5 × 105 trajectories are run for converged results, except at the two lowest energies, 0.0055 and 0.001022 eV for which 106 trajectories are run for convergence of the total reaction probability for O2(v = 0, j = 0). It is found that the difference between the total reaction probabilities at these two energies obtained from 5 × 105 and 106 trajectories appears at the fourth decimal place and therefore all remaining calculations are carried out with 5 × 105 trajectory simulations. As the associated quantum numbers are real-valued, their necessary assignment to integer values was made either by using Gaussian binning (GB)61,64,65 which centers Gaussian weights around the integer values and has a full width at half maximum of 0.1 or from Histogram binning (HB) which rounds to the nearest integer.
The reaction probability at collision energy Ecol is obtained as
![]() | (9) |
Parameter | Value | Description |
---|---|---|
N R/Nr/Nγ | 499/419/240 | Number of grid points for product Jacobi coordinates (RC, rC, γC) |
R min/Rmax (a0) | 2.2/22.0 | Extension of the grid along RC |
r min/rmax (a0) | 1.5/17.9902 | Extension of the grid along rC |
R ∞ (a0) | 12.5 | Location of the dividing surface in the product channel |
R abs/rabs (a0) | 14.0/14.75 | Starting point of the absorption function along RC and rC |
C abs/cabs | 30.0/30.0 | Strength of the absorption along RC and rC |
R 0 (a0) | 13.2 | Center of the initial wave packet in reactant Jacobi coordinate |
E trans (eV) | 0.115 | Initial translational energy in eV |
α | 12.0 | Width of the initial wave packet |
β s | 0.5 | Smoothing of the initial wave packet |
Convergence runs for all parameters were carried out for J = 0 and for reactant O2(v = 0, 1, 2, j = 0). The converged parameters are summarized in Table 1 and demonstrate that extensive grids along the three product Jacobi coordinates are required resulting in computationally expensive calculations. The requirement for such large grids is the result of the heavy masses of the interacting atoms and the presence of the deep potential well.21 The large number of grid points required along the angular Jacobi coordinate reflects the pronounced angular anisotropy of the underlying PES.
Two-dimensional (2-D) contours of the underlying PES are shown in Fig. 2. The top row reports the PES in reactant Jacobi coordinates for ROO = rA = 3.40 a0, i.e. a stretched geometry of the reactant diatom before C-insertion occurs to form O–C–O (minimum at 90°). Because after formation of the triatom the system only spends little time in this state (typically <0.5 ps, see Discussion section) to move towards the product, it is also of interest to provide a representation of the reactant PES in terms of product Jacobi coordinates (RC,γC) which is shown in the top right panel in Fig. 2. Comparison of the two representations clarifies that spatial symmetry is lost for the product state coordinates and the anisotropy of the PES differs considerably. Similarly, the bottom row reports the PES in product state coordinates for RCO = rC = 3.0 a0, i.e. for the system on its way to form products (left panel) and the same PES represented in reactant state coordinates (RA,γA) in the right hand lower panel of Fig. 2. In addition, representative structures for CO2 are reported and their location on the PES are indicated by arrows. Two geometries (see black and purple crosses in the top row panels) are chosen to show their positions in the reactant and equivalent product Jacobi coordinates.
In the production runs, for O2(v = 0, j = 0) the number of product vibrational levels included v′ = 0–18 which was suitably increased to v′ = 19 and v′ = 20 with the reactant in its two vibrationally excited states. The WP was propagated for up to 70000 time steps (corresponding to ∼1.4 ps) in order to obtain converged reaction probabilities for J = 0 and the potential and/or centrifugal cut-off was 0.57Eh.
Fig. 3 shows that starting at low collision energy (Ecol = 0.001 eV ≈ 10 K) the TDWP probability for O2(v = 0, j = 0) rises sharply and oscillates around 0.95 up to Ecol ∼ 0.1 eV. Afterwards, the total reaction probability decreases slowly with undulations and reaches P ∼ 0.7 at Ecol ∼ 0.4 eV. The non-zero probability for Ecol ∼ 10−3 eV (Fig. 3) reflects the barrierless nature of this exoergic reaction.21 The barrierless and exoergic features of this reaction along with the deep potential well (formation of CO2) on the underlying PES (cf., Fig. 1 of ref. 21) and the large masses of the participating atoms make the quantum dynamical calculation computationally challenging. The undulations in the total reaction probability (see Fig. 3) are indicative of formation of intermediate collision complexes inside the deep potential well.66 However, the higher exoergicity of the reaction facilitates breakup of these complexes and thus shortens their lifetimes. Despite a deeper well on the underlying PES, lower exoergicity could change these undulations into sharp and intense resonance oscillations which was found for the S + OH, C + OH, C + H2, and S + H2/D2/HD reactions.32,43,44,67,68
The reaction probabilities from the wavepacket simulations in Fig. 3 reveal that reactant vibrational excitation decreases the reactivity at low and intermediate collision energies, whereas the reactivity is enhanced for O2(v = 2, j = 0) at the highest collision energies considered. Undulations for the total reaction probability are also observed for vibrationally excited reactants and analyzed further below.
Next, the reaction probabilities from the QCT simulations are considered and compared with the TDWP results. QCT reaction probabilities are largest for O2(v = 0, j = 0) and decay with increasing collision energy, see green (GB) and red (HB) circles in the inset of Fig. 3. For O2(v = 1, j = 0) (blue circles) at low Ecol the reaction probability is smaller than for O2(v = 0, j = 0) but approaches it with increasing collision energy. Interestingly, for O2(v = 2, j = 0) (brown circles) the QCT reaction probability starts low at low Ecol, increases – similar to TDWP – and then also decays towards a comparable amplitude as for v = 1 and v = 0. Compared with the results from TDWP simulations the QCT simulations are in excellent overall agreement except for the lowest collision energies considered. In contrast to TDWP, the QCT results decay monotonically for higher collision energies also for O2(v = 0, j = 0), (see inset of Fig. 3). With increasing collision energy the QCT probability for all initial v-states converges to the same value (reaction probability ∼0.6). This is also apparently found for the TDWP simulations except for O2(v = 2). This possibly happens as the opening of energetically higher product (CO) internal levels for vibrationally excited reagent becomes easier at higher collision energies which enhances the reaction probability as a result. Increase of reaction probability (with oscillations) with collision energy is also seen in TDWP investigation of C + OH reaction on its first and second excited states.44 The C + OH reaction on its first (barrierless) and second (small barrier at the exit channel) excited states is also exoergic and proceeds through potential wells on the underlying potential energy surface.
The effect of initial translational/collision energy and vibrational energy on reactivity depending on the shape of the underlying PES has been discussed from the perspective of Polanyi's rules: Translational energy is more effective to promote the reactivity for an early barrier reaction, whereas the effect of reactant vibrational energy is more relevant for a reaction with a late barrier.45 Furthermore, a more recent investigation of a number of atom + diatom reactions46 reveals the importance of a cross over point for the effects of reagent vibrational and translational/collision energy on the reactivity. Therefore, the interplay between initial translational and vibrational energy can play a significant role in reaction dynamics. The present TDWP simulations suggest that with increasing reactant vibrational excitation the reaction probability slightly increases above Ecol ∼ 0.3 eV (which approximately corresponds to the energy of the O2(v = 1) state), also due to the more delocalized nature of the wavepacket which allows for some coupling between translation and vibration. For the QCT simulations such an effect is not observed which is probably due to the short contact time (∼50 fs, see below) which essentially precludes transfer between vibrational and translational degrees of freedom.
A comparison of the state-to-state TDWP and QCT reaction probabilities in Fig. 4 shows that the results are in excellent agreement with each other at all collision energies considered for v′ = 0–8, which is remarkable. However, differences at low collision energies (Ecol < 0.05 eV) appear for v′ = [9, 10, 12–14] in that the TDWP results find reaction probabilities that decay to close to zero for the lowest collision energies followed by a pronounced peak whereas for the QCT simulations the reaction probability does not decay to zero for the lowest collision energy considered. A second observation concerns the reaction probability for v′ = 16 for which the quantum results feature a maximum in the reaction probability between Ecol = 0.05 and 0.1 eV which is not present in the results from QCT simulations. However, this is expected as QCT simulations are unable to capture features related to resonances due to the neglect of coherences.
The two panels of Fig. 4 for CO(v′ = 16, 17) reveal that TDWP find energy thresholds for product formation for Ecol ∼ 0.014 and ∼0.22 eV, respectively. This is in quite good agreement with a value of 0.19 eV for v′ = 17 from recent QCT simulations.21 With histogram binning for analyzing the QCT results there is no threshold for CO(v′ = 16) but generating CO(v′ = 17) requires collision energy. Analyzing the same data with GB finds thresholds for CO(v′ = 16, 17), see Fig. 4. For the reaction with vibrational energy in the reactant state C + O2(v = 1, j = 0), see Fig. S2 (ESI†) the product diatom vibrational level-resolved reaction probabilities from TDWP and QCT simulations compare as well as for O2(v = 0, j = 0). It can also be seen that reactant vibrational excitation has an insignificant effect on product vibrational level resolved reaction probabilities.
The feature that agrees most closely between the TDWP and QCT approaches is the value for v′ beyond which the final state probability P(v′) decays to zero. Also, for low collision energy and most of the distributions emerging from O2(v = 0, j = 0) the TDWP and QCT results for low v′ agree favourably. On the other hand there are also differences in P(v′) from the two approaches which concern P(v′) for low vibrational quantum numbers, see e.g. all initial v-states for Ecol = 0.4 eV or most of the distributions for O2(v = 2, j = 0). In addition, there are features such as the increase in P(v′) for O2(v = 0, j = 0) at Ecol = 0.3 eV around v′ = 10 which are present in the TDWP simulations but absent in the QCT simulations.
It can also be seen from Fig. 5 that with increasing vibrational excitation of the reactant diatom the number of product vibrational state increases. The energy supplied to reactant vibration enhances the total energy of the system and makes the opening of vibrationally excited product channels possible. This suggests that part of the reactant internal energy is transferred/disposed into product vibration. However, a careful look at Fig. 5 reveals that the collision energy and reactant vibrational excitation hardly have any effect on the overall pattern of the distribution and thus it becomes very difficult to find a general/regular trend for the variation of the distribution with these two types of energy.
It can also be seen from each row of Fig. 6 that reactant vibrational excitation does not change the pattern of the distributions. This is observed in both TDWP and QCT results. Thus, similar to the vibrational distributions (Fig. 5), the overall pattern of the rotational distributions is also unaffected by reactant vibrational excitation. However, the peak height of the TDWP distributions at the low j′ gets lowered and the distribution becomes flatter at Ecol ∼ 0.4 eV. This is more prominent when the reactant is in its ground ro-vibrational level. Therefore, at higher collision energies the diatom is formed with lower and moderate j′ with almost equal possibility.
In addition to the forward C + O2 reaction, which was the topic up to this point, the thermoneutral atom-exchange, COA + OB → COB + OA, and the endoergic backward COA + OB → C + OAOB reactions also take place on the 1A′ PES. To complete the picture, total reaction probabilities from QCT simulations for these two reactions were calculated up to 6.5 eV collision energy and the results are shown in Fig. 7 as function of collision energy. It can be seen from Fig. 7 that the COA (v = 0, j = 0) + OB → C + OAOB reaction has a high energy threshold whereas the atom-exchange reaction is barrierless. Comparison of probabilities shown in Fig. 3 and 7 reveals that the reactivity is highest for the exoergic C + O2→ CO + O reaction, followed by the atom-exchange reaction. The TDWP approach can also be applied here, however, the parameters shown in Table 1 need to be converged separately for the atom-exchange and endoergic backward reactions which is outside the scope of the present work. Given the good agreement between TDWP and QCT simulations it is, however, expected that the results from Fig. 7 are representative.
The thermoneutral atom-exchange reaction (COA + OB → COB + OA) has already been considered in the context of stabilizing CO2 on amorphous solid water (ASW) to better understand molecule formation at astrophysical conditions.70,71 For this process it was found that for bulk temperatures of ∼50 K formation and stabilization of CO2 is very likely (60% of the trajectories stabilize). Together with the present results this suggests that quantum effects will only marginally affect this conclusion. This is also consistent with results from QCT and TDWP simulations for O(3P) + O(3P) recombination for which quantum mechanical effects were also found to be far less important than effects of ASW surface roughness on the diffusivity and recombination dynamics.72,73
The undulations in the total reaction probabilities from TDWP simulations (Fig. 3, 4 and Fig. S2, ESI†) are attributed to formation and transient stabilization of the CO2 intermediate inside the potential well. These are largely unaffected by reactant vibrational excitation, and merit some further discussion. Previous work on the N + OH reaction66 has associated a dense pattern of resonances between 0.41 and 0.43 eV with formation of long-lived intermediates in the deep HON well together with a pronounced exoergicity of 1.99 eV of the reaction without, however, providing further support for this interpretation. In the present case the undulations in the total reaction probability are considerably more pronounced and regular and the reaction has an exoergicity of ∼3.8 eV which warrants further exploration. The black-solid curve in Fig. 3 features 19 prominent peaks across the entire range of collision energies considered. Peak separations range from ∼0.009 eV to ∼0.036 eV with an average of ∼0.021 eV. Using the energy/time uncertainty relationship ΔE × Δt = ħ, corresponding lifetimes range from ∼0.018 ps to ∼0.073 ps with an average of ∼0.038 ps.
The three internuclear separations from representative QCT simulations for Ecol = 0.05 eV and Ecol = 0.40 eV are shown in Fig. 8 and the corresponding collision time distributions P(τc) are reported in Fig. S3 (ESI†). Here, τc is defined as the time for which the sum of the three internuclear distances is smaller than 14a0 or 12a0, respectively. The distribution P(τc) was determined from 22000 reactive trajectories at each of the collision energies. The time series in Fig. 8a and d represent trajectories with short collision time τc ∼ 0.03 ps and ∼0.04 ps, respectively, whereas Fig. 8b and e show trajectories with longer τc (∼1.0 ps and ∼1.02 ps, respectively). The results of panels b and e clearly indicate the formation of triatomic intermediate complexes whereas the findings of panels c and f suggest the formation of weakly bound triatomic complexes with short contact times, reminiscent of roaming.
Formation of the collision complex is found from both, TDWP (panel A) and QCT simulations (panels C and D). For this, the deviation ΔA of the C–OA separation from the equilibrium bond length in CO2(OA–C–OB), which is21 2.20 a0, was considered, and similarly for C–OB. From this, the difference ΔA + ΔB was determined and is shown in Fig. 9A for TDWP and in panels C and D for QCT simulations. For small ΔA + ΔB or values close to zero the structure is that of the CO2 intermediate, whereas for large values of ΔA + ΔB the system is in its two asymptotic states COA + OB or COB + OA. The asymptotic region corresponds to ΔA + ΔB ∼ 15a0 whereas ΔA + ΔB ∼ 5 a0 is indicative of CO2 formation.
For the analysis of the quantum simulations the position of the maximum of the wavepacket was extracted. Every 100 time step of the propagation the position (R(m), r(m), γ(m)) of the maximum of the wavepacket in product state Jacobi coordinates was determined from which separations C–OA and C–OB are obtained. For the QCT simulations a total of 100 reactive trajectories with final states COA + OB (50 trajectories) and COB + OA (50 trajectories) were run. Because the reaction time for each trajectory differs, all time series for ΔA + ΔB were first synchronized in time by shifting the time to zero when they reach ΔA + ΔB = 15.0 a0. Then, histograms as a function of time were generated, see Fig. S4 (ESI†), and their centers of gravity of the distributions are reported in Fig. 9C and D.
Fig. 9A shows that the sum of deviations from the CO2 equilibrium structure flattens out between ∼200 fs and 300 fs after launching the TDWP. During this time the van der Waals separation R (lower panel in Fig. 9B) first decreases to 3.5a0 and stabilizes there for the next 75 fs. This is accompanied by a change in the angle to γ ∼ 15° indicative of linear CO2. This linearity is established between 190 fs and 265 fs after launching the wavepacket. From all this data an approximate lifetime of ∼75 fs is assigned to the intermediate from TDWP.
From QCT simulations (Fig. 9C and D) it is first found that with increasing collision energy the collision complex forms more rapidly. Secondly, a collision complex as seen in Fig. 8, corresponds to values of ΔA + ΔB ≲ 5 a0 or shorter. This further confirms that the signatures for R ∼ 3.5 a0 in the TDWP simulations indeed correspond to the collision complex. Due to the large energy (11.22 eV) released upon formation of CO2 from C + O2 the system is formed in a highly non-equilibrium state which is unlikely to sample the CO2 minimum energy structure. Hence, sampling values ΔA + ΔB ∼ 0 is very unlikely. This is found for both, TDWP and QCT simulations. Finally, when comparing results from TDWP and QCT simulations it should be remembered that QCT simulations are always run for one specific value of the collision energy whereas in a TDWP all collision energies considered are propagated together. Hence, only the superposition of results from QCT simulations for a range of collision energies can be expected to give information comparable to that from a TDWP.
The collision time distributions from QCT simulations are strongly peaked with their maxima at τc ∼ 0.07 ps and ∼0.056 ps for Ecol = 0.05 eV and 0.40 eV, respectively, see inset of Fig. S3 (ESI†). The probability for finding such short collision times is 20% and 30%, respectively, for the two collision energies considered. If a shorter distance cutoff of 12a0 is used the maxima of P(τc) are at τc = 0.056 ps whereas the longest lifetimes are 4.5 ps and 3.2 ps for the two collision energies. Hence, the analysis of the lifetimes is only moderately affected by the choice of the cutoff value for the sum of the three interatomic distances.
Collision times of the order of 0.06 ps are in quite good qualitative agreement with the analysis of the recurrences in the reaction probabilities from TDWP simulations which yielded lifetimes of around ∼0.04 ps, see above. On the other hand, the QCT simulations also find a long tail in the collision time distributions extending out to ∼4.8 ps and ∼5.3 ps for the longest collision times for the two simulation conditions considered. Overall, it is found that the implied time scales of the undulations in the TDWP reaction probabilities, which have previously been linked to formation of a collisionally stabilized intermediate,66 are consistent with the collision times from the QCT simulations, and that they constitute an inherent characteristic of the reaction considered.
The present results can also be compared with results from pulsed, crossed, supersonic molecular beams, at relative translational energies between 4.4 and 90 meV.74 Such experiments found that collision energy Ec > 0.04 eV is required to open the product CO(v′ = 17) channel. For CO(v′ = 17), the present QCT investigation shows a threshold of ∼0.15 eV from analysis based on histogram and ∼0.25 eV from Gaussian binning, respectively. An earlier QCT investigation21 on a finer grid for Ec showed that for (v = 0, j < 10) the threshold is 0.18 eV which changes to 0.25 eV if the same data is analyzed for (v = 0, j = 0). This is also consistent with ∼0.22 eV from the TDWP simulations. Including higher rotational states for the reactant progressively shift the threshold to 0.04 eV when all reactant j-states are included which agrees well with the experiments.21,74
In summary, this work provides a stringent comparison for total and product diatom vibrational state resolved reaction probabilities and product state vibrational and rotational distributions from TDWP and QCT simulations for the C + O2 → O + CO reaction on its ground electronic state. Overall, the agreement of a quasi-classical description with results from wavepacket simulations is encouraging except for the lowest collision energies or lowest product rotational states. Expected quantum effects such as undulations in the total reaction probability P(Ec) can be rationalized from analysis of the TDWP and the QCT trajectories. Given the ease with which QCT simulations can be carried out the present work suggests that they are a meaningful approach for understanding the dynamics of atom plus diatom reactions and for larger systems as well. QCT-based approaches have also been used for the reaction of H2 on platinum surfaces,75 the formation of CO2 on amorphous solid water,70,71 CO adsorbed on NaCl,76 and for investigating photodissociation reactions of larger molecules such as acetaldehyde or the syn-Criegee intermediate.77,78
Footnote |
† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp02840a |
This journal is © the Owner Societies 2022 |