Reaction Dynamics for the [NNO] System from State-Resolved and Coarse-Grained Models

The two PESs for the ³A^′ state of the [NNO] system used in the present work were based on MRCI+Q calculations with the (m)aug-cc-pVTZ basis sets.^{18, 19} Surface PES_B (for the Basel PES)¹⁸ was represented as a reproducing kernel Hilbert space^{27, 28, 29} whereas PES_M (for the Minnesota PES)¹⁹ was fit to permutationally invariant polynomials (PIPs). For the CASSCF calculations, 10 active electrons in 9 orbitals CAS(10,9) were used. For PES_B the CASSCF wavefunction was used for the MRCI + Q calculations, while a state-averaged-CASSCF calculation was carried out prior to the MRCI calculation for PES_M. Also, PES_M applied a scaling method, mainly to better describe dissociation of the diatomics when approaching full dissociation of the system.¹⁹ Hence, the two PESs considered in the present work are based on different quantum chemical calculations and use different strategies to represent them.

The two representation schemes (RKHS and PIP) both yield highly accurate representations and allow for the computation of analytical forces required for running MD and QCT simulations. For PES_B, a grid of 30 (O + N₂) and 28 (N + NO) grid points for the radial coordinate between [1.4-12.0] a₀ was used with the diatomic distance covered by 20 and 21 points, respectively, between [1.55/1.50-4.0] a₀. Additionally, 13 Legendre quadrature angles were used to effectively sample the angular component, yielding 7800 and 7644 grid points for the (O + N₂) and (N + NO) channels. The ³A^′ surface for PES_M was constructed from a considerably smaller number of energy evaluations: 592 and 1706 grid points for the two channels at the SA-CASSCF/MRCI level of theory. For the two radial coordinates, the grid covered ranges [1.7–4.9] a₀ and [1.7–7.5] a₀ for and , respectively.

Quasi-Classical trajectory simulations were used to obtain state-to-state cross sections. For PES_B specifics were described previously^{18, 30} and hence only a concise summary is given here together with details of new simulations that were needed to increase the fidelity of the model. For the thermal rates, initial reactant states () were sampled from a Boltzmann distribution at 10000 K. The initial rovibrational quantum states for NO and N₂, using PES_B, were obtained from the semiclassical theory of bound states.³¹ The internal energy of the diatom, , is given by the condition

where

