The puzzle of high lifetime and low stabilization of HO₃˙: rationalization and prediction

Abstract

One of the most important puzzles in atmospheric chemistry is the long-lifetime of HO₃˙ in spite of its low-stabilization energy. In the present work, we have estimated the lifetime of HO₃˙ using classical dynamics simulations by coupling an available neural-network analytical potential energy surface with a chemical dynamics program. The simulation results clearly indicate that at room temperature, the lifetime of HO₃˙ can exceed 1 μs under collision-free conditions. In fact, at 200 K, the lifetime of HO₃˙ can enter the millisecond timescale. This suggests that HO₃˙ is indeed a stable enough intermediate that can affect the outcomes of crucial atmospheric processes.

---

1. Introduction

HO₃˙ is a weakly bound complex between an OH radical and O₂ molecule, and hence, can act as a reservoir for OH radicals. Besides, HO₃˙ is also postulated to be a key intermediate in crucial upper atmospheric reactions.^1–11 For example, the HO₂˙ + O → OH˙ + O₂ reaction has been suggested to occur through the formation of a HO₃˙ transient intermediate rather than via a direct H-atom transfer.³ The ability of HO₃˙ to act as a reservoir, as well as its potential to affect the outcome of atmospheric processes, crucially depends on its lifetime, which in turn depends on its stabilization energy. For example, Meinel band emission¹² occurs due to the formation of vibrationally hot OH radicals by H + O₃ → OH˙ (ν) + O₂ reaction. Interestingly, although this reaction is expected to produce OH˙ in high vibrational states (ν = 8, 9), the actual OH˙ populations are primarily found in lower vibrational levels (ν = 1, 2).¹³ There are various attempted explanations for this observation.^14–17 One of them is that it occurs due to the collisional relaxation of OH˙ (ν) by O₂ in the atmosphere via the formation of a HO₃˙ intermediate.¹⁷ The status of HO₃˙ from a postulated molecule to a true intermediate, was updated in 1999, when Cacace et al.² detected it in the laboratory via neutralization-reionization and neutralization-reionization/collisionally activated dissociation mass spectrometry. The experiment suggests that at 298 K, HO₃˙ can survive for more than one microsecond under collision-free conditions. The microsecond lifetime clearly qualifies HO₃˙ as a reasonably stable intermediate, and thus supports its effectiveness in influencing atmospheric processes. Surprisingly, this observed stability of HO₃˙ is not consistent with recent experimental⁹ and theoretical studies,^10,18,19 which estimate its binding energy to be ∼3.0 kcal mol⁻¹. Based on statistical thermodynamics and this binding energy, the predicted lifetime of HO₃˙ is only a few picosecond, far too short compared to the experimental value. This suggests that there must be other factors beyond stabilization energy that dictate the stability of HO₃˙. To understand the origin of this puzzle, in a recent work,¹¹ on-the-fly classical trajectory simulations were performed. The study found that HO₃˙ dissociation does not follow standard statistical behaviour as described by RRKM theory. Moreover, the non-statistical effect increases as one goes from the higher energy to lower energy region. For example, at 40 kcal mol⁻¹, the lifetime of HO₃˙ due to the non-RRKM effect was found to increase by a factor of ∼2, whereas at 15 kcal mol⁻¹, it increased by more than an order of magnitude. Therefore, one can expect that at atmospherically relevant temperatures (200–300 K), where the true importance of HO₃˙ lies, non-RRKM effects would be even more pronounced and the lifetime would be substantially longer. However, on-the-fly simulations were only able to reach down to 15 kcal mol⁻¹, while the average thermal energy corresponding to 200–300 K is ∼2.4–3.6 kcal mol⁻¹. The major challenge in not going below 15 kcal mol⁻¹ with on-the-fly dynamics was the unrealistic computational cost associated with simulating long timescales at lower energies. This issue can be circumvented using an analytical full potential energy surface (PES). Fortunately, various full analytical PESs are available in the literature.^20–23 Therefore, in the present work, to overcome previous limitations, we employed an analytical PES of HO₃˙ and performed classical trajectory simulation by coupling it with a chemical dynamics program. This approach allows us to perform trajectory simulations across a broad range of energies and timescales. In addition, it also enables us to estimate rate constants as a function of temperature, which is crucial for direct comparison with experimental lifetime. We believe that the present work will further help in bridging the gap between experiment and theory and will enhance the understanding of the fate and impact of HO₃˙ in the atmosphere.

2. Methodology

2.1. Potential energy surface

