Decomposition Kinetics for HONO and HNO 2

This work presents a detailed investigation into the isomerization and decomposition of HONO and HNO 2 . State-of-the-art electronic structure theory is used to compute the HNO 2 potential energy surface. Temperature and pressure dependent rate coeﬃcients are computed using microcanonical rate theory and the master equation. The electronic structure theory properties are optimized against the relevant experimental data. A novel strategy was developed to incorporate uncertainty in the minimum energy pathway into the optimized mechanism. The new mechanism is in excellent agreement with all available experimental data for H + NO 2 → OH + NO and OH + NO → HONO. The calculations identify OH + NO as the dominant products for HNO 2 , which were neglected from all previous mechanisms in the literature


Introduction
Nitrous acid (HONO) and nitryl hydride (HNO 2 ) are key intermediates in many gas-phase systems, from atmospheric chemistry to combustion to energetic materials. The most common source of these two species is through H-abstraction or disproportionation involving NO 2 . Additionally, HONO can be formed via concerted elimination from various NO 2 -containing functional groups, such as nitroalkanes (C−NO 2 ), alkylnitrates (O−NO 2 ), and nitroamines (N−NO 2 ). In contrast, concerted HNO 2 elimination from the analogous nitrite compounds, X−ONO, does not exist as a viable product channel, owing to the enthalpically and entropically more favorable O−NO fission.
Recent developments in advanced engines, such as exhaust gas recirculation (EGR) and the use of alkyl-nitrates as reactivity enhancers, create new environments in which NO 2 chemistry is increasingly important [1][2][3].
One way in which NO 2 contributes to increased reactivity (e.g. decreased ignition delay times) is through various radical + NO 2 reactions that shift the composition from comparatively less reactive radicals to more reactive radicals [4][5][6]: Whereas the reaction sequence (R1-R4), is chain propagating, the reaction sequence (R5-R8) is chain branching. Recent experimental and modeling work has highlighted the sensitivity of ignition delay times to this reaction sequence for H 2 , CH 4 , and C 3 H 8 [7][8][9][10][11]. One of the authors recently published rate coefficients for (R5) obtained from transition state theory (TST) calculations for a broach range of alkanes and alkenes. [12] The HNO 2 potential energy surface (PES), shown in Figure 1.
As seen in Figure 1, two different rotational conformers of HONO exist. The more stable rotamer, anti-HONO (frequently referred to as "trans- HNO 2 potential energy surface and the kinetics of reactions (R4, R6-R8).
The approach will be to use state-of-the-art methods in electronic structure theory and computational kinetics to provide temperature-and pressuredependent rate coefficients. The electronic structure and collisional energy transfer parameters are then optimized against diverse experimental data to provide a rigorously optimized system of rate constants.
As depicted in Figure 1, the H + NO 2 and the OH + NO pathways do not have first-order saddle points in potential energy. The minimum energy paths that are required for the variational transition state theory were treated using two different multi-reference methods CASPT2 [14][15][16] and MRCI+Q [17][18][19]. The active space in the H + NO 2 calculations was 14 electrons in 11 orbitals (14e,11o): 6 orbitals for the delocalized π-system of NO 2 , 4 orbitals for the two pairs of N−O σ, σ * in NO 2 , and 1 orbital for the H atom (this active space is effectively equivalent to the complete valence space for the system, not including the two 2s orbitals for the oxygen atoms). The CASPT2 and MRCI+Q calculations were performed with aug-cc-pVQZ and aug-cc-PV5Z basis sets and extrapolated to the complete basis set limit. To further improve the accuracy of the multi-reference calculations, a separate set of potentials was computed on a higher spin surface [20][21][22][23]. The potential energy of the triplet surface was computed using both the ANL1 method of Equation (1) as well as the CASPT2 and MRCI+Q approaches. The coupled-cluster triplet surface is then added to the multi-reference singlettriple splitting to obtain a more accurate representation of the minimum energy path on the singlet surface: where the superscript indicates the spin multiplicity of the surface, and "MR" implies either CASPT2 or MRCI+Q with aug-cc-pV∞Z basis sets.
The H + NO 2 pathway leading to anti-HONO, despite being a radicalradical reaction, has a well-defined first-order saddle point, depicted by the dotted line in Figure 1. This transition still has a strong bi-radical character, and multi-reference methods were necessary. The transition state optimization was performed using CASPT2(14e,11o)/aug-cc-pVQZ, with the same active space as above. The spin-splitting method of (2) was applied to further improve the accuracy. This local maximum is a consequence of the nodal planes in the NO 2 singularly occupied molecular orbital (SOMO); as the H atom approaches either O atom in the anti configuration, it experiences significant overlap with the lobes of opposite sign on the N atom on syn configuration of the O atoms [12].
For the OH + NO potentials, the active space was (8e,6o): 4 orbitals for the complete π system in NO and 2 orbitals for the π system in OH. The orbitals in these calculations were averaged over four states to account for the spatial degeneracy in both OH and NO on the singlet surface. Because of the low-lying excited states, single determinant methods were not reliable, and the singlet-triplet splitting approach in Equation (2) could not be applied.
All wavefunction-based calculations were performed in Molpro. [24] The density functional theory calculations used for the anharmonic correction were done in Gaussian09 [25].

Computational Kinetics
Reactions (R4) and (R6)-(R8) proceed through one or more unimolecular intermediates, and thus the rate constants are a function of pressure. The temperature-and pressure-dependent rate constants are computing using the microcanonical rate theory and the master equation (RRKM/ME) code Mess [26,27], which is part of the computational kinetics package Papr developed by Argonne National Laboratory [28]. In the time-dependent master equation, the eigenvalues of the transition matrix can be separated into chemically significant eigenvalues (CSE) and internal energy relaxation eigenvalues (IERE), with the former typically smaller than the latter by several orders of magnitude [30,31]. When a CSE approaches the smallest IERE, it indicates that a pair of wells (or a well and bimolecular products) will equilibrate on a timescale that approaches the rate of ro-vibrational energy relaxation. When this merging of timescales occurs, it no longer makes sense to think of these species as being chemically distinct in a macroscopic model [32]. An important feature of Mess is the ability to determine when two wells equilibrate rapidly (e.g. at high temperature) and combine them accordingly [26].
This "well-merging" capability is important for the HNO 2 reactive system. At low temperatures, the two HONO rotamers remain chemically distinct. At higher temperatures, however, the chemically significant eigenvalue that corresponds to the interconversion between anti-HONO and syn-HONO blends into the IERE continuum. Preliminary calculations suggest that this merging occurs at 650 K and 0.01 bar, and at 800 K and 100 bar. Above these temperatures, they exist only as a single HONO. Below these temperatures, when the two rotamers are distinct, the rate constant for their interconversion is large (e.g. 2 × 10 9 s −1 ), and they will rapidly equilibrate in the conventional sense. Precisely how species merging should be handled at high temperatures remains an open challenge in reactive flow simulations, since it represents a reduction in the number of species (and accordingly the state vector) [33]. Since our primary interest is in high-temperature kinetics, we lump these two species into a single HONO isomer. In this two-well model, with HONO and HNO 2 as distinct isomers, the interconversion between anti- are relaxed, V relaxed (r). The difference in these two V (r) is the 1D geometry relaxation correction. Additionally, single-point calculations are performed using larger basis sets for the fully relaxed case, which can then be used to provide a basis set correction. As will be detailed below, these 1D corrections can be manipulated to explore how the uncertainty in the interaction For the H + NO 2 reactions, the H atom can add to either the N or the O atoms. Within the VRC-TST methodology, the flux to these different products (HNO 2 and HONO, respectively) must be determined separately.
Accordingly, a plane of infinite potential was imposed on the NO 2 , normal to the plane [39,41]. For H + NO 2 → HNO 2 , the plane is defined such that bonding to the O atoms is impossible, and for H + NO 2 → HONO, the plane is in the same location, but with a sign change, such that bonding to the N atom is impossible. The location of this dividing surface was optimized variationally so as to minimize the total H + NO 2 flux.
The rate constants are computed at temperatures 200 ≤ T ≤ 2500 K and pressures 0.01 ≤ P ≤ 100 bar.

Optimization Methodology
A key aspect of the present work is the inclusion of uncertainty in each of the model parameters. Each parameter can be adjusted within its upper and lower bound as a means of providing global uncertainty quantification.
Moreover, this ensemble of models can then be compared to diverse experimental data to determine the optimum set of model parameters that best reproduce the data. This form of optimization is similar to the Multiscale Informatics (MSI) of Burke and coworkers [42][43][44][45]. Previous implementations of the MSI methodology focused on stationary points on the PES, as well as vibrational frequencies and energy transfer models. The present work is the first extension of the MSI framework to include the potential energy surface for barrierless reactions. The MSI variables are the eight stationary points in Table 1, with the upper and lower limits as expected by the ANL1 method; the collisional energy transfer pre-factor, 100 ≤ ∆E down ≤ 500 cm −1 , and the scaling factors, f , for H + NO 2 → HNO 2 , H + NO 2 → HONO, and OH + NO → HONO.
Two sets of experimental data are used to constrain the model parameters. The first set of experiments are for the H + NO 2 → OH + NO reaction, which includes flow cell data at lower temperatures [46], high-pressure flow cells at intermediate temperatures [47], and shock tube measurements using two different diagnostics at high temperatures [48]. These data primarily constrain the two H + NO 2 interaction potentials. The second set of experimental data are the high-pressure flow cell work by Troe and coworkers on OH + NO → HONO [49]. These data primarily constrain the OH + NO interaction potential and the collisional energy transfer coefficient, ∆E down .  Table 1: key stationary points on the HNO 2 potential energy surface, relative to H + NO 2 .

Results and Discussion
The ANL1 results for the stationary points on the PES are presented in Table 1. Also included in the table are the ATcT results. The root mean square deviation between the two methods was 0.7 kJ/mol. The results of the H + NO 2 → HNO 2 potential are shown in Figure 3.

H + NO
The four colors correspond to the four different methods that were used to compute the minimum energy pathway: pure CASPT2(14e,11o), cyan; pure MRCI+Q(14e,11o), magenta; ANL1 on the triplet surface, with CASPT2(14e,11o) singlet-triplet splitting, blue; and ANL1 on the triplet surface, with MRCI+Q(14e,11o) singlet-triplet splitting, red. Each of these four methods was computed using both aug-cc-pVQZ and aug-cc-pV5Z basis sets and then extrapolated to the CBS limit.
As can be seen in the figure, all four methods are well converged with respect to the basis set and with respect to each other. Figure 3b plots the difference between each of the four methods and the geometric mean of the four (at aug-cc-pV∞Z). Of greater significance is the percent deviation of the methods from the average, Figure 3c, since the percent deviation is a better reflection on the final uncertainty of the transition state. Note that all twelve curves agree to within 5%. In principle, it is difficult to say exactly what the "true" value for the potential is. However, given the degree of convergence, we estimate that the average value for the minimum energy path is likely within 25% of the correct value.   HONO; the blue-shaded region represents an uncertainty envelope of 25%, and the yellow-shaded region represents 50% uncertainty in the minimum energy path. The degree of convergence for this potential was the same as for the competing H + NO 2 → HNO 2 pathway, shown in Figure 3. Figure   4b illustrates the 25 V new (r) that are created from the scaling factors f in dividing surface at 300 K dividing surface at 400 K The results for the OH + NO interaction potential are shown in Figure 5.
Proper calculation of the rate constant for OH + NO → HONO is complicated by non-adiabatic effects, with four states on the singlet surface (blue curves) and four states on the triplet surface (red). For r ≥ 3.5Å, all eight curves are attractive, but only the ground state singlet remains attractive for r ≤ 2.4Å. A detailed analysis of the contribution of these higher states would require non-adiabatic dynamics simulations, which is beyond the scope of the present work. Further complications are spin-orbit effects and rovibronic coupling at lower temperatures. Instead, we note that these effects become increasingly irrelevant at higher temperatures.
To quantify the temperature at which these effects become negligible, we consider a limiting case in which the crossing from the excited electronic and spin states to the ground state singlet is rapid. In this approximation, the sum of the contributions to the interaction potential should be included prior to the variational optimization of the dividing surface. From the VRC-TST analysis, which used only the singlet ground state (solid blue line in Figure 5), we can obtain the location of the optimum dividing surface as a function of temperature. At 300 K, it is approximately r = 3.6Å, and the contribution of the other surfaces is not insignificant (nearly a factor of three). As the temperature is increased, the optimum dividing surface shifts to shorter separation distances. By 400 K, the dividing surface is located at r = 2.8Å; under these conditions, the contribution of the seven additional states is a mere 4%. This value likely is an over-estimate, since our strongcoupling assumption may be appropriate for the exited singlet states, but is less likely to be appropriate for the triplet surface. (The other limiting model is to assume that the crossing is slow, in which case the variational analysis is performed on each surface independently. Given that the other surfaces have local minima between −3.8 and −0.5 kJ/mol, their flux through the dividing surface will be comparatively negligible.) For these reasons, we do not include experimental data below 400 K in the optimization procedure, since it would bias the result towards more attractive minimum energy paths.
Instead, we restrict ourselves to the 400 K data of Fulle et al.  value (black line) corresponds to the unscaled result, which is higher than the least-squares fit by slightly less than a factor of two. As part of the optimization process, the barrier height for HNO 2 → anti-HONO was reduced by 2.5 kJ/mol from -93.9 to -96.4 kJ/mol (relative to H + NO 2 ). In terms of the scaling factors for the barrierless entrance channels: for H + NO 2 → syn-HONO, the optimized potential was 25% more attractive at close range, but unperturbed at long range (e.g. close-dash black in Figure 2); for H + NO 2 → HNO 2 , the optimized potential was unperturbed at close range, but was 25% less attractive at long range (e.g. wide-dash green in Figure   2  The optimized results for OH + NO → HONO are shown in Figure 7. The high-pressure data of Fulle [49] are an excellent compliment to the H + NO 2 data, since they provide constraints for both collisional energy transfer and OH + NO minimum energy path. The optimized result for the collisional energy transfer prefactor was ∆E down = 300 (T /298[K]) 0.85 cm −1 . The OH + NO interaction potential was shifted up by a constant of 25% (e.g. dash-dot green in Figure 2). Note that even though only the data from 400 K (black symbols) were included in the optimization process, the lower termperature data are still well reproduced by the optimized model. The optimized rate constant slightly under-predicts the experimental data at 250 and 298 K, which is consistent with our expectation that non-adiabatic effects will be important at lower temperatures.   Figure 8a, the only significant product is OH + NO, with isomerization to HNO 2 being four orders of magnitude smaller. This result is unsurprising, since, as seen in Figure 1, the OH + NO products are lower in energy than the barrier for isomerization to HNO 2 , and thus the barrierless decomposition is both enthalpically and entropically favorable at all temperatures. For HNO 2 decomposition, Figure 8b,    To assess the impact of the newly optimized HONO/HNO 2 mechanism on model predictions, we consider six different mechanisms that have been published in the past ten years: Daguat et al. [51], Konnov [52], Mathieu et al. [53], Ahmed et al. [54], Zhang et al. [10], and Glarborg et al. [6]. These mechanisms, along with the present work, were imported into Cantera. [55] In all cases, the literature mechanisms represent the reaction in the association direction, OH + NO → HONO, with values taken from either Tsang and Herron [56] (Dagaut, Konnov, Mathieu) or Fulle et al. [49] (Ahmed, Zhang, Glarborg). In some cases, the Troe broadening factor F cent [57] was modified by the authors. The largest source of discrepancy between the literature models, however, is due to the thermochemistry of OH, NO, and HONO, since the rate constant for HONO → OH + NO is computed using the equilibrium constant to maintain thermodynamic consistency. Nonetheless the agreement of the literature models with each other and with the present work is good, typically differing by less than an order of magnitude. The best agreement is with the Glarborg mechanism, which uses the ATcT thermochemistry, and is within a factor of two of the present work. As demonstrated in Figure 8, the most important product channel for HNO 2 is the well skipping reaction to form OH + NO, but none of the previous mechanisms include this reaction. Instead, they tend to rely on early predictions of Dean and Bozzelli for HNO 2 , which used QRRK and modified strong collision [5]. Importantly, that earlier work did not include OH + NO as a possible channel and thus overestimated the rate collisional stabilization of HONO from HNO 2 . The set of unimolecular rate constants are listed in

Conclusions
The HNO 2 potential energy surface is computed using post-coupled cluster methods, with an uncertainty of ± mechanism, coupled with the uncertainty analysis, is expected to provide important constraints on future mechanism development for gas-phase nitrogen chemistry.

Supplemental Material
The optimized mechanism is available in CHEMKIN and Cantera format.

Conflicts of interest
There are no conflicts to declare.