Using PES_B, Eq. 1 was solved numerically using the Newton-Raphson method. With an in-house QCT program the 6329 initial states for NO, and 8733 for N₂ were determined up to the corresponding dissociation limits.³² The and quantum numbers reach maximum values of NO(, and N₂(, for PES_B.³³ For both the N + NO and the O + N₂ reactions a total of trajectories were run for each initial condition. The impact parameter was sampled using stratified sampling and divided into six strata with a₀. Final rotational and vibrational states of the product diatom were assigned using a Gaussian binning scheme,³⁴ which results in faster convergence. The state-to-state cross section was then computed from , where is the ratio of trajectories that reacted for a particular initial condition.

For PES_M, all necessary STS information had been determined in previous work using the CoarseAIR^{35, 11} code which is an updated version of the original VVTC code developed at NASA Ames Research Center.³⁶ These simulations also employ stratified sampling over the impact parameter with binning of 0.25 a₀. The maximum simulation time was au with an initial time step of 5 au in the adaptive-time step Adams-Bashforth-Moulton integrator³⁷ for each trajectory. For the thermal rates, trajectories were simulated for the initially sampled reactant states out of the Boltzmann distribution. The maximum internal quantum numbers for PES_M are NO(, and N₂(, for PES_M.¹⁹ For the STS rates, QCT simulations were performed for both, the NO+N and N₂+O reactions, using 50000 trajectories for each initial condition. It should be noted that this number of trajectories does not suffice to converge all state-to-state cross sections.

For PES_B, a neural network-based model to predict the full set of STS cross sections was used. Evaluation and analysis of the existing STS2019 model³³ for the NO() + N(⁴S) (forward, downhill) reaction indicated that extensions with respect to the state-space covered were required. For the highest state of NO, the predicted cross sections did not match those from explicit QCT simulations, as discussed further below. As a consequence, STS2019 needed to be improved and retrained.

The grid of initial conditions (IC) on which QCT simulations were carried out for STS2019 consisted of , 12 rotational states , and 15 collision energies eV which amounts to 1800 different ICs. This was extended by running additional QCT simulations for the and identical and grids as before to improve the ST2019 model. For these 360 ICs, QCT simulations were carried out each to determine the necessary STS-information.

For the reverse (”uphill”), i.e. considerably lower probability reaction N) + O(³P) NO() + N(⁴S) a new STS model was also trained from reference QCT simulations using PES_B. For the QCT simulations, initial conditions were generated from a predefined grid: ; and eV due to the energy difference for the uphill process. A total of 160000 trajectories for each initial condition were run with impact parameters a₀. Overall, 20790481 state-to-state rates were determined for the reverse reaction.

For training the NN model, 12 diatomic properties were used as input features: internal energy, vibrational energy, vibrational quantum number, rotational energy, rotational quantum number, angular momentum of the diatom, relative translational energy, relative velocity, turning periods of the diatom, rotational barrier height, and vibrational time period of the diatom.³³ Together with the STS cross sections, this constituted the input to train the updated forward (STS2025) and reverse NN-based models required for the present work. From the energy dependent STS cross sections, for given , the respective kinetic STS reaction rate dictionary was built using the predicted NN cross-section according to , where and K .

In addition to the atom-exchange reactions, it was found that the dissociation channels needed to be included for completeness in the coarse grained modeling¹². Consequently, further QCT simulations for the two dissociation channels N) + O(³P) N(⁴S) + N(⁴S) + O(³P) and NO() + N(⁴S) N(⁴S) + N(⁴S) + O(³P) were carried out. For the N+NO and O+N₂ dissociation N + N + O reactions, 50000 trajectories for each initial condition (8733 and 6329 for N₂ and NO, respectively) were run using PES_B. For each initial condition, the collision energy was sampled from a Maxwell-Boltzmann distribution at K. The initial distance between the reactants was 20 a₀ with a maximum impact parameter a₀. The atoms were labeled as N_A, N_B, and O so that (in)elastic and atom exchange reactions can be distinguished. The number of dissociative trajectories was counted independently for each initial state.

For PES_M, previous work¹² determined a total of 12’742’744 STS kinetic rates for the NO() + N N₂() + O reaction based on QCT simulations. The same trajectories were also used to extract the relevant information for full dissociation to atomic products which constitutes the University of Illinois state-to-state model (STS-UI).

Master equation simulations were carried out for an isothermal heat bath condition to study the influence of the forward and reverse kinetics obtained from the two different PESs. For the forward and reverse heterogeneous exchange processes, the set of master equations reads

where is the number density of the species/level, denotes time, and and are multi-indices that refer to the NO and rovibrational states, respectively. is the net mass production rate due to the forward and reverse heterogeneous exchange reactions (superscript ), NO() + N(⁴S) N) + O(³P). denotes the state-to-state forward exchange reaction rate. The reverse reaction rate, , can be determined by employing either the micro-reversibility (MR) (i.e., where denotes the equilibrium constant) or through direct calculation (DR) by means of QCT simulations or evaluating the NN-trained models to yield . If microreversibility is used, was determined from the species’ partition functions.¹²

The master equations (Eqs. 3 to 6) were then numerically integrated using plato (PLAsmas in Thermodynamic nOnequilibrium) ²⁰, an object-oriented library for nonequilibrium plasmas developed within the Center for Hypersonics and Entry Systems Studies (CHESS) at the University of Illinois at Urbana-Champaign. In the present study, the state-to-state master equations in Eqs. 3 to 6 are solved to obtain time evolution of the microscopic population distributions, and .

The sums in Eqs. 3 and 4 over the molecular internal states define the macroscopic reaction rate coefficients, and , for the forward and reverse heterogeneous exchange reactions. In this case, the set of master equations governs the species’ macroscopic concentrations (i.e., , , , and ) and and are defined using the Arrhenius parameters in Table 1. In addition to the chemical rates, a model is required to describe the internal (i.e. rovibrational) energy relaxation of the molecular species, and NO. In the present study, the concept of the conventional two-temperature (2-T) model ³⁸ is adopted for simplicity, resulting in that the classical Landau-Teller formulation³⁹ was employed:

where and are the average vibrational energy in the equilibrium and the non-equilibrium states, respectively. denotes the average vibrational energy of species . The vibrational-translational (VT) relaxation times of the collision pairs considered in the present study, and , were taken from the work of Park ³⁸. It is noted that for the VT energy transfer the STS and Arrhenius treatments differ. In the Arrhenius treatment describes the bound-bound transition of the same species (e.g. inelastic + homogeneous exchange), for example, N+ON+O instead of product of NO()+N. Conversely, for the STS treatment, the rovibrational-translational (RVT) transfer occurs all simultaneously based on the STS rate coefficients, which might include the heterogeneous exchange, unlike the Arrhenius treatment. In the discussion section, a somewhat more balanced approach is discussed.

To handle detailed balance for multiple competing reactions the individual kinetic processes were balanced separately.⁴⁰ Hence, for the forward reaction N₂+ONO+N the reaction N₂+ONO+N is used for detailed balance, whereas for dissociation N₂+ON+N+O the reverse reaction N₂+ON+N+O is relevant for detailed balance. In other words, each chemical kinetic process has a corresponding reverse process, and the different reaction pathways remain independent of one another, even when detailed balance is enforced.

The master equation simulations were started from equal fractions of NO and N unless otherwise stated. For the initial gas pressure and internal temperature Pa and K were used, respectively, and the temperature of the surrounding heat bath was 10000 K.

The topography of the two ³A^′ PESs used in the present work, PES_B and PES_M, is reported in Figure 1. Both PESs are drawn at identical values of the isocontours with the zero of energy taken as the global minimum of each PES, respectively. The general shapes and critical points of features of the two PESs are similar to each other. These include the location of the global minimum at , the transition state at around , and the distant transition state at around . The N₂O PIP-PES is generally steeper around the global minimum and the two PESs differ in their anisotropy in the long-range part of the PES ( a₀). Still, given the different levels of electronic structure theory and their representation, the similarities between the PESs are rather striking.

For the forward (downhill, see Figure 2) reaction NO(X) + N(⁴S) N₂(X) + O(³P) on the ³A^′ PES two pathways were reported.^{19, 18} The first one amounts to N-attack at the oxygen-side of NO, compression of the NON angle, formation of N₂ and ejection of the oxygen atom, see Figure 2. An alternative pathway, not shown in Figure 2, starts at NO(X) + N(⁴S) and directly leads to minimum M1 over a single, lower-lying, transition state and involves a collision of the incoming nitrogen atom with the nitrogen-side of NO. TS4 is 32.9 kcal/mol above the entrance channel for PES_B which compares with 37.3 kcal/mol for PES_M. For the second pathway the single TS is 9.6 kcal/mol and 10.5 kcal/mol above the NO() + N(⁴S) entrance channel for PES_B and PES_M. PES_B has a 73.1 kcal/mol energy difference between NO + N and N₂ + O, whereas for PES_M this value is 75.8 kcal/mol.

Computing complete and converged dictionaries for state-to-state rates even for atom+diatom reactions is a daunting task. This motivated the development of machine learning-based approaches (statistical models) that require only a fraction of explicitly QCT-determined STS-rates and predict the remaining rates from a trained neural network.³³ Evaluating the trained model on the ³A^′ ground state of the NO() + N(⁴S) N₂(X) + O(³P) reaction for all accessible state NO pointed to some deficiencies of the STS2019 model.³³ The distribution of all state-to-state cross sections (Figure S1) using STS2019 (green trace in Figure 3) features a non-negligible fraction of unusually large cross sections. This prompted the explicit validation of such cross sections compared with results from explicit QCT simulations. For this, trajectories were run for for all and eV ( K) each initial condition at K. Upon analysis it was found that unreliable predictions invariably concerned high states with , i.e. near the NO dissociation threshold, see red circles in Figure 3A. Consequently, the STS2019 model needed to be extended with STS rates starting from high quantum number and subsequent retraining of the NN.

The additional QCT simulation required to improve and retrain the STS2019 model yielded an additional 1261 non-zero STS cross sections for the NO + N N₂ + O reaction which were added to the training set. After retraining the NN, model STS2025 was found to be considerably improved, see the black and blue circles in Figure 3B. Furthermore, the distribution of cross sections does not contain unusually large values, red trace in Figure S1. It is also relevant to note that the distribution of cross sections from STS2025 and those determined from QCT simulations using PES_M are largely consistent, see Figure S1.

The distribution of these STS cross sections (blue) and those from models STS2019 (green) and STS2025 (red) is reported in Figure S1. For STS2019 the fraction of large STS cross sections around a is clearly visible. On the other hand the range of the retrained STS2025 rates is closer to the STS cross sections from simulations using PES_M whereas the shapes of the blue and red distributions differ to some extent. This is understandable as PES_B and PES_M differ in their topography and the number of bound states they support and the number of QCT simulations that were used for STS2019/STS2025 and for STS-UI differ as well.

The performance of the new STS2025 model for the reverse reaction is shown in Figure S2. For a test set the model preforms satisfactorily for cross sections in the range [0-0.09] . Further improvement of the model, for example through hyperparameter optimization, is in principle possible but was not attempted in the context of the present study.

In the following, the chemical evolution of the NO() + N(⁴S) N) + O(³P) system will be evaluated using different models and approximations. PLATO determines the time evolution of the chemical composition of the system given initial populations of the species involved and either a thermal rate expression such as an Arrhenius fit or a dictionary of STS rates from a machine learned model (STS2025) or from explicit QCT simulations (STS-UI). The Arrhenius parameters (Arr_B and Arr_M) from QCT simulations using the two PESs (PES_B and PES_M) are reported in Table 1.