There are two main components in performing dynamics simulation; the selection of an appropriate potential energy surface (PES), and carrying out classical trajectory simulations by coupling the chosen PES with a chemical dynamics software. To the best of our knowledge, four PESs for HO₃˙ are currently available in the literature. The previous study¹¹ suggests that there are two crucial factors that are key in governing the dynamics of HO₃˙: cis–trans stability and cis–trans isomerization barrier. Therefore, for the suitability of the PES, we have compared these parameters for all four PESs. A comparative analysis of these two parameters for all four PESs is illustrated in Fig. 1. To assess the reliability of these PESs, we have used zero-point energy (ZPE) corrected values obtained at the CCSDTQ(P)/CBS level of theory as our benchmark (shown in black).¹⁰ In microcanonical sampling, by definition, zero-point energy cannot be explicitly included. The total energy is defined relative to the bottom of the harmonic oscillator potential for each vibrational mode.^24,25 Therefore, we need to choose ZPE-uncorrected PES values for comparison with the values of the ZPE-corrected benchmark PES and experimental results. The PES in purple (superscripted by d) was developed by Hua-Gen and Varandas²⁰ in 2001 at the QCISD(T)/CBS level of theory, while the PES in violet (superscripted by e) was constructed by Braams and Hua-Gen²¹ in 2008, employing the HCTH/aug-cc-pVTZ level of theory. It is evident from Fig. 1 that the PES of Braams and Hua-Gen clearly overestimates the stabilization energy of both cis and trans conformers. The Hua-Gen and Varandas PES is better at predicting the stabilization and dissociation energy of HO₃˙, but it suggests that the cis-conformer is more stable than the trans one, a fact that contrasts with benchmark results. On the other hand, the PESs in red and blue were recently developed using higher-level theory and neural-network fitting.^22,23 In fact, the blue PES (superscripted by c) was developed by Xixi Hu et al.²² in 2019 by fitting 2087 energy points using the permutation invariant polynomial-neural network (PIP-NN) approach. These energy points were obtained using the MRCI(Q)-F12/VDZ-F12 level of theory. The red PES (superscripted by b) is the most recent one and was constructed by Zuo et al.²³ in 2020, employing PIPs in the input layer of a neural-network (NN). This global multichannel PES was fitted using more than 21 000 ab initio points at the MRCI(Q)-F12/VDZ-F12 level of theory. Fig. 1 clearly suggests that these PIP-NN-PESs are more closely aligned with the benchmark results. In particular, the most recent PES (Zuo et al.) shows excellent agreement with both benchmark theoretical results and experimental measurements. For example, the ZPE-uncorrected binding energy of the Zuo et al. PES is ∼3.8 kcal mol⁻¹, which is reasonably closer to the experimental dissociation energy (D₀) of ∼3.0 kcal mol⁻¹. Similarly, the ZPE-uncorrected cis–trans and trans–cis isomerization barriers of the PES are ∼0.9 and 1.1 kcal mol⁻¹, respectively, which are in reasonable agreement with the ZPE-corrected values obtained from benchmark theory, i.e., 0.3 and 0.8 kcal mol⁻¹. Therefore, in the present work, we have used the PES by Zuo et al. to carry out the trajectory simulation by coupling it with the classical dynamics program VENUS.^26,27

Fig. 1 Comparison of the four potential energy surfaces for HO₃˙ (ZPE-uncorrected) developed by: (e) Braams and Hua-Gen Yu,²¹ (d) Hua-Gen Yu and Varandas,²⁰ (c) Xixi Hu et al.,²² and (b) J. Zuo et al.,²³ with benchmark calculations from (a) Bartlett et al.¹⁰ The benchmark results of Bartlett et al. are ZPE-corrected.

2.2. Classical dynamics simulation

In the present work, we have performed dynamics simulations using classical microcanonical sampling. In microcanonical normal mode sampling, it is assumed that the system consists of n separable harmonic oscillators. The total energy is distributed randomly among the harmonic oscillator modes. After the energy distribution, the energy of each harmonic oscillator is converted into the corresponding coordinates and momenta (see SI for details). To perform dynamics simulation, the starting structures were the cis-conformers of HO₃˙. Although we have started the simulation with the cis-conformer, we believe it does not bias anything. The cis–trans interconversion barrier is quite low, and hence most of the time conformational transition happens on the femtosecond time scale.¹¹ This means that during a typical dissociation trajectory, the cis–trans interconversion happens many times. In addition, as we are using microcanonical sampling, the initial energy is distributed randomly among all normal modes. Consequently, these low-barrier conformational transitions may occur even in the very first step of the dynamics simulation. This behaviour is a natural outcome of the sampling process (see SI for the mathematical formulation). The sampling energies were chosen as 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 20, 30, 40, and 50 kcal mol⁻¹. It is worth mentioning that the threshold energy for HO₃˙ dissociation is ∼3.8 kcal mol⁻¹ in the present PES. This means that the lowest energy at which dissociation can happen is ∼3.8 kcal mol⁻¹. This suggests that the minimum energy for our trajectory simulation should be greater than ∼3.8 kcal mol⁻¹. We have taken the lowest energy for the simulation at 5 kcal mol⁻¹ rather than 3.8 kcal mol⁻¹ for practical reasons. First, the velocity Verlet integrator was used for integration that inherently introduces small energy conservation errors of ∼ ± 0.5 kcal mol⁻¹. Second, usually during the simulation it is very unlikely that 100% of the total energy goes into the dissociation mode and hence to make dissociation feasible, one needs energy greater than the threshold. The rotational temperature was set at 300 K. At each energy, 5000 trajectories were propagated with a time step of 0.5 fs. The time at which the HO–O₂ bond distance reaches 6 Å is defined as the lifetime of the trajectory.

3. Results and discussion

3.1. Rationalization

As mentioned in the introduction, the previous on-the-fly dynamics study clearly indicates the presence of the non-statistical effects in the dissociation dynamics of HO₃˙. This hypothesis is further confirmed by our present dynamics simulations. To demonstrate this, we have illustrated the relative population of HO₃˙ remaining (N(t)/N(0)) as a function of time, along with the corresponding ln[N(t)/N(0)] vs. time plots, at different microcanonical sampling energies in Fig. 2. Here, N(t) represents the number of trajectories in which the molecule remains undissociated at time t, while N(0) is the total number of initial trajectories. It is important to note that at all sampling energies, 100% of the trajectories result in complete dissociation into products. According to RRKM theory, which assumes ergodic redistribution of energy among all available degrees of freedom, the N(t)/N(0) vs. time plot should exhibit a mono-exponential decay. Therefore, deviation from this, i.e., a multi-exponential decay, indicates the presence of non-RRKM (i.e., non-statistical) dynamics.^28–30 Although our trajectory simulations were performed over a broad range of sampling energies (5–50 kcal mol⁻¹), we have selected three representative microcanonical sampling energies (5, 10, and 20 kcal mol⁻¹) to depict N(t)/N(0) vs. time behaviour (time-dependent relative population plots at other energies are shown in Figure S2 and S3 of the SI). These three energies correspond to low, intermediate, and high energy regimes. It can be seen from Fig. 2 that, at each of these energies, the N(t)/N(0) profile requires a bi-exponential fitting function of the form f₁e^−k₁t + f₂e^−k₂t to represent the dynamics. This function includes two decay components; a very fast one and a very slow one, determined by the rate constants k₁ and k₂, respectively. The relative magnitudes of coefficients f₁ and f₂ serve as an indicator of non-RRKM behaviour. Specifically, a larger value of f₂, which is associated with the slower decay or higher-lifetime trajectories, suggests a stronger deviation from RRKM dynamics. It is evident from Fig. 2 that the contribution of f₂ increases as we move from high to intermediate to low sampling energies. For example, at 20 kcal mol⁻¹, f₂ is ∼0.33, which increases to ∼0.37 at 10 kcal mol⁻¹, and further rises substantially to ∼0.65 at the lowest sampling energy of 5 kcal mol⁻¹. This trend suggests that lower sampling energies result in a larger fraction of long-lived trajectories.

Fig. 2 Relative population of HO₃˙, N(t)/N(0), as a function of time, along with the corresponding ln[N(t)/N(0)] vs. time plots, at microcanonical sampling energies of E = 5, 10, and 20 kcal mol⁻¹.

To better understand this behaviour, a phase-space bottleneck representation of the biexponential fitting is shown in Fig. 3. The presence of two exponential terms in N(t)/N(0) decay implies a partitioning of the phase space into two distinct regions, which indicates that the vibrational modes of HO₃˙ can be divided into two groups: one leading to fast decay and the other contributing to slow decay. At lower sampling energies, the majority of trajectories fall into the slow decay category, whereas at higher energies, most trajectories follow the fast decay path. One possible reason for this behaviour can be the slower intramolecular vibrational energy redistribution (IVR) at lower energies, which hinders smooth energy flow among vibrational modes.

Fig. 3 Phase space bottleneck representation of biexponential fittings of N(t)/N(0) at excitation energies of E = 5, 10, and 20 kcal mol⁻¹.

Fig. 2 also contains the average lifetime (τ) of the trajectories at each energy. It is evident from Fig. 2 that the average lifetime increases dramatically with decreasing sampling energy. For example, at 20 kcal mol⁻¹, the average lifetime is ∼6 ps, whereas at 5 kcal mol⁻¹, it reaches the nanosecond timescale (∼1.1 ns). This finding is consistent with previous on-the-fly predictions, which suggest that the lifetimes should increase substantially at lower sampling energies.¹¹

3.2. Prediction

The experimental lifetime was measured at a particular temperature (298 K), not energy. Therefore, to connect our results to experimental lifetime, we need rate constants (inverse of lifetime) as a function of temperature rather than energy, i.e., we need to switch from microcanonical rate constant to canonical rate constant. This transformation can be achieved using the following equation:

In this equation, E₀ is the threshold energy below which HO₃˙ dissociation does not take place. k(T) is the canonical rate constant at temperature T, and it is calculated by averaging the microcanonical rate constants k(E), weighted by Boltzmann distribution.³¹ We define k(E) as the inverse of the average lifetimes at each sampling energy. To perform the integration, we need a smooth, continuous function that represents k(E) across all relevant energies not just at the sampled points. In other words, this continuous function allows us to estimate k(E) at any energy value needed for the integration. It is important to mention that one single function for k(E) does not fit well across the entire energy range, because the system behaves differently at low, intermediate, and high energies. Therefore, we have divided the energy range into three parts; low (5–7 kcal mol⁻¹), intermediate (8–12 kcal mol⁻¹), and high (13–50 kcal mol⁻¹), and fitted k(E) separately in each region. In the low-energy region, RRK (Rice–Ramsperger–Kassel) expression is commonly used to fit the energy dependence of classical microcanonical unimolecular rate constants:^32,33

Here, ν is the frequency factor, E₀ is the threshold energy required for the reaction (taken as 3.8 kcal mol⁻¹ in this work), and s represents the effective number of uncoupled oscillators. Both ν and s are treated as fitting parameters. The RRK expression is particularly suitable for the low-energy regime for two main reasons. First, it offers a convenient analytical form for describing unimolecular rates. Second, at lower energies, the impact of anharmonicity is minimal, making the harmonic assumption behind the RRK model more valid. In the high-energy region, the rate of IVR becomes comparable to the unimolecular dissociation rate, and hence, an IVR-based expression³² is sufficient to describe and fit the energy dependence of the rate constant in this region:

k_IVR(E) = AEⁿ + B (3)

Here, A, n, and B are the fitting parameters. In the intermediate energy region, both RRK behaviour (seen at low energies) and IVR effects (important at high energies) play a role. To account for both, we use a diffusional model³² given by the following equation:Here, k_RRK(E) is the same as defined in eqn (2), and k_IVR(E) is a simple expression of the form AEⁿ, where A and n are fitting parameters. This model combines both rates in a way that smoothly connects the low-energy RRK behaviour with the high-energy IVR behaviour. In simple terms, this model gives a good overall fit in the intermediate energy range, where neither RRK nor IVR alone was enough to describe the system accurately. In principle, diffusional expression can also be used to fit the entire energy range. However, we used separate fittings to reduce the fitting error, particularly in light of the limited data available (5–50 kcal⁻¹). Fig. S6 contains the mean unsigned deviation (MUD) for the separate fittings, along with the MUD obtained from fitting the entire dataset using the diffusional energy expression. Figure S6 clearly shows that the MUD obtained with the diffusional model over the full range is comparatively large (MUD ∼ 14%) compared to the same obtained when the data are fitted in three separate regions (MUD ∼ 4.5%).

Using eqn (2)–(4), we have fitted the energy-dependent rate constants k(E) within low, intermediate, and high energy regions. The results of these fits are shown in Fig. 4. It is evident from Fig. 4 that each expression provides a good fit within its respective energy range. The maximum mean unsigned deviation (MUD) is found to be ∼12% for RRK expression. This level of uncertainty corresponds to ∼0.1 μs for 1 μs lifetime, which is reasonable for the purposes of this study. Finally, to compute the temperature-dependent rate constant k(T), we have used a piecewise integration approach based on fitted k(E) expressions in their respective energy ranges, i.e., k(T) is calculated using the following equation:

Here, , where k_B is Boltzmann constant and T is temperature. E₀, E₁, E₂, and E_max define the integration limits of the RRK, diffusional, and IVR regions, respectively. Quantitatively, E₀ is the threshold energy (3.8 kcal mol⁻¹) below which the reaction probability is zero. The value of E₁ is 7 kcal mol⁻¹, and E₂ is 12 kcal mol⁻¹. Although our maximum sampling energy was 50 kcal mol⁻¹, we set E_max to 500 kcal mol⁻¹ to make sure the integral in the denominator converges properly.

Fig. 4 Fitting of the logarithm of the microcanonical rate constants, k(E), separately in the low, intermediate, and high-energy regimes using RRK, diffusional, and IVR-based energy-dependent expressions, along with the fitting of k(E) in the entire energy range using diffusional energy-dependent expression.

The calculated lifetimes are depicted in Fig. 5 for the temperature range of 200–500 K. Fig. 5 clearly shows that the predicted lifetime at 300 K is ∼1.8 μs, which agrees well with the experimental recommendation of greater than 1 μs. Fig. 5 also shows that the lifetime increases drastically as the temperature decreases. For example, the lifetime is ∼0.2 μs at 400 K, which increases to ∼1.8 μs at 300 K, and becomes ∼116 μs at 200 K. These results suggest that at temperatures ≤200 K, the lifetime of HO₃˙ can enter the millisecond timescale. This means that at low temperatures and under collision-free conditions, HO₃˙ is indeed a stable intermediate.