When running the PLATO simulations, the rates for the reverse (“uphill”) reaction N) + O(³P) NO() + N(⁴S) can either be determined from assuming detailed balance using the available forward rates for NO() + N(⁴S) N) + O(³P), or the reverse rates need to be evaluated explicitly. This can be accomplished from either using the STS2025 NN-based models (for PES_B) or from explicit QCT simulations (for PES_M).

First, the assumption of microreversibility for the reverse reaction is tested for PES_B. For this, PLATO simulations were carried out using Arr_B parameters for the forward reaction (subscript “f”) with reverse rates from assuming microreversibility (dashed lines in Figure 4A). Using Arrhenius rates , for the forward and reverse reaction (solid lines in Figure 4A) leads only to insignificant changes. Hence, the results reported in Figure 4A indicate that the final Nconcentration does not depend on whether reverse rates are determined from assuming microreversibility or an explicit Arrhenius rate is used. Also, the “ignition point” (time at which [) is s when assuming microreversibility compared with s using explicit rates for forward and reverse reactions. Results using Arr_M parameters from PES_M are reported in Figure S3 and confirm that PLATO simulations using explicit reverse rates and those assuming microreversibility yield essentially identical equilibrium concentrations and ignition points.

Next, the STS2025 model with and without assuming microreversibility for the reverse rate - dashed and solid lines in Figure 4B - was used to follow the species’ kinetics. The equilibrium species’ concentrations only differ by 0.03 units whether or not explicit reverse rates were employed. Also, the ignition point for both types of simulations are at s. This is an increase by s compared with the simulations based on Arrhenius-rates. The equilibrium product concentrations between Arrhenius- and STS2025-based simulations differ, however, by about 20 %.

In summary, the species’ kinetics and ignition times do not differ whether Arrhenius- or STS-based rates are employed but the equilibrium amount of product generated differs by %. This point will be further discussed below.

The two PESs - PES_B and PES_M - considered in the present work support different numbers of bound rovibrational states. This is due to the fact that a) the two PESs were determined at different levels of theory (MRCI+Q/aug-cc-pVTZ for PES_B and MRCI/maug-cc-pVTZ with dynamical scaling (DSEC) for PES_M, respectively) and b) the reference energies were represented as a RKHS and PIPs, respectively.^{18, 19} Therefore, the influence of the PES used in the QCT-simulations, i.e. PES_B vs. PES_M, was probed next.