Fig. 5 Average lifetimes estimated using canonical rate constants k(T) in the temperature range of 200–500 K.

After establishing that HO₃˙ is a stable transient species, it is important to discuss the limitations of this study and quantify the uncertainty associated with the present computation of the lifetime of HO₃˙. One approximation in the present work is the negligence of tunneling effects, as microcanonical classical dynamics simulations do not account for quantum tunneling. To understand the role of tunneling, we have calculated the crossover temperature³⁴ (T_c), which is generally used as an indicator of the temperature below which tunneling becomes important:

Here, ν is the imaginary frequency associated with the transition state, h is Planck's constant, and k_B is the Boltzmann constant. The two main bottlenecks that can contribute to the lifetime of HO₃˙ via tunneling are cis–trans isomerization, and the direct unimolecular dissociation path. The potential energy surface containing both of these pathways was obtained at the B3LYP/aug-cc-pVDZ level of theory (shown in Fig. S1 of the SI). This method was chosen due to its good agreement with benchmark CCSDTQ(P)/CBS energetics. Firstly, we have estimated the crossover temperatures (T_c) for both pathways to assess the possible role of tunneling. The calculated T_c turns out to be ∼67 K for HO₃˙ dissociation and ∼40 K for cis–trans isomerization. These values suggest that tunneling may become crucial only below ∼70 K. Secondly, we have also estimated Eckart tunneling corrections for both processes within 200–400 K (Table S1 of the SI). It is evident from Table S1 that as expected, the tunneling contributions are negligible for both the pathways. Therefore, we believe that tunneling is not going to affect the estimated lifetimes of HO₃˙ within the temperature range (200–400 K) considered in the present work.

Another factor that may influence our results is the uncertainty in the threshold energies of three crucial phenomena in our system. As mentioned earlier, the PES employed in this study has a threshold energy of ∼3.8 kcal mol⁻¹ for unimolecular dissociation from cis-HO₃˙, whereas the experimentally measured binding energy is ∼3.0 kcal mol⁻¹.⁹ Similarly, the cis–trans and trans–cis isomerization thresholds are ∼0.9 and 1.1 kcal mol⁻¹, respectively, whereas the benchmark CCSDTQ(P)/CBS results suggest that they should be 0.3 and 0.8 kcal mol⁻¹, respectively. To assess the impact of this discrepancy, we employed a two-state coupled phase space model to estimate the lifetime in the temperature range of 200–500 K. The model consists of two distinct phase space regions (see Fig. 6). The first region in Fig. 6 represents cis/trans-HO₃˙, which is directly connected to the products, i.e., OH˙ + O₂, while the second region represents trajectories that remain trapped in a long-lived intermediate due to rapid cis–trans isomerization. Although products can form from either the cis or trans conformer, we are using the two-state model (assuming product formation only from cis-HO₃˙) to estimate the effect of threshold energy for two practical reasons. The first reason is that most of the trajectories dissociate from cis because it has a slightly lower dissociation energy compared to trans. Second, if one includes the product formation channel from trans, the system becomes a three-state model, for which analytical solutions are not available. On the other hand, our goal of the kinetic model is to study the effect of threshold energy for which the cis conformer alone is sufficient. The kinetic scheme for the two-state model is given as:

k₁ is the rate of product formation, while k₂ and k₃ represent isomerization rates for the cis-to-trans and trans-to-cis conversions, respectively. Assuming an initial microcanonical ensemble, the effective rate constant of dissociation in this model is:

Fig. 6 A schematic of the two-state coupled phase space model of HO₃˙.

All three thresholds that can affect the lifetime of HO₃˙ are present in this model. First, the energy needed for the cis-HO₃˙ molecule to dissociate directly is present in k₁. Second, the energy required for cis-HO₃˙ to convert into trans form is included in k₂. Third, the energy for trans-HO₃˙ to convert into cis is contained in k₃. This model enables us to adjust each of these energy thresholds, making it easier to understand how each one affects the overall reaction rate. Let N₁(t) and N₂(t) denote the population of the trajectories in region I and region II, respectively. Under the initial conditions, i.e., N₁(0) = N(0), and N₂(0) = 0, solving the coupled differential equations corresponding to the kinetic scheme gives an analytical expression for the time-dependent dissociation probability from region I:

Here, λ₁ + λ₂ = k₁ + k₂ + k₃, and λ₁λ₂ = k₁k₃. The values of λ₁, λ₂, and k₁ can be obtained by comparing eqn (10) with our bi-exponential fit of the form f₁e^−λ₁t + f₂e^−λ₂t. Once these values are obtained, they are used to calculate k₂, k₃, and k₀(E). The resulting k₀(E) values are again fitted across the three energy regions using eqn (2)–(4). The fitted rate constants k₀(E) across the three energy regions are shown in Fig. S4 of the SI. Here also, we found a good fitting with maximum MUD of 7%. The canonical rate constants k(T) are then calculated using eqn (5). Fig. 7 compares the lifetimes estimated using both the methods, i.e., earlier average rate constant approach and two state model approach. Quantitatively, the lifetime calculated using the average rate constant approach is ∼1.8 μs, while the two-state model approach gives a lifetime of ∼1.5 μs, only ∼0.3 μs shorter than the average approach. This close agreement confirms that quantifying the effect of the threshold energies using the two-state model is reasonable.

Fig. 7 Comparison of lifetimes obtained using the two-state model approach, the average approach, and the threshold-corrected model approach in the temperature range of 200–500 K.

To incorporate the change in threshold energy, we first fitted k_1(E), k_2(E), and k_3(E) separately using the three energy-dependent expressions (eqn (2)–(4)) corresponding to the low, intermediate, and high energy regions. The fitted rate constants k₁(E), k₂(E), and k₃(E) across the three energy regions are shown in Fig. S5 of the SI. Then, k₀(E) is calculated using eqn (9), followed by the calculation of the canonical rate constants using eqn (5). In the next step, using the same fitted parameters for k₁, k₂, and k₃, we replaced the threshold energies of PES with the same computed at the CCSDTQ(P)/CBS level of theory and the experimental value. Specifically, the dissociation energy of ∼3.8 kcal mol⁻¹ was replaced with ∼3.0 kcal mol⁻¹, which corresponds to the experimental value. Similarly, the cis–trans isomerization thresholds of ∼0.9 and ∼1.1 kcal mol⁻¹ were replaced with ∼0.3 and ∼0.8 kcal mol⁻¹, respectively, which are the benchmark computational values. As discussed in Section 2, this benchmark-level PES provides the most accurate threshold energies available to date for HO₃˙.¹⁰ This procedure was applied over the 200–500 K temperature range. The resulting threshold correction factors (defined as the ratio of rate constants calculated using the PES value and benchmark value) were then multiplied by the average lifetimes across this temperature range to obtain the threshold-corrected lifetimes. The average lifetimes, two-state model lifetimes, and threshold-corrected lifetimes, along with the corresponding threshold correction factors (TCFs), are provided in Table S2 of the SI. The same results are also illustrated in Fig. 7. It can be seen from Fig. 7 that the lifetime does not change much at all temperatures when the threshold energy corrections are incorporated. For example, at 300 K, the average lifetime is ∼1.8 μs, and the threshold-corrected lifetime is ∼1.9 μs, while at 200 K, the average lifetime is ∼116 μs, and the threshold-corrected lifetime increases only slightly (to ∼117 μs, see Table S2). This suggests that the lowering of the HO₃˙ unimolecular dissociation threshold from ∼3.8 to ∼3.0 kcal mol⁻¹ is nullified by the simultaneous lowering of the other bottleneck, i.e., the cis–trans isomerization barrier, from ∼0.9 to ∼0.3 kcal mol⁻¹. Therefore, it is the cis–trans isomerization barrier that regulates the dissociation dynamics of HO₃˙. In summary, the effect of threshold energy does not alter our conclusion that HO₃˙ is a stable intermediate. At 300 K, a more reasonable lifetime incorporating all possible effects is ∼1.8 ± 0.1 μs, which is in good agreement with experimental recommendations.

In addition, it is important to mention that, in a previous study, Varner et al.³⁵ found a small exit-channel barrier (∼3 kcal mol⁻¹) along the trans-HO₃˙ → OH˙ + O₂ minimum-energy path (computed at the EOMIP-CCSD*/cc-pVQZ level of theory). In fact, this was key in explaining the discrepancy between theory and early experiments conducted before 2010, where the dissociation energy was found to be ≤5.3 kcal mol⁻¹.⁸ In the present PES, we have not found any exit barrier from trans-HO₃˙. One reason could be the level of theory used in the present PES, i.e., MRCI(Q)-F12/VDZ-F12 compared to EOMIP-CCSD*/cc-pVQZ used by Varner et al. Usually, in bond dissociation processes, dynamic correlation plays a key role. Although the EOMIP-CCSD*/cc-pVQZ level of theory does account for some dynamic correlation, MRCI(Q)-F12/VDZ-F12 is certainly more accurate in dealing with bond dissociation processes. Nevertheless, if this exit barrier does exist in the real system, we believe it may slightly increase the lifetime of HO₃˙.

Finally, it is important to mention that the lifetime obtained from unimolecular dissociation is an upper limit. In the real atmosphere, there will always be bimolecular collisions that will further reduce the lifetime. Consequently, these results are more relevant to the upper atmosphere, where both pressure and density are low and the probability of bimolecular collisions is minimal. In contrast, in the lower atmosphere, where pressure and density are higher, fast bimolecular relaxation can further reduce the lifetime.

4. Conclusion