Figure 5A compares the concentration profiles from Arrhenius-rates for the forward reaction with reverse rates based on assuming MR. Using PES_B the Arrhenius-rates Arr_B yield the same limiting concentration for (solid line) as does Arr_M determined from running QCT simulations on PES_M (dashed line). The ignition point from Arr_M occurs at s compared with s from simulations using Arr_B. When using the STS information from QCT simulations on the two PESs (rates from STS2025 or STS-UI, respectively) the limiting concentrations converge to , see Figure 5B. Again, and consistent with using Arrhenius rates, the ignition point from STS-UI is about 1 order of magnitude earlier than that from using STS2025.

Finally, Figure 6 compares the kinetics from using MR-reverse rates (panel A) with those using explicit STS-reverse rates (panel B) based on simulations with STS2025 and STS-UI. The data in Figures 6A and 5B is identical but reported again for direct comparison. Figure 6B shows that with explicit rates for the reverse reaction the kinetics for both, STS2025 and STS-UI, yields asymptotically on the s time scale. This compares with if reverse rates are obtained from assuming MR, see Figure 6A. The ignition points when using STS2025 and STS-UI also remain unchanged whether reverse rates are from assuming microreversibility (panel A) or explicitly given from the STS dictionaries (panel B).

In summary, the species’ kinetics is surprisingly insensitive to whether input data for the coarse grained simulations was obtained from PES_B or PES_M, except for a slightly different ignition point. This pertains to both, whether or not microreversibility for the reverse rates was assumed. Also, the results from using Arrhenius-parameters instead of STS-information are consistent and provide additional reassurance as to the robustness of the findings independent on which of the two PESs was used.