In the present work, using classical dynamics simulations, we have shown that at 300 K, the average lifetime of HO₃˙ is ∼1.8 μs, which agrees well with the experimental recommendation of greater than 1 μs. We have also found that the lifetime increases dramatically as the temperature decreases, reaching the millisecond timescale below 200 K. This suggests that under low-temperature and collision-free conditions, HO₃˙ is indeed a stable intermediate. Finally, we have discussed the limitations of our results, which can arise due to tunneling and uncertainty in threshold energy. The effect of tunneling was negligible, and at the same time, uncertainty due to threshold energy does not change much the lifetime of HO₃˙ at all temperatures. This suggests that it is not the unimolecular dissociation threshold which only decides the fate of HO₃˙; rather, it is the combined effect of both the bottlenecks, i.e., unimolecular dissociation as well as the cis–trans isomerization barrier, which governs the lifetime of HO₃˙.

Author contributions

PKR: conducted the investigation, writing – original draft, formal analysis, curated the data. AG: helps in coupling neural-network potential with VENUS code. PK: provided supervision, resources, and methodology; conceptualized the study; acquired funding; and contributed to the review and editing of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data available within the article or its SI. Supplementary information: Eckart tunneling corrections, Average lifetimes, two-state model lifetimes, and threshold-corrected lifetimes along with threshold-correction factors (tcf) within 200–500 K, relative population of HO₃˙ remaining, N(t)/N(0) vs. time, potential energy surface for HO₃˙, fitting of microcanonical rate constants k₀(E), fitting of rate constants k₁(E), k₂(E), and k₃(E), and mathematical formulation of microcanonical sampling. See DOI: https://doi.org/10.1039/d5cp02134k.

Acknowledgements

P. K. R. acknowledges MNIT Jaipur for financial assistance. A. G. acknowledges the PMRF scheme, Govt. of India, for the financial support, and P. K. acknowledges DST, Govt. of India, for the financial support through the sanctioned project No. EEQ/2023/000351.

References

1. M. Speranza, Structure, Stability, and Reactivity of Cationic Hydrogen Trioxides and Thermochemistry of their Neutral Analogs. A Fourier-Transform Ion Cyclotron Resonance Study, Inorg. Chem., 1996, 35, 6140–6151.
2. F. Cacace, G. De Petris, F. Pepi and A. Troiani, Experimental Detection of Hydrogen Trioxide, Science, 1999, 285, 81–82.
3. U. Sridharan, F. Klein and F. Kaufman, Detailed course of the O + HO₂ reaction, J. Chem. Phys., 1985, 82, 592–593.
4. K. Suma, Y. Sumiyoshi and Y. Endo, The Rotational Spectrum and Structure of the HOOO Radical, Science, 2005, 308, 1885–1886.
5. C. Murray, E. L. Derro, T. D. Sechler and M. I. Lester, Stability of the Hydrogen Trioxy Radical via Infrared Action Spectroscopy, J. Phys. Chem. A, 2007, 111, 4727–4730.
6. E. L. Derro, C. Murray, T. D. Sechler and M. I. Lester, Infrared Action Spectroscopy and Dissociation Dynamics of the HOOO Radical, J. Phys. Chem. A, 2007, 111, 11592–11601.
7. E. L. Derro, T. D. Sechler, C. Murray and M. I. Lester, Observation of ν₁ + ν_n combination bands of the HOOO and DOOO radicals using infrared action spectroscopy, J. Chem. Phys., 2008, 128, 244313.
8. C. Murray, E. L. Derro, T. D. Sechler and M. I. Lester, Weakly Bound Molecules in the Atmosphere: A Case Study of HOOO, Acc. Chem. Res., 2009, 42, 419–427.
9. S. D. Le Picard, M. Tizniti, A. Canosa, I. R. Sims and I. W. Smith, The Thermodynamics of the Elusive HO₃ Radical, Science, 2010, 328, 1258–1262.
10. M. A. Bartlett, A. H. Kazez, H. F. Schaefer and W. D. Allen, Riddles of the structure and vibrational dynamics of HO₃ resolved near the ab initio limit, J. Chem. Phys., 2019, 151, 094304.
11. P. K. Rai and P. Kumar, Role of non-statistical effects in deciding the fate of HO₃˙ in the atmosphere, Phys. Chem. Chem. Phys., 2024, 26, 24785–24790.
12. A. B. Meinel, OH Emission Bands in the Spectrum of the Night Sky. II, Astrophys. J., 1950, 112, 120.
13. K. S. Kalogerakis, A previously unrecognized source of the O2 Atmospheric band emission in Earth's nightglow, Sci. Adv., 2019, 5, eaau9255.
14. K. J. Rensberger, J. B. Jeffries and D. R. Crosley, Vibrational relaxation of OH (X²Π_i, ν = 2), J. Chem. Phys., 1989, 90, 2174–2181.
15. J. A. Dodd, S. J. Lipson and W. A. Blumberg, Vibrational relaxation of OH (X²Π_i, ν = 1–3) by O₂, J. Chem. Phys., 1990, 92, 3387–3393.
16. L. D'Ottone, D. Bauer, P. Campuzano-Jost, M. Fardy and A. J. Hynes, Vibrational deactivation studies of OH X²Π (ν = 1–5) by N₂ and O₂, Phys. Chem. Chem. Phys., 2004, 6, 4276–4282.
17. D. McCabe, I. Smith, B. Rajakumar and A. Ravishankara, Rate coefficients for the relaxation of OH (ν = 1) by O₂ at temperatures from 204–371 K and by N₂O from 243–372 K, Chem. Phys. Lett., 2006, 421, 111–117.
18. P. K. Rai and P. Kumar, HO₂˙ + O₃ → OH˙ + 2O₂ reaction: A potential source of vibrationally hot OH radicals in the atmosphere, Int. J. Chem. Kinet., 2023, 55, 619–628.
19. P. K. Rai and P. Kumar, Accurate determination of reaction energetics and kinetics of HO₂˙ +O₃ → OH˙ +2O₂ reaction, Phys. Chem. Chem. Phys., 2023, 25, 8153–8160.
20. H. Yu and A. Varandas, Ab initio theoretical calculation and potential energy surface for ground-state HO₃, Chem. Phys. Lett., 2001, 334, 173–178.
21. B. J. Braams and H.-G. Yu, Potential energy surface and quantum dynamics study of rovibrational states for HO₃ (X²A′′), Phys. Chem. Chem. Phys., 2008, 10, 3150–3155.
22. X. Hu, J. Zuo, C. Xie, R. Dawes, H. Guo and D. Xie, An ab initio based full-dimensional potential energy surface for OH + O₂ HO₃ and low-lying vibrational levels of HO₃, Phys. Chem. Chem. Phys., 2019, 21, 13766–13775.
23. J. Zuo, Q. Chen, X. Hu, H. Guo and D. Xie, Theoretical Investigations of Rate Coefficients for H + O₃ and HO₂ + O Reactions on a Full-Dimensional Potential Energy Surface, J. Phys. Chem. A, 2020, 124, 6427–6437.
24. W. L. Hase and D. G. Buckowski, Monte Carlo sampling of a microcanonical ensemble of classical harmonic oscillators, Chem. Phys. Lett., 1980, 74, 284–287.
25. G. H. Peslherbe, H. Wang and W. L. Hase, Monte Carlo sampling for classical trajectory simulations, Adv. Chem. Phys., 1999, 105, 171–201.
26. X. Hu, W. L. Hase and T. Pirraglia, Vectorization of the general Monte Carlo classical trajectory program VENUS, J. Comput. Chem., 1991, 12, 1014–1024.
27. W. L. Hase, R. J. Duchovic, X. Hu, A. Komornicki, K. F. Lim, D.-H. Lu, G. H. Peslherbe, K. N. Swamy, S. Linde and A. Varandas, et al., A general chemical dynamics computer program. Quantum Chem, Program Exch. Bull, 1996, 16, 671.
28. M. Paranjothy, R. Sun, A. Kumar Paul and W. L. Hase, Models for Intrinsic Non-RRKM Dynamics. Decomposition of the S_N2 Intermediate Cl⁻–CH₃Br, Z. Phys. Chem., 2013, 227, 1361–1379.
29. A. Kumar Paul, S. Kolakkandy, S. Pratihar and W. L. Hase, Reaction Rate Constant Computations: Theories and Applications, The Royal Society of Chemistry, 2013.
30. S. Malpathak, X. Ma and W. L. Hase, Direct dynamics simulations of the unimolecular dissociation of dioxetane: Probing the non-RRKM dynamics, J. Chem. Phys., 2018, 148, 164309.
31. J. I. Steinfeld, J. S. Francisco and W. L. Hase, Chemical Kinetics and Dynamics, Prentice Hall, Upper Saddle River, NJ, 1999, vol. 26.
32. D. V. Shalashilin and D. L. Thompson, Intrinsic non-RRK behavior: Classical trajectory, statistical theory, and diffusional theory studies of a unimolecular reaction, J. Chem. Phys., 1996, 105, 1833–1845.
33. K. Song and W. L. Hase, Fitting classical microcanonical unimolecular rate constants to a modified RRK expression: Anharmonic and variational effects, J. Chem. Phys., 1999, 110, 6198–6207.
34. M. Gillan, Quantum-classical crossover of the transition rate in the damped double well, J. Phys. C: Solid State Phys, 1987, 20, 3621.
35. M. E. Varner, M. E. Harding, J. Vázquez, J. Gauss and J. F. Stanton, Dissociation Energy of the HOOO Radical, J. Phys. Chem. A, 2009, 113, 11238–11241.

---