Up to this point the results from PLATO simulations using STS-data generated from PES_B and PES_M are surprisingly similar, see Figures 5 and 6. However, modeling of reactive flow without accounting for the possibility of the reaction products to fully decompose into atomic fragments at temperatures K or above is not entirely realistic. This is reflected, for example, in the equilibrium product concentrations, see Figures 4 to 6: in all cases the diatomic species’ concentrations converge to finite values at the end of the simulation. Given that all simulations were carried out at 10000 K and in an isothermal heat bath, all molecules should dissociate in thermodynamic equilibrium. This, however, is not possible without including the dissociation channel. Therefore, PLATO simulations were also carried out by including fragmentation N) + O(³P) N(⁴S) + N(⁴S) + O(³P) and NO() + N(⁴S) N(⁴S) + N(⁴S) + O(³P), see Figure 2. For PES_B the required initial state-dependent dissociation rates were determined from explicit QCT simulations. For PES_M the same procedure was followed.

Figure 7 reports the species concentrations within isothermal heat bath conditions (see Methods) as a function of time over 15 orders of magnitude, until all reagents’ concentrations were stationary, which occurs after s. Simulations using information from PES_B are reported as solid lines and those using PES_M as dashed lines. The overall behaviour of both models is comparable over the entire time range. However, the ignition point for PLATO simulations from PES_B is later by about half an order or magnitude compared with PES_M. In other words, the generation of atomic oxygen (black) is faster with data from PES_M. The equilibrium compositions reached (1/3 for [O] and 2/3 for [N], as required) are virtually identical and the temporal evolution of [N₂], [NO], [N], and [O] match rather closely given the different origins of the STS-rates, the dissociation rates and differences in constructing and representing the PESs underlying the QCT simulations.