Towards a hybrid semi-quasi-classical trajectory method for the accurate calculation of vibrationally inelastic probabilities in atom–diatom collisions

Abstract

The quasiclassical trajectory method (QCT) for molecular dynamics is routinely used for the calculation of complete sets of cross sections and rate coefficients in vibrationally inelastic collision processes concerning air molecular species. Such processes are key for kinetic models of interest in accurate simulations of practical relevance, such as green fertilizer production by non-equilibrium (or cold) plasmas. This popularity is due to the general good compromise of QCT between accuracy and computational effort required. However, QCT at low total energy often fails to accurately describe low-probability vibrationally inelastic processes. The major impact of QCT data inaccuracy in cold plasma modeling is caused not by poor statistics associated with low probability events nor barrier tunneling, but by the inadequacy of QCT to capture some vibrational distribution features in the final products. The recently introduced collision remoteness concept allows a clear-cut distinction between purely nonreactive and quasireactive regimes. The first, where transition probabilities are often small, is easily treated accurately by semiclassical (SC) methods. The quasireactive regime is well modelled by QCT. This suggests a possible merging of QCT with a suitable SC method, with a reasonable capture of the most accurate partial result of each method. A first proposal of a merging procedure is presented here, with an excellent level of agreement with both accurate quantum-mechanical (QM) time-independent and time-dependent close-coupling calculations for vibrationally inelastic O + N₂ collisions. This level of agreement is far better than the one found using the QM coupled states approximation in the same conditions. Accurate QM treatment of air species with the level of vibrational detail required by kinetic simulations, including air species, is prohibitively far from computational feasibility, so development of accurate approximations remains crucial. The details of the hybrid method and its possible impact on modelling of non-equilibrium plasmas of technological interest are discussed.

---

1. Introduction

Vibrationally inelastic molecular collision processes are at the core of studies concerning “cold” non-equilibrium plasmas.^1,2 Cold plasmas (CP in the following) are of interest for problems of large impact like CO₂ dissociation for utilization in e-fuels³ and green production of fertilizers from waste recycling.⁴ Under nonequilibrium conditions, descriptions of species concentrations and reaction pathways must use kinetic models since neither thermal nor chemical equilibrium relations are trustworthy. In cold air plasmas electron–molecule, molecule–molecule and molecule–surface processes all play important roles. It is often necessary to treat different internal states of the atomic or molecular plasma components separately, both because rate coefficients may depend on internal state and because internal excitations may be substantial energy reservoirs in the technological system.⁵ Models of cold plasmas intended to assist process improvement efforts⁶ should therefore include large libraries of microscopic dynamical data.

It is impractical to populate such a library entirely with experimentally determined rate coefficients or cross sections, even for just the neutral atomic and molecular species present in a cold air plasma. The numbers of species and states are too large and the relevant range of collision energies too wide. Empirical estimates or theoretical calculations are necessary.

Theoretical calculations of reactive or inelastic rates require potential energy surfaces describing the interactions between colliders and dynamical methods to model the resulting motions. Potential surfaces of adequate accuracy are already available for many collider pairs important in air plasmas, including many of the small molecules containing H, C, N and O atoms. The principal barrier to practical computation of state-to-state kinetic data is the effort required in the dynamics calculations. “Exact” quantum mechanical methods, in the sense that no physical approximations are used beyond the Born–Oppenheimer approximation and truncated basis sets, are available for both inelastic and reactive collisions of small molecules. Their computational cost increases rapidly as molecular energy level spacings decrease and as higher rotational states are included in a calculation. These exact methods are not currently practical for generating full kinetic data sets for plasma modeling. Like experimental state-to-state measurements, they can be used to provide reliable validation benchmarks against which more approximate but more practical methods may be tested.

Several types of more practical dynamical methods are available. Approximate quantum methods are often built from the exact methods by neglect of specific terms in the collision Hamiltonian that make the calculations more complicated, but in some circumstances may not affect the results very much. One of the best-known such approximations is the coupled states method.^7,8 Semiclassical methods combine some kind of wavefunction picture with a classical description of the overall motion. Early semiclassical methods often used simple models of the forces and trajectories,^9–11 but more recent ones often use detailed classical trajectories to provide “backbone” dynamics.^12–14 Quasiclassical methods use classical mechanics to model each collision, but select their initial conditions and bin their final results to reflect quantization.¹⁵

Temperatures involved in cold plasma processing are never lower than 300 K, and often significantly higher.⁵ Approximate quantum methods can still be expensive for air species at those temperatures. The computationally efficient semiclassical and quasiclassical approaches seem promising for CP applications, since the quantum effects they neglect decrease in importance as atomic masses and collision energies increase.

Computational effort can also be reduced using suitable neural networks on sparse and even heterogeneous data from different sources to obtain a sort of high level interpolation.^16–18 The approach is being exploited using QCT, SC and QM methods for calculating a small subset of state-to-state data, then constructing the whole set of transitions using artificial intelligence (AI) techniques, exploiting the redundancy of the data. The approach is appealing for drastically reducing the computational time typically associated with the compilations of large databases. However, vibrational energy exchange in molecular collisions is not necessarily treated with sufficient accuracy by SC and QCT methods. The subset used as input for artificial intelligence reconstruction can be rough, giving origin to a poor final output. Whether AI reconstructions are used or not, it is important to study the accuracy of dynamical simulations.

One of the authors (F.E.) has started a systematic investigation activity in several collaborations in the past years,^19–23 comparing QCT, accurate QM calculations and experiments of light and heavy collisional systems, with the purpose of clarifying the limits of application of QCT, in particular for vibrational energy transfer, and possibly integrating it with different methods. While a light system like reactive H + HeH⁺ has been simulated with high accuracy by QCT,²⁰ even considering detailed vibration, the vibrational energy transfer of heavy O + N₂ collisions computed by QCT shows important defects at energies of relevance for technological simulations.²¹ O + N₂ enters all air plasmas. Furthermore, many other combinations of oxygen and nitrogen species in air plasmas, as well as carbon and oxygen in CO₂ treatments, are probably affected by the same issue, even if it is difficult to assess the problem because accurate QM calculations are difficult for those systems. SC methods can partially solve this problem in specific conditions. Studying QCT and SC methods and their possible combinations may lead to promising advances in provision of detailed and accurate sets of molecular dynamics data.

The collision remoteness (CR) concept^21,24–29 can be useful for distinguishing different regimes in which different trajectory methods (QCT and SC) can be used. This simple computational tool is the basis for categorizing vibrationally inelastic collisions as purely nonreactive (PNR) or quasireactive (QR), in a sense that will be introduced in the following sections. A general description and some details of all the “standard” methods used in this work is in the next Section 2. In Section 3 quasireaction is defined and discussed, while in Section 4 a description of trajectory behavior around the quasireactive limit is presented. In Section 5 a hybrid semi-quasi-classical method is described, with comparisons of semiclassical and quasiclassical probabilities with accurate QM results, obtained in this work with the time-independent close coupled (TICC) and time-dependent or wavepacket close-coupled (TDCC) methods. The comparisons identify strengths and weaknesses of the various computational methods for the vibrational energy transfer problem in molecular collisions. In Section 6, full cross sections calculated with both SC and QCT methods, and a possible combination of the two, are compared with one another and with TICC calculations where possible. The framework that can be envisaged by these first comparisons seems to indicate that suitably mixing different methods could be extremely useful for getting accurate data in total energy ranges that at the present time are unattainable by any single method and still unexplored.

2. Details of molecular dynamics methods used in this work

In this work five different standard methods for molecular dynamics calculations have been used. All the methods have been applied to the same problem, the calculation of probability of translational-to-vibrational energy exchange in O + N₂ (v = 0, j = 0) → O + N₂ (v′ = 1, j′) and a few other initial or final states on the ground state potential energy surface A″ from ref. 30. The O + N₂ collisional system has a high threshold to reaction (3.26 eV), and this feature has important consequences on the inelastic channel, as will be clear in the following.

Vibronic transitions towards other electronic states can affect low energy O + N₂ collisions, as shown in ref. 31 and 32. These aspects, which do not affect the critical collision energy region around the reaction threshold of special interest here, are outside the scope of the present work.

In the following, a description of each dynamical method is presented.

2.1. QCT

The first method to consider is the well-known QCT.¹⁵ Our implementation was initially described in ref. 33 and then improved over the years with grid-distributed, parallel and multithreaded versions. The current version includes the possibility of calculating J-dependent probabilities (J the total angular momentum quantum number, not to be confused with j the diatomic rotational angular momentum) as described in ref. 29 and 34 and the value of collision remoteness associated with each trajectory,^25,27,29 as detailed in the following section. A good compromise between accuracy and computational efficiency is obtained with step-by-step backintegration of trajectories with halved time step and comparison with the forward calculation, with typical tolerance of 10⁻⁹ Å on space coordinates and 10⁻⁹ Å fs⁻¹ on velocities (see for example ref. 35). At each time step, if the test fails, the calculation is repeated with a halved time step, up to a maximum of 6 times. This procedure guarantees a computational high speed when interaction is low, typically in asymptotic regions, and very good accuracy wherever interaction is strong. The maximum interaction distance is fixed at 15 Å. The number of trajectories calculated (over 96 millions) allows a statistical noise at least ten times lower than any non-zero probability value, often much better. Initial and final rotations are calculated in a standard way,¹⁵ using the relation j_cl² = j(j + 1)ħ² between the classical rotational angular momentum j_cl and the j quantum number. Initial and final vibrations are calculated as described in ref. 15, using the WKB approximation for final vibrational actions. Standard histogram binning is used (the final quantum number is assigned by the nearest integer procedure). It is worth noting that Gaussian³⁶ or any other binning scheme in this work are not of help, because the problem, as will be clear in the following sections, is the treatment of extremely narrow final distributions, for which binning cannot be used. Conservation of rotational angular momentum parity in inelastic processes involving N₂ is enforced by assigning a weight zero to forbidden and 2 to allowed transitions, respectively.

Probabilities are calculated as in ref. 34, i.e.:

where l is the orbital angular momentum, with P(J,l) = N_r/N_tot, N_tot being the total number of trajectories for given J,l values and N_r the number of trajectories satisfying a given criterion (for example reaction, quasireaction, inelastic events, etc.).

The cross section was obtained by summing probabilities over J values according to

with k² = 2μE_coll/ħ², interpolating linearly among J values (one in four J values are calculated by dynamics). The [2 min(J,j) + 1] factor produces a conventional cross section when used with the probability definition in eqn (1).

Probabilities and cross section calculations described in these last two equations extend also to the two other trajectory-based methods described below (tFHO and CMM).

2.2. tFHO

tFHO stands for trajectory-forced-harmonic-oscillator. This acronym describes a trajectory-based molecular collision method, based on the DECENT-INDECENT model in ref. 13 and with code written from scratch by one of the authors (F.E.). The diatomic target is a harmonic oscillator, modeled on the diatom of interest, excited by the atomic projectile. The small “t” in the acronym stresses the difference with the FHO method,^11,29 which is based on the same forced harmonic oscillator model, but uses purely analytic dynamics on a simplistic model interaction potential. In tFHO, the trajectory is computed using the same propagator on the same PES as in QCT. The initial conditions are also selected like those used in QCT except that the initial diatom interatomic distance is always set to the equilibrium bond length, and the vibrational momentum to zero, instead of the distance and momentum corresponding to a randomly determined vibrational phase. There is therefore no dependence on initial vibration in the tFHO trajectory dynamics. The trajectory analysis consists of obtaining the diatomic excitation by measuring the final vibrational energy or, more accurately, the final vibrational action by a WKB method as in ref. 37. Much higher accuracy (5 significant digits) is required in this vibrational analysis than in QCT, in order to get the level of agreement with quantum-mechanical results shown in the results sections. Using final vibrational actions (instead of energies) not only implies much higher accuracy, but solves the problems of vibrational energy definition and the relation between the reactant diatomic potential and the harmonic potential of the tFHO model. This is because the actions are calculated using only the “true” diatomic potential, not the harmonic one, requiring also a negligible computational time with respect to trajectory computation. The state-to-state vibrational energy exchange probability between initial and final quantum states v and v′ is then computed in the forced harmonic oscillator approximation³⁸ with

P = v!v′!e^−εε^v+v′S_vv′², (3)

wherewith k_max = min(v,v′). The final vibrational motion from the trajectory is characterized by

ε = ΔE_ν/ħω (5)

or, as in ref. 37 and in this work,

ε = Δν_cl, (6)

where ΔE_ν and noninteger, classical Δν_cl are the changes (the final values, indeed) of vibrational energy and vibrational action, respectively.

It is worth noting that the initial vibrational quantum number v is introduced only at the analysis stage through eqn (3), while the initial classical vibration in the trajectory is exactly zero. Detailed balance would not be satisfied if the trajectories were computed at the collision energy of interest, since the change in translational energy during the trajectory is usually different (often much smaller) than the energy difference between the initial and final diatomic quantum states. We therefore adopt a simple energy shift technique very similar to the one adopted by several semiclassical approaches^39,40 to enforce detailed balance. We associate the transition probability obtained with eqn (3) for a particular pair of initial and final quantum states from a trajectory computed with kinetic energy E_kcl with a physical collision energy (see appendix for details):

E_coll(tFHO) = E_kcl + ΔE/2 (7)

where ΔE = E(v′,j′) − E(v,j) (it is worth noting that initial and final roto-vibrational energies are used, differently from ref. 38). Eqn (7) implies that a set of trajectories should be calculated for each pair of initial and final states for each collision energy value required. In this work, probabilities and cross sections were calculated with tFHO for a single set of E_kcl values, giving different collision energy lists E_coll(tFHO) for each pair of initial and final states. Then, the probabilities or cross sections for each pair were interpolated on collision energy to obtain values corresponding to physical collision energies used in the corresponding QCT calculation. tFHO calculations are performed considering full dimensional collisions, in a “quasiclassical fashion”, i.e. by randomly and uniformly varying all the other quantities not related to molecular vibration: initial atom-molecule relative position, rotational phase and axis (for given values of total angular momentum and collision energy). Probability results are then obtained by Monte Carlo averaging over all these quantities.

tFHO, as used in this work, is computationally cheaper than QCT, because with QCT one has to compute a new trajectory set, including averaging over vibrational phase, for each initial state. In addition, with tFHO even very low inelastic probabilities, 10⁻¹⁰ or lower, can be calculated easily, in contrast with QCT. The physical model behind tFHO is described for example in ref. 38, in which a harmonic potential moves according to the classical evolution of the center of mass of the diatomic oscillator during the interaction with the atomic projectile. This interpretation permits an analytical solution, as explained in ref. 38 and 41. From this analytical solution of the motion, as in eqn (5) of ref. 41, it is easy to recognise that the forced oscillation induced by the projectile passage near the diatom generates a rigid body oscillation of the initial vibrational wavepacket (which is determined as the harmonic vibrational state v of the reactant diatom). This wavepacket rigid body oscillation, equivalent to the well-known coherent state as in ref. 42, can be shown to strictly follow the oscillation of the classical diatom subject to the atomic projectile interaction. This is possible because of the special feature of the harmonic potential, in which the wavepacket center of mass exactly oscillates in a classical fashion.^38,41 This ideal interaction will only be realistic if two conditions are fulfilled: the diatomic potential is harmonic, and the force f between the projectile A and the diatom BC does not depend (significantly) on the interatomic distance BC, but only on time f = f(t) through the evolution of the other classical coordinates.^38,41 In real cases, the first condition is easily satisfied considering sufficiently low-lying states. The second one is more subtle and generally overlooked in the literature, but important. A weak dependence of the A–BC interaction on BC distance can be well approximated when the atomic species involved are equal and the two distances AB and AC remain large in comparison with the BC distance during the collision (which is another way of saying that the BC interaction prevails on AB and AC interactions during the whole collision). In order to satisfy this condition, the distance between A and the BC center of mass must remain significantly greater than BC during the entire collision. This immediately implies that no rearrangement can take place even momentarily during the collision. This will be the basis for the “collision remoteness” concept. If the atomic species involved are different, a suitable distance scaling is needed and will be described in the following section.

At the end of the interaction, the tFHO inelastic final analysis corresponds to model wavepacket final analysis, i.e., a projection of the wavepacket onto the harmonic states of the initial diatom. This projection is obtained analytically, assuming that the final wavepacket is the initial one plus a rigid body oscillation acquired during the dynamics. It is exactly this rigid body oscillation that is modeled by classical trajectories. This model remains valid as long as the rigid body hypothesis is well approximated, i.e. when the A–BC interaction is sufficiently weak. This explains why low-energy vibrational energy transfer probabilities are very well approximated by tFHO, as will be clear in the following sections.

Because tFHO is based on a harmonic approximation of the diatom potential, only low-lying vibrational states and non-reactive events with relatively weak interactions can be treated accurately. Inelastic probabilities arising from strong, reactive-like interactions can be computed by this method, but are likely to be inaccurate, because this kind of collision is beyond the scope of the tFHO physical model. In a reactive-like interaction during a vibrationally inelastic collision, the diatom oscillator is substantially distorted or completely destroyed and then re-formed at the end of the collision, so there is no hope that it can be accurately modeled with a fixed harmonic potential throughout the dynamics. Some possibilities of extending tFHO validity to reactive-like inelastic processes could be studied with the use of a time-dependent harmonic potential, as in ref. 43.

2.3. CMM

The classical moment method^44–46 is also used in this work in a very simple form as a comparison of another classical trajectory method more closely related to QCT method than tFHO. Strictly speaking it is QCT, with a final analysis based on suitably manipulating final moments of the vibrational distribution to calculate final probabilities, instead of using the well-known binning. For pratical reasons, it will be briefly called CMM. It has the great advantage that the same trajectories can be analysed as QCT or as CMM. In this work we extract only the first moment and use it to get the probability of the v = 0 → v′ = 1 transition. We estimate this probability as the difference between the final and initial classical vibrational actions normalized to Planck constant (exactly the ε of eqn (6) for tFHO, but for the initial vibrational action, which in this case depends on initial vibrational state), determined to 5 significant digits with the WKB method, as in tFHO. A Monte Carlo average is then applied to probability results, as in tFHO. No energy shift is required. The use of CMM, at least in this primitive form, is limited to relatively weak interactions, excluding reactive-like collisions.

2.4. TICC

The time-independent scattering calculations were carried out with version 14 of the Molscat program^47–49 or its PMP Molscat parallel variant. The close-coupled method was used with incorporation of centrifugal distortion in the vibrational wavefunctions (ITYPE = 7 in Molscat).

For cross section calculations, N₂ vibrational wavefunctions were computed from the accurate N₂ potential curve of Le Roy et al.,⁵⁰ evaluated with code from the LEVEL package.⁵¹ For transition probabilities, the asymptotic N₂ curve incorporated into the O + N₂ potential was used for consistency with the time-dependent wavepacket calculations. In both cases N₂ vibrational wavefunctions ϕ_v,j(r) were determined from the potential curve using the Fourier grid Hamiltonian method of Marston and Balint-Kurti.⁵² The grid was adjusted for different rovibrational basis sets; the most extensive diatomic calculation used 1600 grid points spanning the range 0.75–1.9 Å in N–N distance. Vibrational matrix elements of the interaction potential 〈ϕ_v,j(r)|V(R,r,γ)|ϕ_v′,j′(r)〉 were computed with Gauss–Hermite quadratures of up to 36 points, using vibrational wavefunction values interpolated from the grid with piecewise cubic Hermite interpolation routines from the SLATEC package. These matrix elements were then expanded in even Legendre polynomials up to λ = 20 with Gauss–Legendre quadrature using 32 points to provide radial strength functions V_vjv′j′λ(R) for the scattering calculations.

The coupled equations were solved with the hybrid log-derivative Airy propagator of Manolopoulos and Alexander.⁵³ The STEPS parameter was chosen to give 25 radial steps per half-wavelength at the highest scattering energy in each run. Propagations extended to O–N₂ distances of at least 15 Å. Partial wave sums were usually carried out by propagating equations for only every fourth value of total angular momentum J and multiplying the resulting cross sections by 4. The maximum total angular momentum was selected to ensure convergence of vibrationally inelastic cross sections. The values used in each energy range are shown in Table 1. These maximum J are not adequate to converge purely rotationally inelastic cross sections. Collisions cannot change N₂ rotational states from odd to even or vice versa. All time-independent calculations included only basis states with even j.

Table 1 Basis sets used in TICC calculations for different ranges of total energy (collision energy plus internal energy of diatomic initial state). The last basis listed was used only in J = 0 transition probability calculations

Energy range, cm^−1	vmax	Nlevel	Description	Maximum J	Maximum channel count
2500–6500	4	132	72/62/54/42/24	100	3830
6500–10 500	5	152	76/68/58/48/34/8	100	4620
10 500–14 500	6	219	86/80/72/64/54/42/26	148	7539
14 500–18 500	8	235	78/78/72/64/54/42/26/20/18	180	7343
4000–28 000	13	564	120/116/110/104/100/94/88/82/74/66/58/46/12/8	0	564

---

Basis set selection required particular care. These “inelastic-only” quantum calculations represent the motion of the system as a superposition of functions that are products of isolated N₂ rovibrational states with functions representing the relative motion of O and N₂. The collection of diatomic rovibrational functions must therefore be flexible enough to describe any important motion of the two N atoms during the collision. Calculations were usually done in groups of several energies for each basis set, with common energies at the overlaps between groups for convergence evaluation. Basis sets were constructed by beginning with all asymptotically open channels at some energy slightly above the highest energy intended for the group, adding or deleting states on the basis of initial convergence calculations done at J = 0, and then checking convergence with intergroup comparisons of completed cross sections. The resulting basis sets contained the largest number of rotational states at v = 0, and fewer rotational states at successively higher v. Table 1 shows the sets used in the study. In the table, v_max is the highest N₂ vibrational quantum number used in the basis and N_level is the number of N₂ rovibrational states it includes. The description shows the highest rotational level j associated with each vibrational level, beginning at v = 0. The channel basis is constructed by combining each rovibrational level with all values of the orbital angular momentum l consistent with |J − j| ≤ l ≤ J + j. The resulting set of coupled differential equations in the atom–diatom distance R breaks into two independent sets on the basis of parity; the maximum channel count gives the dimension of the largest such set that was solved in each energy range during computation of cross sections.

The 564-level set was used only for J = 0 transition probability calculations.

Inelastic cross sections were computed in the standard way⁴⁸ using S. Green's Molscat S-matrix postprocessing program sig_save.f. The program was modified to compute transition probabilities using eqn (1) of ref. 54,

where j is the initial rotational quantum number, g and g′ label initial and final internal states of the transition, l and l′ are channel orbital angular momenta, and is a T matrix element.

We estimate that all our reported TICC cross sections are accurate to 3%; most are better. Similar accuracy was obtained in the transition probability calculations up to just below 3 eV total energy, where basis set convergence becomes much more difficult to obtain.

2.5. TDCC

A quantum time-dependent wave packet coupled-channel method was used to compute state-to-state probabilities. The simulations were performed using the MAD-WAVE3 code, described in previous works.^55,56

The initial wave packet for a given reactant state defined by quantum numbers α ≡ J, M, ε, v, j, Ω₀ is expressed in reactant Jacobi coordinates (r,R,γ) as:

where ε, the parity of the total wave function under inversion, can be +1 or −1; is a symmetry-adapted Wigner function with total angular momentum J with projections M in the space-fixed frame and Ω₀ in the body-fixed frame; X_vj(r) are numerically calculated ro-vibrational eigenfunctions of the N₂(v,j) reactant; Y_jΩ0(γ) are associated Legendre polynomials; and is a Gaussian function centered at large R, representing the translational motion.

Wave packet propagation was performed using a modified Chebyshev propagator,^57–63 with approximately 9 × 10⁴ iterations for J = 0, ensuring convergence even at collision energies around threshold. The number of iterations decreases for higher J due to increased centrifugal barriers. Special absorbing potentials⁶⁴ at the radial grid boundaries were used to prevent unphysical reflections. Table 2 summaries numerical parameters used for the calculations. Although the transitions examined in this work are nonreactive and inelastic, the dynamics still accesses configuration-space regions typically associated with reactive processes due to the quasireactive features that arise in the narrow 2.7–3.1 eV interval. For this reason, the TDCC configuration space adopted here cannot be safely reduced: even modest truncations of the radial or angular grids lead to appreciable errors in the energy range where quasireactive dynamics and PNR interference dominate. The present grid choice is therefore essential to ensure a quantitatively reliable benchmark for the hybrid methodology.

Table 2 Parameters used in the TDCC calculations (all distances are given in Ångstroms)

Product scattering coordinate range	Rmin = 0.01, Rmax = 23
Number of grid points in R	960
Diatomic coordinate range	rmin = 0.01, rmax = 12
Number of grid points in r	420
Number of angular basis functions	200
Center of initial wave packet	R0 = 14.0
Initial translational kinetic energy (eV)	Ec = 1.5
Position of the analysis line	R∞ = 4.0
Number of Chebyshev iterations	90 000

---

Exact calculations for J > 0 would require summation over all helicities Ω′, which is computationally expensive. Therefore, the basis was truncated as where corresponds to the maximum allowed by the available memory on our supercomputer for J = 80.

3. Collision remoteness

The “collision remoteness” metric has been introduced in some recent works^21,24–29 as a tool for studying vibrational energy transfer of molecular collisions in trajectory-based methods, for potentially reactive systems (like practically all those containing air species). The name “collision remoteness” (CR) has been chosen because it is a measure of how far the projectile remains from the target during the whole collision (the higher the CR, the lower is the interaction between the projectile and the target, and vice versa). The basic concept was probably introduced for the first time by Ian W.M. Smith in the 1970s.^65–68 It consists of monitoring the internuclear distances along each collisional event of the atom A with molecule BC, detecting whether the BC distance ever becomes longer than one of the other distances AB and AC. This definition is appropriate in the presence of a symmetric system, like O + O₂, in which any possible product molecular species are the same as the reactant one. For nonsymmetric systems like O + N₂, Smith suggested⁶⁵ that the XY internuclear distance can be replaced by its counterpart that has been normalized with respect to its equilibrium distance XY_eq (i.e. ). For conciseness, it is useful to consistently use the normalized version, which is equivalent to the original for symmetric systems. This is the basis on which collision remoteness (CR) is built, by defining:²⁵where means minimum along the whole time evolution of each trajectory, while the internal minimum is between the shorter of the two distances, step-by-step. If the normalized distances between the incoming atom and the diatomic remain large compared to the diatomic bond length throughout a trajectory, is large compared to 1, and this corresponds to the purely non-reactive (PNR) condition.

Previous work^{21,24–29,65} has suggested that CR can be a useful metric for characterizing vibrational dynamics in inelastic collisions. The physical meaning of this metric is simple for symmetric systems, because in this case the condition coincides with crossing of the dividing surface to reaction during the encounter. In such a “quasireactive” (QR) encounter the trajectory necessarily enters a region where the minimum energy path is curved, so the vibrational degrees of freedom become coupled. The BC bond therefore suffers a transient weakening, or also breaking and successive re-formation.

Smith found for nonsymmetric Br + HCl collisions⁶⁵ that if somewhere along the trajectory, the distance becomes larger than or , the subsequent dynamics is distinctly different from the case in which it remains smaller throughout the whole trajectory. In particular, the final vibrational distribution is significantly wider and smoother in the first case than in the second one. He stated: “crossing of this surface [the one defined by ] is, of course, necessary for reaction, but this criterion also served to separate the non-reactive trajectories into those where little energy was transferred, and those where […] a wide distribution of final vibrational energies in BC was found”. We demonstrate in the next sections that the criterion also separates our nonsymmetric O + N₂ trajectories into QR and PNR classes in a very accurate way.

In the QR regime the collisional system tends to lose memory of the initial internal molecular state and give a final vibrational distribution that is relatively flat and smooth. In contrast, PNR collisions produce a final distribution that is highly peaked around the initial state, because there is only a slight classical modification of state of vibrational motion with respect to the initial condition.

Final vibrational distributions of reactive and inelastic rate coefficients from a given initial vibrational state for the symmetric systems H + H₂, O + O₂ and N + N₂ have been shown²⁴ to be almost coincident, with the exclusion of a few final states in the numerical neighborhood of the initial state. This observation can be explained by attributing the origin of the distribution to reactivity and quasireactivity respectively for reactive and inelastic rates, plus the contribution of the PNR highly peaked distribution around the initial state for the inelastic part.

Quasireaction, also called “frustrated reaction” or “recrossing”, also appears in theoretical and experimental studies of the differential cross sections for vibrationally inelastic scattering.^69–73 Vibrationally inelastic collisions are not always backscattered, as in the traditional model of compression of the initial diatom bond, but are also forward scattered, in a way that can be interpreted as a missed (“frustrated”) reactive process.

CR can easily be plotted for a set of trajectories as a function of final classical vibrational “quantum” number, the continuous quasiclassical analog of the quantum number. (In the following, it will be called quantum number for conciseness). Such plots display the differences in final distributions corresponding to different CR values.

Fig. 1 shows the two-dimensional distribution of collision remoteness and final vibrational quantum number for a sample of trajectories for O + N₂ (v = 0, j = 0) collisions at fixed total angular momentum value J = 0, for several collision energy values. The color map on the right of each panel shows the density of points, as defined in SI. The horizontal gray line in each panel of Fig. 1 is the quasireactive limit. It is exactly 1 for normalized CR, but in this work it is more appropriate to show CR as obtained without normalized distances, for reasons clarified below. CR is shown considering that for the case of O + N₂ collisions, for which AB_eq = AC_eq, eqn (1) can be rewritten

withconsidering N_2eq = 1.09 Å, NO_eq = 1.15 Å. Therefore, when , , so the gray line representing the QR limit is plotted at that vertical position.

Fig. 1 Distribution of unnormalized CR and vibrational quantum number in a small set of trajectories for O + N₂ (v = 0, j = 0) collisions at total angular momentum J = 0 and collision energies specified in the panels. Each point represents one trajectory. The color scale is proportional to trajectory density as detailed in the text. The horizontal gray line is the QR limit. By definition, trajectories are QR under this line, and PNR above it. Two completely different distributions are clearly present at the higher energies. The QR distribution under the gray line shows wide regions of low density and becomes wider with increasing collision energy. The PNR distribution, above the line, remains compact across the energy range and does not exceed the initial bin limits of v_initial ± 0.5.

The focus here is on the vibrational distributions over and under the QR line. QR trajectories, whose collision remoteness is lower than this limit by definition, are distributed over a relatively wide final range of values, definitely covering more than one final vibrational bin as collision energy is increased. PNR trajectories, which by definition are over the QR line, on the contrary, show a distribution that remains compact and narrow even at high collision energy. CR may therefore be useful in separating trajectories into groups that can be analyzed by different approaches appropriate for different final vibrational distributions. If the final distribution is peaked and narrow, as in the PNR case, with a typical width of the order of only one bin, the QCT result can be inaccurate²⁵ or even erroneously zero for classically forbidden transitions.⁷⁴ On the contrary, QR trajectories show a distribution that tends to be wide in comparison with bin width, so in this case QCT binning is able to correctly capture the result. It is worth noting that the QCT description of the QR distribution can be improved by using more trajectories, while its description of the narrow PNR distribution, even if apparently converging with few trajectories, shows poor accuracy. The wide final vibration distribution is a necessary condition for accuracy of QCT binning. This means that in a QCT calculation of a cross section for a vibrationally inelastic process, generally including both PNR and QR events, the PNR part cannot be accurate at low energy, whatever binning is used. This is because at sufficiently low energy, only the initial bin (the one corresponding to an elastic transition) collects some trajectories, the others are practically empty, independently of the binning scheme adopted. From what observed in Fig. 1, PNR tends to be predominant at low total energy, a range critically important in CP modellisation. Indeed, other quantum effects, for example reaction barrier tunneling, can make the QCT result inaccurate,^25,26 but this is a quite limited problem in the framework of air species collisions, because the species involved are quite heavy and the gas temperatures considered in typical applications of air plasmas start at least from 300 K.

For example, a recent thorough analysis²³ of tunneling in the N + O₂ reaction stressed the presence of the tunneling effect, but confined to energy ranges with limited overlap with typical CP applications. Generally speaking, issues related to vibrationally forbidden transitions in the field of CP simulations of practical applications are much more important than other QCT inaccuracies.^28,75

Many semiclassical methods for treating vibrationally inelastic collision processes are designed for weak interactions (typically transitions of few quanta, generally not of QR nature), so they give the best results for PNR transitions by construction. For strong interactions, these methods generally give approximate results and sometimes require a large computational effort that rapidly increases with the amount of energy exchanged. A good strategy is to keep the best of QCT and of a suitable semiclassical method for a given collisional problem, considering that in general QR and PNR collisional events are simultaneously present in the same cross section calculation. In SI is reported Fig. S1 similar to Fig. 1 but in which collision energy is fixed at 3 eV and J is varied in the 0–120 interval, stressing the simultaneous presence of PNR and QR trajectories from low to intermediate J values, of importance in the cross section calculation.

4. Trajectory behavior around QR limit

4.1. An abrupt change of trajectory character

It is useful to introduce a representation of trajectory evolution and the associated potential energy on a plane of normalized internuclear distances , the same used for the definition of CR in eqn (12). While the normalized distance is directly reported on the vertical axis of the plane, the minimum of and is reported on the horizontal axis. This presentation was chosen to focus the representation on the most relevant interactions. The potential energy of each trajectory along its evolution is represented with a color scale shown in eV on the right side of each panel. The thick line in the figure is the minimum energy path (MEP) of the PES, with colors according to the same scale. The MEP color at each point indicates the local potential energy minimum. Trajectories close to the minimum show a color similar to that of the MEP, and higher trajectories show brighter colors. To simplify the trajectory view, two panels are shown, on the left from the trajectory start to a given intermediate time (150 fs in the case of Fig. 2), while on the right panel the subsequent trajectory history is followed from the same time towards the products. This intermediate time is chosen to stop the evolution when the trajectories are near the diagonal of the figure plane along the MEP path (final trajectory points in the left panel are starting points in the right panel). The gray diagonal represents the separation of purely non-reactive trajectories completely on the right of the diagonal from quasireactive trajectories. Trajectories passing this line during their evolution are quasireactive by definition (see eqn (12)), so we call it the QR diagonal. The region around the intersection of the MEP with this diagonal will be called the strong coupling region (SCR).

Fig. 2 Evolution of QR trajectories represented on a simplified plane (see text for details). In the left panel trajectories are shown from reactants to strong coupling region (SCR), on the right panel the same trajectories are shown from SCR to products. Colormap indicates the potential energy associated with each trajectory during time evolution and with the MEP (the thick line). It is worth noting the complex dynamics after the passage through the gray QR diagonal, which is quite far from the reaction barrier for this collisional system.

In Fig. 2, the collision energy is fixed at 4.5 eV and J = 0. In the left panel, the reactant BC diatom motion is clearly visible in oscillating curves. These are very similar to each other but for their vibrational phases, which are uniformly distributed according to the standard QCT method. The curves proceed towards decreasing normalized or distance with similar trends, down to about 1.5–1.6. This first part is characterized by an increase of potential energy as the atom approaches the molecule, with trajectories oscillating at an energy close to the MEP. There is no significant difference between trajectory colors and MEP color along the MEP in this approach region. For values lower than 1.5 of or , the dynamics changes, with brighter trajectory colors than the MEP, not only near the turning points of the molecular vibration, but even when trajectories are passing over the MEP. This means that trajectories are exploring parts of the PES energetically far from the MEP. It is important to recognize that the 3D trajectory view offered by this representation is necessarily partial; away from the asymptotes, the invisible longest distance is comparable with the other two. In the right panel there is a much more complex dynamics from the SCR towards products. Given the high collision energy (4.5 eV, much higher than the reaction threshold of ≈3.26 eV and the dynamical classical threshold at about 4.25 eV at J = 0) and the small number of trajectories considered, all trajectories are quasireactive, i.e. pass the QR diagonal line during their evolution. The minimum potential energy required for reaction is about −6.6 eV (the classical minimum of the NO diatomic potential with respect to the 3-body dissociation limit that defines V = 0; see also Fig. S2 in SI). Practically all the trajectories represented reach peaks higher than this lower limit, but in regions far from the reaction barrier, which is located outside the figure, along the MEP near value 2.5 of normalized coordinate. It is worth noting that the MEP is always increasing from N₂ reactant up to the barrier towards NO formation, with no local minimum, as can be appreciated in Fig. S2. Only one trajectory of the bunch represented in Fig. 2 is reactive, and terminates its evolution by smoothly oscillating around the final part of the MEP and the value 1 on the normalized coordinate , as expected. All the other trajectories show an inelastic quasireactive behaviour, returning to the initial channel after passing the QR diagonal an even number of times. This return in the right panel is much more complex than the approaches in the left panel. The final vibrational elongation is much larger for some trajectories, and globally is much more varying in both amplitudes and phases, demonstrating quasireactive behaviour.

It is now easy to appreciate the crucial difference between left and right panels of Fig. 3, in which the collision energy is 3.5 eV, slightly over the reaction threshold, and only the SCR to products evolution is shown. (The reactants to SCR part is very similar to the left panel of Fig. 2.) The difference between the two panels is in the selection of PNR and QR trajectories for the left and right panels, respectively. In the left panel the returning dynamics at 3.5 eV is as smooth as in the reactants to SCR left panel of Fig. 2. On the contrary, in the right panel the QR trajectories are definitely more unevenly distributed, occasionally with larger final vibrations. From Fig. 2 and 3 it should be clear that the sharp discriminator between complex (QR) and smooth (PNR) behaviour is whether the diagonal of the normalized plane of the figures is passed or not. As already noted, the reaction barrier in this system is quite far from the diagonal. The quasireactivity observed for the current collisional system does not involve any passage through the reaction barrier, but only QR diagonal passing. It is a sort of cis-quasireactivity, in the sense that the quasireactive behaviour is manifested without passing the reaction barrier, with a 2.75 eV QR threshold that is much lower than the 3.26 eV reaction threshold. It is quasireactivity nonetheless, because the final effect on the vibrational distribution is much more similar to the effect of reactivity (wide final vibrational distribution) than to that of PNR collisions (narrow distribution). The maximum potential energies reached by both PNR and QR subsets of trajectories are very similar (same yellow color in Fig. 3) and close to the minimum to reaction (−6.6 eV), so an energy criterion is unable to make any distinctions between PNR and QR features.

Fig. 3 Similarly to Fig. 2, but for two different sets of trajectories, on the left panel PNR only, on the right panel QR only, shown only from SCR to products, at a collision energy of 3.5 eV. It is worth noting the smooth dynamics, very similar to initial propagation (cfr. the left panel of Fig. 2), for PNR trajectories, and the complex dynamics of QR trajectories on the right panel. In fact, the main difference between PNR and QR dynamics consists in the wide final vibrational distribution in the last case as opposed to the very small vibrational variation in the first case.

4.2. Rationalisation of the CR discrimination efficiency

The problem remains of attempting a possible explanation of why collision remoteness appears to be so sharp in vibrational behavior discrimination. In ref. 71 quasireactivity through the reaction barrier in H + D₂ collisions is explained as a sort of tug-of-war between contrasting attractive inner forces of the initial and possible final reactive diatoms. The interpretation is insightful and can be adapted even in the present case, which differs from H + D₂ in the fact that the reaction barrier is never reached, and the whole quasireactive dynamics takes place around the quasireactive diagonal rather than the reaction barrier. In the left panel of Fig. 3 trajectories pass very near the diagonal without going through it, but still manifest a PNR behaviour, while in the right panel trajectories reaching just a bit past the diagonal show a clear QR behaviour. It is this abrupt transition character associated with QR diagonal passage that strongly suggests the presence of a tug-of-war of contrasting forces, in which a very slight advantage of one force over the other determines the final fate of the challenge. The present “vibrational” tug-of-war could be due to the competition between the attractive force into the reactant diatom BC in contrast with the possible attractive force between the initially free atom A and one of the atoms in the BC diatom, the one (temporarily) closer to A. In the following brief discussion, the closest atom to A will be called Q, the other T. In the normalized plane of Fig. 2 and 3, this atomic struggle can be approximately visualized by following the collisional paths. On the QR diagonal, the normalized internuclear distances are equal, so forces acting on the temporary central atom are approximately similar but opposite. A slight movement beyond the QR diagonal determines a significant net attraction between A and Q, temporarily higher than the one between Q and T, and this transient attraction determines the vibrational behavior of QR.

It is worth noting that when the PES is chemically symmetric, as in the H + D₂ collisions of ref. 71 (irrespective of tiny isotope effects), the reaction barrier and the QR diagonal coincide, so the whole quasireactive effect is attributed to the dynamics beyond the barrier and back to reactants (a sort of trans-quasireactivity, which is normally indicated in the literature as quasireactivity). However, given the evidence in this work, the symmetric PES appears more as a special case, while in general the reaction barrier (or more generally the dividing surface) and the QR diagonal are two distinct concepts useful for studying molecular dynamics in different conditions.

5. Probability results

The conditions mainly studied in this work are relevant to the collision process O + N₂ (v = 0, j = 0) → O + N₂ (v′ = 1, j′ ≤ 10) unless otherwise indicated. This is because the interest is in studying vibrational energy transfer by different MD methods, without neglecting final rotation (which is important in chemical kinetics simulations^3,76). On the other hand, effects stemming from very high rotational states should be avoided in this context, for the computational difficulty of treating these states in accurate QM methods and for possible issues in tFHO (in this last case, vibrational energy exchange is essentially independent of rotation, at variance with QCT. When both low and high rotational states are involved in a transition, this can become a problem).

5.1. Application of the standard MD methods to O + N₂ collisions

In Fig. 4 a comparison of J = 0 transition probabilities computed using different dynamics methods is shown. The trajectory-based results (QCT, tFHO and CMM) were obtained as described in Section 2. The two quantum mechanical methods (TICC and TDCC) agree very well from 0.6 eV up to about 3 eV of collision energy. Beyond that limit, the TICC calculation, which is fundamentally limited to the purely inelastic arrangement, does not agree with the TDCC result. Tests indicated secure TICC convergence up to 2.9 eV on the left side of the dip but the basis sets used may not be adequate beyond that collision energy.

Fig. 4 Comparison of different methods used for the calculation of J = 0 probability for the process O + N₂ (v = 0, j = 0) → O + N₂ (v′ = 1, j′ ≤ 10). CMM and tFHO show good to excellent agreement with accurate time-dependent and time-independent QM results up to the dip around 3.1 eV. On the right of the dip, the TICC result may not be converged, while QCT is the only trajectory-based method able to follow the wavepacket probability.

The QCT result appears divided into two disjointed parts, on the left and the right of 3 eV. In the middle, the QCT result is zero. The left branch agrees poorly with the QM data, while the right branch agrees well with the wavepacket (TDCC) result above 3.05 eV. It is worth noting that the absence of continuity along the energy axis is not due to poor statistics. On the contrary, it is a typical manifestation of a classically forbidden process as studied in ref. 74, associated to PNR dynamics. This is exactly the reason for which it is convenient to treat QR events quasiclassically and PNR events semiclassically.

The tFHO probability agrees remarkably well with QM results from low energy up to about 2.75 eV. Fig. S3 in the SI shows that the agreement extends to lower energy and quite low probability values. tFHO is able to reproduce the accurate TICC calculation at low energy because the wavepacket rigid body oscillation hypothesis described in Section 2 is accurate under these conditions. Fig. S3 also displays the well-known convergence problems of the wavepacket calculation at sufficiently low energy. On the other hand, for collision energies higher than about 2.75 eV, tFHO departs from QM. It reproduces the shape of the TDCC curve, but shifted towards higher energy.

The CMM result follows the TDCC probability fairly well from a probability value near 10⁻⁵ up to the collision energy of 3 eV, then shows a deep cuspidal minimum around 3.1 eV followed by a steep increase up to values definitely higher than all the other methods. Because CMM, like QCT, relies on the frequency of events in a trajectory set, it cannot be used for very low probabilities without requiring an unreasonable number of trajectories. It is not expected to work correctly at high energy, at least in the form used here. Preliminary calculations indicate that including the second moment⁴⁴ for application to Δν = 2 does not achieve the same level of accuracy, at least not with the same computational cost used for Δν = 1 in this work.

Summarizing, no trajectory-based method is able to maintain accuracy over the whole energy range explored. The region beyond about 2.75 eV of collision energy appears difficult for all the MD methods used here but the TDCC. QCT is able to reproduce the high energy TDCC trend, but starting from a threshold of 3 eV, while its behavior at lower energy appears totally unreliable.

It is clear that a process different from reaction takes place in the energy region between 2.75 eV and the reaction threshold that disrupts the smoothly increasing low-energy trend. To illuminate this point, Fig. 5 displays the TDCC transition probability as a function of collision energy for specific values of final rotation. Here the change in trend is dramatic in the region with E_coll ≥ 2.75 eV, indicated by the dotted vertical line in the figure, with bumps and dips. The reaction threshold is at 3.26 eV, as shown by the dashed vertical line. The rapid probability oscillations suggest a strong interference effect, probably due to the superposition of the QM QR contribution (intended here as a sort of fast complex formation dynamics) with the QM PNR one (i.e. direct inelastic scattering). Characteristic collision times for PNR/QR/reactive dynamics have been studied in ref. 77. The results are that QR dynamics shows times definitely longer than PNR dynamics and very similar to reactive ones, supporting the hypothesis of QM interference between the two regimes. This abrupt passage to a different dynamics into the “QR region” between 2.75 and 3.26 eV can be interpreted as the emergence of quasireactive dynamics. This emergence can be confirmed from the quasiclassical trajectory dynamics. The same figure shows the sum of the quasireactive, reactive and dissociation quasiclassical probability from the collision O + N₂ (v = 0, j = 0), labeled QRD. (QRD reduces only to the QCT QR probability in the QR region; the reason for adopting this sum will be explained in the next subsection). It is noteworthy that the QRD threshold is exactly 2.75 eV. QRD increases rapidly from its onset up to the reaction threshold, and then increases again beyond this threshold. Quasireaction as defined here is a classical concept linked to well defined internuclear distances, but it appears that its effect can be spotted also in quantum mechanical results, even if it is difficult to define a QM recrossing.⁶⁹

Fig. 5 Probability at J = 0 of the process: O + N₂ (v = 0, j = 0) → O + N₂ (v′ = 1, j′) as a function of collision energy, calculated by TDCC method, compared with the cumulative quasi-reactive, reactive and dissociative (QRD) probability from the same initial state, as obtained from QCT (indeed, dynamical classical reaction threshold is at 4.25 eV for J = 0 and dissociation is at 9.75 eV, so in this specific case QRD reduces to only quasi-reaction). It is insightful the correspondence between the QRD threshold, indicated by a vertical dotted line, and the appearance of the complex whirling of QM probabilities. This effect is likely due to quantum interference effects that can be essentially explained by analogy with the classical direct inelastic scattering from the barrier (without traversing it) and quasi-reactive scattering. However, quasi-reaction does not involve necessarily a reaction barrier crossing (and subsequent recrossing). It is sufficient to cross the QR line as defined in the previous section. In fact, the reaction threshold in the present case is at a much higher value of 3.26 eV, as indicated by the vertical dashed line in the figure.

5.2. Merging SC and QCT

The previous section presented the state of the art concerning standard methods for MD, applied to the problem of interest. Now it is time to introduce the hybrid method, consisting in suitably mixing the results from the trajectory-based methods already presented. This mixing is based on the idea of retaining only QR trajectories in QCT results, obtaining a QCT QR partial result by dividing the number of QR trajectories by the total number of trajectories. Then, something similar will be done with a suitable SC method in a way complementary to QR dynamics, i.e. retaining this time only PNR dynamics. However, more attention is needed in this last case than in the QCT one. We follow ref. 78, extending its assumptions about inelastic processes to PNR events. For QCT, we write the sum of all outcome probabilities as

1 = P_QCT(n₁) = P_QCT-R(n₁) + P_QCT-D(n₁) + P_QCT-QR(n₁) + P_QCT-PNR(n₁) (15)

where n₁ represents the set of initial conditions (vibrational state, rotational state, total angular momentum) while R and D subscripts stand for reactive and dissociative processes, respectively. The first three terms being accurate in QCT in the framework already specified, it is possible to extract the last term as

P_QCT-PNR(n₁) = 1 − (P_QCT-R(n₁) + P_QCT-D(n₁) + P_QCT-QR(n₁)) = 1 − P_QCT-QRD(n₁) (16)

Final vibration is intended as summed wherever it is missing in the argument of P. P_QCT-QRD(n₁) is exactly the same quantity used in Fig. 5. It is the P_QCT-PNR(n₁) term that should be substituted with the SC analogous result. In order to enforce QCT normalization, it is necessary to multiply the normalized probability coming from the SC method times P_QCT-PNR(n₁). In fact, the SC method has its own obvious probability normalization w.r.t. the final vibration n₂:

This normalization has to be valid for any possible outcome of the considered method. For example, tFHO is able to provide a numerical result both for PNR and QR events, but not for reactive and dissociation processes, which in the method simply do not exist (and possible reactive or dissociative trajectories should be dropped). Therefore, the normalization would be

P_SC(n₁) = 1 = P_SC-QR(n₁) + P_SC-PNR(n₁),

where the event frequencies that define the two QR and PNR probabilities are calculated with

P_SC-PNR(n₁) = N_PNR/N_tot (18)

and similarly for P_SC-QR, N_PNR being the number of PNR trajectories and N_tot the total number of trajectories calculated excluding reaction and dissociation. This is not convenient, because the QR part of tFHO is unreliable, being without a physical basis. So, the probability calculation will be operated only including PNR trajectories at the denominator of (18), enforcing:
as if QR trajectories did not exist at all. This probability is normalized by construction relatively to the SC method. When this probability has to be summed to QCT probability in the hybrid method (HY), however, it must be newly normalized within QCT probability, as explained above. Therefore, following and extending⁷⁸ to PNR processes gives:

P_HY-PNR(n₁,n₂) = P_SC-PNR(n₁,n₂) × P_QCT-PNR(n₁) = P_SC-PNR(n₁,n₂) × (1 − P_QCT-QRD(n₁)) (19)

taking P_HY-PNR(n₁,n₂) as the definition of PNR probability in the hybrid semi-quasi-classical method. This definition uses only P_SC-PNR and P_QCT-QRD, which are exactly the most reliable results that can be got for PNR and QR events, respectively. It is important to underline that, with this definition, any SC method could be used if the PNR condition can be correctly imposed.

It is now possible to examine Fig. 6, in which the PNR and QR normalized probability results as just defined in eqn (19) are presented relative to the trajectory-based methods, in comparison with TDCC and TICC. QRD probability is also reported because, as clear from eqn (19), it gradually extinguishes SC probabilities, because it approximates the value 1 around 4 eV, while its threshold is the QR region lower limit. TICC agrees very well with TDCC calculations up to about 3 eV. At 3.1 eV (the dip bottom) TICC appears lower than TDCC. The source of this discrepancy is unclear; TICC basis sets including maximum vibrational quantum numbers of 11 and 13 gave very similar results. For higher values of collision energy, it is not clear that the TICC basis sets are converged.

Fig. 6 Probability results as obtained by QCT-QR, tFHO-PNR and CMM-PNR trajectory-based methods, in comparison with TICC and TDCC results in the same conditions as Fig. 4.

The QCT-QR probability in Fig. 6 reproduces only the right branch of total QCT in Fig. 4, retaining only the result in good agreement with TDCC, showing that the filtered-out QCT-PNR part is totally responsible for the unreliable behavior at lower energy. On the other hand, tFHO-PNR maintains the very good agreement with both quantum calculations up to the QRD threshold, then it is largely supressed by the normalization in eqn (19). Beyond the QRD threshold the QM result is not reproduced very accurately by tFHO-PNR. On the contrary, the CMM-PNR result correctly reproduces the QM result including the QR region up to about 3 eV. This different behavior in the QR region could be due to the slightly different way in which the QR limit is applied in the two cases of tFHO and CMM. CMM uses exactly the same trajectories as the QCT calculation, so the QR limit is applied exactly in the same way. tFHO trajectories start with no vibrational energy, while QCT trajectories start with zero-point vibrational energy. All the other conditions are exactly the same. This difference in initial conditions produces slightly different trajectory dynamics. However, even a slight difference in the QR limit neighbourhood can have a measurable effect, as shown in the previous section. Whether this difficulty using tFHO with QCT can be solved could be the object of further studies. For the time being, it is possible to recover a good agreement of tFHO-PNR with accurate QM probabilities by empirically increasing the QR limit from 1.055 (eqn (13)) to 1.33 (see Fig. S4 in SI). This higher QR limit could be rationalized by invoking the more limited vibrational elongation during tFHO dynamics; QCT dynamics starts with a larger value of vibrational limits. A lower QR limit value has the effect of including QR events into tFHO-PNR calculation, with inaccurate contributions, as already stressed. This point is important but not fundamental in the present discussion, because it pertains to the features of the specific SC method chosen, not to the hybrid method per se. One could improve tFHO-PNR performance in this context or propose another SC method, more similar to QCT. CMM is just QCT with a different final analysis, so it is perfect in this sense, but it cannot be as accurate as tFHO with very low probabilities, while its accuracy at high total energy has not been assessed yet. Other attractive possibilities could be the classical S-matrix⁷⁴ and the frozen Gaussian⁴² semiclassical methods.

In the following, the QR limit is always the standard one, equal to 1.055. In Fig. 7 the merging of QCT-QR with tFHO-PNR and CMM-PNR are shown in comparison with TDCC probability. “QCTFHO” is the sum of QCT-QR with tFHO-PNR probabilities (the small “t” of tFHO has been substituted with capital “T” of QCT, stressing the use of trajectories in both methods), and QCTCMM is the analogous sum for QCT-QR and CMM-PNR results. The very good level of agreement is the same already seen in the previous figure for each partial result in a specific range, because with the specific conditions studied it happens that the two regimes, PNR and QR, are completely disjoint on the left and on the right of the 3 eV value of collision energy. As will be clear in the following, this is not a rule; on the contrary, it is more likely that the two contributions are similar and overlapping. Similar rather than disjoint contributions are practically a rule in cross section calculations that include many different values of J. The novelty in Fig. 7 is in the critical point of junction around 3 eV between SC and QCT methods. The QM dip is approximately reproduced by QCTFHO, which suffers for the approximation introduced by tFHO-PNR in the QR region, while it is approximated more accurately by QCTCMM, for the reasons already explained. It is worth noting that this level of accuracy, stressed also in Table 3, is relative to J = 0 probability, which does not have a great weight in a complete cross section calculation. It will be important in future studies to understand the behavior at much higher J values in this critical region. The general expectation is an easier approximation by classical methods. It is worth noting that the level of agreement between QCTFHO and TDCC results extends also for the transition to v′ = 2, as can be appreciated in Fig. S5 in SI.

Fig. 7 Same as in Fig. 6, but with the hybrid method results of QCTFHO and QCTCMM.

Table 3 Some relevant values of probabilities for J = 0 as obtained using QCTFHO, QCTCMM and TDCC methods

Coll. energy (eV)	QCTFHO	QCTCMM	TDCC
0.5	1.3626 × 10^−8	0	2.56099 × 10^−8
0.6	2.2008 × 10^−7	0	1.85593 × 10^−7
1	9.899 × 10^−5	7.6299 × 10^−5	8.31495 × 10^−5
1.5	0.0021455	0.0022311	0.00253062
2	0.013023	0.01543	0.0151283
2.5	0.032975	0.040323	0.0374921
2.75	0.0405295	0.0461667	0.0404487
3	0.0246301	0.0192888	0.0147775
3.1	0.024216	0.0156354	0.00872165
3.5	0.123841	0.147766	0.109098

---

Fig. 8 shows a comparison of QCTFHO probabilities with QM results at total angular momentum J of 40, 80, and 120. The QR region is stressed by showing the QRD probability in the three cases. It is noteworthy that at J = 120 the QR threshold is around 3.1 eV, much higher than in the other cases examined. The level of agreement with TICC probabilities is excellent, but unfortunately these TICC results are available only at energies lower than 3 eV, so a direct comparison with the full range of QCTFHO data is not possible.

Fig. 8 These three panels show the high level of agreement of QCTFHO with accurate TDCC and TICC calculations for increasing total angular momentum J = 40, 80, 120. Wavepacket method with the coupled states approximation (TDCS), for the cases J = 40 and 80 does not achieve the same level of accuracy in common ranges, so it cannot be taken as a reference for comparison at high energy. The presence of QRD probability marks the QR region in the three cases. It is worth noting the increasing threshold of this probability with J.

In the same figure some TDCC results are shown, as calculated with for J = 80 and for J = 40, because of the huge memory size required by these calculations. The level of agreement with TICC result is excellent at low energy, just slightly degrading at high energy, probably due to the limitation. The general level of agreement of TDCC with QCTFHO results is very good in the common ranges for all the cases examined.

The approximate time-dependent coupled states (TDCS) result is shown in the same figure. It is generally higher than all the QM and trajectory-based results, particularly at high energy. It is clear that TDCS cannot be used as a benchmark, because the discrepancies of QCTFHO with accurate QM results are much smaller than the differences observed using TDCS in the same conditions. The problem of comparison with accurate QM results, which are difficult to obtain in these conditions, becomes the most important limit in the assessment of the new hybrid method. Consideration of lighter collisional systems, more easily treated by accurate QM methods, would not be of help. In those cases QM effects could become more important and provide more difficulty for QCT comparison interpretation, but even successful results would not be useful in the context of cold plasmas including air species.

6. Cross section results

Converged calculations of cross sections have been obtained in this work with QCT, tFHO, QCTFHO and TICC methods (the resources required for the TDCC probability calculation at higher total angular momentum J make it extremely difficult to obtain converged cross sections with that method). The calculations with QCT, tFHO and TICC include J values up to convergence (i.e. from a minimum of 100 to a maximum of 220 of maximum value of J, depending on the collision energy).

6.1. Cross section comparisons with TICC calculations

In Fig. 9 a comparison of QCT, tFHO and TICC cross sections from v = 0, j = 0 to v′ = 1, j′ = 10 is shown as a function of collision energy. In the common ranges, tFHO and TICC show an excellent agreement over many orders of magnitude of cross section values. The QCT result starts to be different from zero only at the QR threshold of 2.75 eV. QCT and tFHO in the high energy region starting from about 3 eV are in fairly good agreement. However, this fact could be accidental, due to the low reliability of tFHO when quasireactivity is in action as already discussed.

Fig. 9 Comparison of QCT, tFHO and TICC cross sections from v = 0, j = 0 to v′ = 1, j′ = 10. The level of agreement between TICC and tFHO is excellent in the common ranges.

Other cases including more final rotational states and for different final vibration is shown in SI, with similar very good comparisons, in Fig. S6 (j′ ≤ 64) and Fig. S7 (v′ = 2).

6.2. Higher vibrational states treated with the hybrid method: expected impact on cold plasmas

In Fig. 10 QCT, tFHO and QCTFHO results from v = 5, j = 0 to v′ = 4, j′ ≤ 64 are presented. The tFHO result shows the typical ascending trend already seen in the previous examples, very steep at low energy, with the slope decreasing gradually up to a plateau around 3.5 eV. The QCT cross section starts from about 1.3 eV and rapidly becomes significantly higher than tFHO. The two curves then converge around 4 eV. QCTFHO strictly follows tFHO up to about 2.5 eV, then it becomes somewhat lower. QCTFHO is the sum of QCT QR and tFHO PNR, so it can occasionally be lower than QCT and tFHO results, as in this case. (In this section tFHO-PNR is always QRD-normalized, as explained in the previous section (eqn (19).) The coincidence of QCTFHO with tFHO for a large interval, in particular when QCT is higher than tFHO, means that QCT in that region is essentially given by QCT PNR, which is filtered out in the QCTFHO sum by construction. The clear prevalence of tFHO in this case is due to considering a small vibrational gap Δv = v′ − v from a relatively high initial vibration (v = 5), where normally PNR effects are expected to be the highest. On the other hand, a QCT result null up to a threshold and then significantly higher than more accurate calculations is typical in these conditions, as shown by different authors.^45,78

Fig. 10 Comparison of QCT, tFHO and QCTFHO cross section results from v = 5, j = 0 to v′ = 4, j′ ≤ 64. When the vibrational gap is low, QCTFHO tends to follow tFHO result.

However, QCTFHO behavior becomes quite different with increasing Δv, as shown in Fig. 11. Here for Δv = 2 QCTFHO is very similar to QCT in the energy interval between 2.25 and 3 eV, then becomes closer to tFHO. When Δv = 3 from v = 5, as in Fig. 12, the result is more similar to Fig. S7, i.e. QCTFHO tends to be mainly given by tFHO at low energy (and low probability) and QCT at high energy (and high probability). This explains also the relative success concerning vibrational energy exchange of standard QCT at high energy, where the maximum part of the probability is given by the accurate QCT QR contribution rather than the unreliable QCT PNR.

Fig. 11 Comparison of QCT, tFHO and QCTFHO cross section results from v = 5, j = 0 to v′ = 3, j′ ≤ 64. Here QCTFHO is intermediate and can hardly be deduced only on the basis of QCT and tFHO trends.

Fig. 12 Comparison of QCT, tFHO and QCTFHO results from v = 5, j = 0 to v′ = 2, j′ ≤ 64. For v–v′ difference sufficiently high, as in this case, the two contributions of QCT and tFHO are approximately summed in QCTFHO, with the exclusion of the high energy region.

In the same Fig. 12, QCT-QR and tFHO-PNR are also shown. The straight line trend for tFHO-PNR from 2 to 4 eV contrasts with the rapid increase of the corresponding tFHO in the same interval. This increase is entirely due to the tFHO-QR contribution that brings tFHO result close to QCT only at high energy, showing once again that tFHO-QR is generally not reliable. The QCT result, however, is not due only to QCT-QR: starting from about 3.5 eV on, even a QCT-PNR contribution is visible (by difference between QCT and QCT-QR). This means that PNR contributions are not limited to low energy. Beyond a first threshold region in which QR regime dominates, the PNR contribution cannot be ignored at high energy. This aspect could also qualitatively explain the success of FHO calculations (different from tFHO, see below in this section) at high energy, as in ref. 10 and 11.

These examples demonstrate that QCTFHO cannot be easily deduced from QCT and tFHO trends, and that QCT-QR and tFHO-PNR cross sections are needed instead. The QCT-QR threshold, which appears to have an important role in QCTFHO result, cannot be easily deduced without specific calculations; it can be very different from reaction threshold, as in this case. Preliminary calculations show that this behavior continues, and is even increasingly important, for higher initial vibration. This means that the complex interplay between PNR and QR kinds of behavior studied in this work depends on total energy, not simply on collision energy.

The foreseeable consequences in CP kinetics modeling are relevant. In fact, the input data of the most accurate air CP models essentially are taken from state-to-state QCT and FHO calculations.^76,79 FHO^10,11 is an analytic method based on a forced harmonic oscillator model with simplified hypothesis with respect to the more general tFHO, but in a first approximation the two methods should provide at least qualitatively similar results. The results presented in this work suggest that both QCT and tFHO (and as a consequence even FHO) are only partially correct over large total energy ranges of vibrational energy exchange processes in molecular collisions, exactly those ones fundamental in CP kinetics. By carefully selecting partial results of those methods, however, it seems possible to extract accurate data, that are very difficult and extremely expensive to calculate with the most accurate methods available. Of course, this work is just scraping the surface of the problem. Research is now needed to assess the range of validity of this hybrid method in the wide variety of conditions of interest in CP, including the treatment of different species combinations and high-lying rovibrational states.

7. Conclusions and perspectives

In non-equilibrium plasmas, a central role is played by vibrational energy transfer in molecular collisions. Modeling those plasmas requires accurate data as input, with a complete level of detail concerning fluence of vibrational energy among the molecular species involved, including also the highest vibrational levels, while rotation must be considered at least in an approximate way. This in turn requires a reliable method for calculating these input data, because experiments cannot in general provide such a level of vibrational detail. Accurate QM methods are available, but hardly can provide all the data needed with reasonable amounts of computational resources, in particular when heavy particles like air species are involved. In this work, quasiclassical and semiclassical trajectory results have been merged in order to overcome the limitations of each method, and get a result that achieves significant level of agreement with theoretical calculations from the most accurate QM methods available. This merging originates from studying the dynamical characteristics associated with a quantity, the collision remoteness, that appears to provide a useful guide to the trajectory behavior in vibrational energy transfer. With this metric it is possible to discriminate trajectories with a behavior that can accurately be treated with QCT from others that can be successfully treated with a suitable semiclassical method. The results at the level of probabilities and cross sections are very encouraging, with very accurate comparisons, paving the way to the possible application of the “hybrid” semi-quasi-classical method to the extensive calculations required by CP data for modeling. Some minor inaccuracies are attributed to the slightly different treatment of initial conditions in molecular collisions in the semiclassical and quasiclassical methods. On this point more research is needed concerning the most suited SC method to adopt. Other SC methods different from tFHO and CMM used in this work can be proposed for the purpose, with the only requirements that the method permit the accurate calculation of collision remoteness similarly to QCT, and accurately treat purely non-reactive processes. For example, the classical S-matrix theory⁷⁴ or a suitable modification of the frozen gaussians⁴² method could likely be adapted to work in these conditions, overcoming the limitations of tFHO and CMM. The possibility of accurate, detailed and affordable calculations for vibrational energy transfer can be a game-changer in CP modeling, where only approximate data are routinely used,²⁸ with inaccuracies that can be as large as orders of magnitude, in particular with air species.

Quasireaction is normally associated with a collisional system in which there is a missed reaction event, i.e. an event in which the reaction barrier is passed an even number of times, with the system eventually returning back to the initial channel after a reactive-like dynamics. In this work it has been shown that in the more general case of the non-symmetric system O + N₂, the collisional system shows a quasireactive dynamics by just passing the QR diagonal as defined in this work, rather than the reaction barrier. The QR diagonal coincides with the reaction barrier for symmetric systems.

The sum of quasireactive, reaction and dissociation probabilities from classical trajectories has been shown to be in strict connection with the emergence of a distinctly complex behavior in quantum-mechanical probabilities, which can be identified with the QM appearance of quasireaction, which is intrinsically difficult to define exactly in QM terms. This correspondence can be useful in correctly attributing such complex behavior in different contexts.

Other possibilities to explore concern the use of final vibrational distributions “sliced” into more than two groups using CR, for assignment to multiple different methods. For example, tFHO can be very accurate in weak interactions (i.e. high collision remoteness) involving low-lying vibrational states, but a more appropriate SC method could be used in the intermediate regime between tFHO PNR and QCT QR regimes, where some inaccuracies have been shown. Interestingly, final vibrational distributions show complex structures, partially visible in Fig. 1 and Fig. S1, but more clearly visualized as functions of rotation and vibration that were not shown in this work. They deserve further study in order to understand their origin and possibly exploit them to help select the most appropriate dynamical methods.

Author contributions

F. Esposito: conceptualization, formal analysis, investigation, methodology, resources, software, supervision, visualization, writing – original draft, review & editing; N. Bulut: formal analysis, investigation, resources, writing – review & editing; P. Żuchowski: resources, review & editing; G. McBane: formal analysis, investigation, resources, software, visualization, writing – original draft, review & editing.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp03853g.

Appendices

Microscopic reversibility applied to tFHO probabilities/cross sections

Microscopic reversibility requires that the probability of a given forward transition must be the same as the corresponding reverse one, so for a given total energy E_tot:³⁹

P_forward(E_tot) = P_reverse(E_tot) (A1)

Now, E_tot can be written:

E_tot = E_ki + E_i = E_kf + E_f (A2)

where E_ki is the kinetic energy for the forward collision, with E_kf final kinetic energy, and analogously E_i is the initial internal energy, E_f the final internal energy. However, trajectories in tFHO start with a given classical kinetic energy E_kcl. In order to recover microscopic reversibility, E_kcl should be the average of initial and final kinetic energies:

E_kcl = (E_ki + E_kf)/2 (A3)

As a consequence:

E_kcl = (2E_tot − E_i − E_kf)/2 (A4)

This can be solved for E_ki and E_kf respectively as:

E_ki = E_kcl + ΔE/2 (A5)

as in eqn (7), and

E_kf = E_kcl − ΔE/2

The values of E_ki and E_kf can be used for substituting E_tot on the left and on the right of eqn (A1). The solution proposed in the original INDECENT model¹³ and in many others semiclassical approaches in which vibration is treated quantum mechanically while translation is classical, considers the average of the relative speeds of reagents and products. In the case of ref. 39 and 40, the relation between E_ki, E_kf and E_kcl adopted is similar to (A5). We tried that relation, but it does not bring to any significant difference with respect to eqn (7) in the comparisons with accurate QM calculations, so we reverted to the simpler eqn (5), which has the great advantage of simplifying the interpolation described in Section 2.

Acknowledgements

F. E. thanks Qhizen Hong, Cecilia Coletti and Dmitri Babikov for useful discussions on the semiclassical method of G. D. Billing (propaedeutic to further studies in connection with the present work). N. B. is grateful for the support provided by the Polish National Agency for Academic Exchange (NAWA) and the Nobelium program ‘Joining Gdańsk Tech Research Community,’ grant number DEC-1/12025/IDUB/I.1a/No: 038122. The ADEP.24.09 project is also gratefully acknowledged. P. S. Z. is grateful to National Science Centre (NCN) of Poland for funding the grant “Sonata Bis 9” No. 2019/34/E/ST4/00407. We also would like to thank COST action (CA21101) “Confined Electronic Systems”. G. C. M. is grateful to Heriot-Watt University for hospitality in 2024 and 2025.

References

1. M. Capitelli, R. Celiberto, G. Colonna, F. Esposito, C. Gorse, K. Hassouni, A. Laricchiuta and S. Longo, Fundamental Aspects of Plasma Chemical Physics, Springer New York, New York, NY, 2016.
2. Plasma Modeling: Methods and applications, ed. G. Colonna and A. D’Angola, IOP Publishing, 2nd edn, 2022.
3. L. D. Pietanza, O. Guaitella, V. Aquilanti, I. Armenise, A. Bogaerts, M. Capitelli, G. Colonna, V. Guerra, R. Engeln, E. Kustova, A. Lombardi, F. Palazzetti and T. Silva, Advances in non-equilibrium CO₂ plasma kinetics: a theoretical and experimental review, Eur. Phys. J. D, 2021, 75, 237,  DOI:10.1140/epjd/s10053-021-00226-0.
4. D. Aceto, P. F. Ambrico and F. Esposito, Air cold plasmas as a new tool for nitrogen fixation in agriculture: underlying mechanisms and current experimental insights, Front. Phys., 2024, 12, 1455481,  DOI:10.3389/fphy.2024.1455481.
5. A. Fridman, Plasma Chemistry, Cambridge University Press, 2008.
6. A. Hurlbatt, A. R. Gibson, S. Schröter, J. Bredin, A. P. S. Foote, P. Grondein, D. O’Connell and T. Gans, Concepts, Capabilities, and Limitations of Global Models: A Review, Plasma Process. Polym., 2017, 14, 1600138,  DOI:10.1002/ppap.201600138.
7. R. T. Pack, Space-fixed vs body-fixed axes in atom-diatomic molecule scattering. Sudden approximations, J. Chem. Phys., 1974, 60, 633–639,  DOI:10.1063/1.1681085.
8. P. McGuire and D. J. Kouri, Quantum mechanical close coupling approach to molecular collisions. j_z-conserving coupled states approximation, J. Chem. Phys., 1974, 60, 2488–2499,  DOI:10.1063/1.1681388.
9. R. D. Sharma and C. A. Brau, Energy Transfer in Near-Resonant Molecular Collisions due to Long-Range Forces with Application to Transfer of Vibrational Energy from ν₃ Mode of CO₂ to N₂, J. Chem. Phys., 1969, 50, 924–930,  DOI:10.1063/1.1671145.
10. I. V. Adamovich, S. O. Macheret, J. W. Rich and C. E. Treanor, Vibrational Energy Transfer Rates Using a Forced Harmonic Oscillator Model, J. Thermophys. Heat Transfer, 1998, 12, 57–65,  DOI:10.2514/2.6302.
11. M. Lino da Silva, J. Loureiro and V. Guerra, Nonequilibrium dissociation and recombination rates in nitrogen: From shock waves to discharge conditions, Chem. Phys., 2009, 358, 123–131,  DOI:10.1016/j.chemphys.2009.01.007.
12. G. D. Billing, Classical path method in inelastic and reactive scattering, Int. Rev. Phys. Chem., 1994, 13, 309–335,  DOI:10.1080/01442359409353298.
13. C. F. Giese and W. R. Gentry, Classical trajectory treatment of inelastic scattering in collisions of H⁺ with H₂, HD, and D₂, Phys. Rev. A:At., Mol., Opt. Phys., 1974, 10, 2156.
14. D. Babikov and A. Semenov, Recent Advances in Development and Applications of the Mixed Quantum/Classical Theory for Inelastic Scattering, J. Phys. Chem. A, 2016, 120, 319–331,  DOI:10.1021/acs.jpca.5b09569.
15. D. G. Truhlar and J. T. Muckerman, Reactive scattering cross sections: quasiclassical and semiclassical methods, Atom-Molecule Collision Theory: A Guide for the Experimentalist, ed. R. B. Bernstein, Plenum Press, New York, 1979.
16. D. E. Mihalik, R. Wang, B. H. Yang, P. C. Stancil, T. J. Price, R. C. Forrey, N. Balakrishnan and R. V. Krems, Accurate machine learning of rate coefficients for state-to-state transitions in molecular collisions, J. Chem. Phys., 2025, 162, 024116,  DOI:10.1063/5.0242182.
17. D. Bossion, G. Nyman and Y. Scribano, Machine learning prediction of state-to-state rate constants for astrochemistry, Artif. Intell. Chem., 2024, 2, 100052,  DOI:10.1016/j.aichem.2024.100052.
18. J. C. San Vicente Veliz, J. Arnold, R. J. Bemish and M. Meuwly, Combining Machine Learning and Spectroscopy to Model Reactive Atom + Diatom Collisions, J. Phys. Chem. A, 2022, acs.jpca.2c06267,  DOI:10.1021/acs.jpca.2c06267.
19. S. Akpinar, I. Armenise, P. Defazio, F. Esposito, P. Gamallo, C. Petrongolo and R. Sayós, Quantum mechanical and quasiclassical Born–Oppenheimer dynamics of the reaction N₂ + O → N + NO on the N₂O a³A″ and b³A′ surfaces, Chem. Phys., 2012, 398, 81–89,  DOI:10.1016/j.chemphys.2011.05.005.
20. F. Esposito, C. M. Coppola and D. De Fazio, Complementarity between Quantum and Classical Mechanics in Chemical Modeling. The H + HeH⁺ → H₂⁺ + He Reaction: A Rigorous Test for Reaction Dynamics Methods, J. Phys. Chem. A, 2015, 119, 12615–12626,  DOI:10.1021/acs.jpca.5b09660.
21. F. Esposito and I. Armenise, Reactive, Inelastic, and Dissociation Processes in Collisions of Atomic Oxygen with Molecular Nitrogen, J. Phys. Chem. A, 2017, 121, 6211–6219,  DOI:10.1021/acs.jpca.7b04442.
22. R. Celiberto, M. Capitelli, G. Colonna, G. D’Ammando, F. Esposito, R. Janev, V. Laporta, A. Laricchiuta, L. Pietanza, M. Rutigliano and J. Wadehra, Elementary Processes and Kinetic Modeling for Hydrogen and Helium Plasmas, Atoms, 2017, 5, 18,  DOI:10.3390/atoms5020018.
23. F. Esposito, P. Gamallo, M. González and C. Petrongolo, Exploring quantum tunneling in heavy atom reactions using a rigorous theoretical approach to the dynamics: formation of NO + O from the N + O₂ atmospheric reaction, Phys. Chem. Chem. Phys., 2025, 10.1039.D5CP00539F,  10.1039/D5CP00539F.
24. M. Capitelli, R. Celiberto, G. Colonna, F. Esposito, C. Gorse, K. Hassouni, A. Laricchiuta and S. Longo, Reactivity and Relaxation of Vibrationally/Rotationally Excited Molecules with Open Shell Atoms, Fundamental Aspects of Plasma Chemical Physics, Springer, New York, 2016, pp. 31–56.
25. F. Esposito, Reactivity, relaxation and dissociation of vibrationally excited molecules in low-temperature plasma modeling, Rend. Lincei Sci. Fis. Nat., 2019, 30, 57–66,  DOI:10.1007/s12210-019-00778-9.
26. F. Esposito, R. Macdonald, I. D. Boyd, K. Neitzel and D. A. Andrienko, Heavy-particle elementary processes in hypersonic flows, Hypersonic Meteoroid Entry Physics, IOP Publishing, 2019, pp. 16-1–16-32.
27. F. Esposito and I. Armenise, Reactive, Inelastic, and Dissociation Processes in Collisions of Atomic Nitrogen with Molecular Oxygen, J. Phys. Chem. A, 2021, 125, 3953–3964,  DOI:10.1021/acs.jpca.0c09999.
28. F. Esposito, On the relevance of accurate input data for vibrational kinetics in air cold plasmas: the case of nitrogen fixation, Plasma Sources Sci. Technol., 2022, 31, 094010,  DOI:10.1088/1361-6595/ac9082.
29. I. V. Adamovich, F. Esposito and S. O. Macheret, Rate coefficients for energy transfer and chemical reactions in heavy particle collisions, in Plasma Modeling Methods and applications, ed. G. Colonna and A. D’Angola, IOP Publishing, 2nd edn, 2022, pp. 23-1–23-30.
30. P. Gamallo, M. González and R. Sayós, Ab initio derived analytical fits of the two lowest triplet potential energy surfaces and theoretical rate constants for the N(⁴S) + NO(X ²Π) system, J. Chem. Phys., 2003, 119, 2545,  DOI:10.1063/1.1586251.
31. Q. Hong, M. Bartolomei, F. Esposito, C. Coletti, Q. Sun and F. Pirani, Reconciling experimental and theoretical vibrational deactivation in low-energy O + N₂ collisions, Phys. Chem. Chem. Phys., 2021, 23, 15475–15479,  10.1039/D1CP01976G.
32. Q. Hong, M. Bartolomei, F. Pirani, F. Esposito, Q. Sun and C. Coletti, Vibrational deactivation in O(³P) + N₂ collisions: from an old problem towards its solution, Plasma Sources Sci. Technol., 2022, 31, 084008,  DOI:10.1088/1361-6595/ac86f3.
33. F. Esposito, Dinamica quasiclassica di processi collisionali inelastici e reattivi in sistemi H + H₂ e N + N₂ rotovibrazionalmente risolti, PhD thesis, Università degli Studi di Bari, Bari (Italy), 1999.
34. F. J. Aoiz, V. Sáez-Rábanos, B. Martínez-Haya and T. González-Lezana, Quasiclassical determination of reaction probabilities as a function of the total angular momentum, J. Chem. Phys., 2005, 123, 094101,  DOI:10.1063/1.2009739.
35. F. Esposito and M. Capitelli, Quasiclassical trajectory calculations of vibrationally specific dissociation cross-sections and rate constants for the reaction O + O₂(v) → 3O, Chem. Phys. Lett., 2002, 364, 180–187,  DOI:10.1016/S0009-2614(02)01329-5.
36. L. Bonnet and J.-C. Rayez, Gaussian weighting in the quasiclassical trajectory method, Chem. Phys. Lett., 2004, 397, 106–109,  DOI:10.1016/j.cplett.2004.08.068.
37. N. C. Blais and D. G. Truhlar, Ab initio calculation of the vibrational energy transfer rate of H₂ in Ar using Monte Carlo classical trajectories and the forced quantum oscillator model, J. Chem. Phys., 1978, 69, 846,  DOI:10.1063/1.436600.
38. C. E. Treanor, Vibrational Energy Transfer in High-Energy Collisions, J. Chem. Phys., 1965, 43, 532,  DOI:10.1063/1.1696777.
39. G. Billing, The semiclassical treatment of molecular roto/vibrational energy transfer, Comput. Phys. Rep., 1984, 1, 239–296,  DOI:10.1016/0167-7977(84)90006-6.
40. A. Semenov, M. Ivanov and D. Babikov, Ro-vibrational quenching of CO (v = 1) by He impact in a broad range of temperatures: A benchmark study using mixed quantum/classical inelastic scattering theory, J. Chem. Phys., 2013, 139, 074306,  DOI:10.1063/1.4818488.
41. E. H. Kerner, Note on the forced and damped oscillator in quantum mechanics, Can. J. Phys., 1958, 36, 371–377.
42. E. J. Heller, Time-dependent approach to semiclassical dynamics, J. Chem. Phys., 1975, 62, 1544–1555,  DOI:10.1063/1.430620.
43. H.-D. Meyer, On the forced harmonic oscillator with time-dependent frequency, Chem. Phys., 1981, 61, 365–383,  DOI:10.1016/0301-0104(81)85155-5.
44. R. J. Gordon, A comparison of exact classical and quantum mechanical calculations of vibrational energy transfer, J. Chem. Phys., 1976, 65, 4945–4957,  DOI:10.1063/1.432971.
45. R. J. Gordon, A comparison of exact classical and quantum mechanical calculations of vibrational energy transfer. II. The effect of long-lived collision complexes, J. Chem. Phys., 1977, 67, 5923–5929,  DOI:10.1063/1.434799.
46. D. G. Truhlar, B. P. Reid, D. E. Zurawski and J. C. Gray, Vibrational energy transfer and an improved information-theoretic moment method. Comparison of the accuracy of several methods for determining state-to-state transition probabilities from quasiclassical trajectories, J. Phys. Chem., 1981, 85, 786–791.
47. J. M. Hutson and S. Green, MOLSCAT computer code, version 14, distributed by Collaborative Computational Project No. 6 of the Engineering and Physical Sciences Research Council, Swindon, UK, 1994.
48. J. Light, Inelastic Cross Sections: Theory, in Atom-Molecule Collision Theory: A Guide for the Experimentalist, ed. R. B. Bernstein, Plenum, New York, 1979.
49. J. M. Hutson and C. R. Le Sueur, MOLSCAT: A program for non-reactive quantum scattering calculations on atomic and molecular collisions, Comput. Phys. Commun., 2019, 241, 9–18,  DOI:10.1016/j.cpc.2019.02.014.
50. R. J. Le Roy, Y. Huang and C. Jary, An accurate analytic potential function for ground-state N₂ from a direct-potential-fit analysis of spectroscopic data, J. Chem. Phys., 2006, 125, 164310.
51. R. J. Le Roy, LEVEL: A computer program for solving the radial Schrödinger equation for bound and quasibound levels, J. Quant. Spectrosc. Radiat. Transfer, 2017, 186, 167–178,  DOI:10.1016/j.jqsrt.2016.05.028.
52. C. C. Marston and G. G. Balint-Kurti, The Fourier grid Hamiltonian method for bound state eigenvalues and eigenfunctions, J. Chem. Phys., 1989, 91, 3571,  DOI:10.1063/1.456888.
53. M. H. Alexander and D. E. Manolopoulos, A stable linear reference potential algorithm for solution of the quantum close-coupled equations in molecular scattering theory, J. Chem. Phys., 1987, 86, 2044–2050,  DOI:10.1063/1.452154.
54. M. Alagia, N. Balucani, L. Cartechini, P. Casavecchia, G. Gualberto Volpi, F. Javier Aoiz, L. Bañares, T. C. Allison, S. L. Mielke and D. G. Truhlar, Dynamics of the Cl + H₂/D₂ reaction: a comparison of crossed molecular beam experiments with quasiclassical trajectory and quantum mechanical calculations, Phys. Chem. Chem. Phys., 2000, 2, 599–612,  10.1039/a908829f.
55. S. Gómez-Carrasco and O. Roncero, Coordinate transformation methods to calculate state-to-state reaction probabilities with wave packet treatments, J. Chem. Phys., 2006, 125, 054102,  DOI:10.1063/1.2218337.
56. A. Zanchet, T. González-Lezana, A. Aguado, S. Gómez-Carrasco and O. Roncero, Nonadiabatic state-to-state reactive collisions among open-shell reactants with conical intersections: the OH(²Π) + F(²P) example, J. Phys. Chem. A, 2010, 114, 9733.
57. Y. Huang, D. J. Kouri and D. K. Hoffman, General energy separable Faber polynomial representation of operator functions: Theory and application in quantum scattering, J. Chem. Phys., 1994, 101, 10493.
58. V. A. Mandelshtam and H. S. Taylor, A simple recursion polynomial expansion of the Green's function with absorbing boundary conditions. Application to the reactive scattering, J. Chem. Phys., 1995, 103, 2903.
59. Y. Huang, S. S. Iyengar, D. J. Kouri and D. K. Hoffman, Further analysis of solutions to the time independent wave packet equations of quantum dynamics. II. Scattering as a continuous function of energy using finite, discrete approximate Hamiltonians, J. Chem. Phys., 1996, 105, 924.
60. G. C. Kroes and D. Neuhauser, Performance of a time-independent scattering wave packet technique using real operators and wave functions, J. Chem. Phys., 1996, 105, 8690.
61. R. Chen and H. Guo, Evolution of quantum system in order domain of Chebyshev operator, J. Chem. Phys., 1996, 105, 3569.
62. S. K. Gray and G. G. Balint-Kurti, Quantum dynamics with real wave packets, including application to three-dimensional (J = 0)D + H₂ → HD + H reactive scattering, J. Chem. Phys., 1998, 108, 950,  DOI:10.1063/1.475495.
63. T. González-Lezana, A. Aguado, M. Paniagua and O. Roncero, Quantum approaches for the insertion dynamics of the H⁺ + D₂ and D⁺ + H₂ reactive collisions, J. Chem. Phys., 2005, 123, 194309.
64. O. Roncero and P. del Mazo-Sevillano, MADWAVE3: A quantum time dependent wave packet code for nonadiabatic state-to-state reaction dynamics of triatomic systems, Comput. Phys. Commun., 2025, 308, 109471,  DOI:10.1016/j.cpc.2024.109471.
65. I. W. Smith, Reaction and relaxation of vibrationally excited molecules: a classical trajectory study of Br + HCl(v′) and Br + DCl(v′) collisions, Chem. Phys., 1977, 20, 437–443.
66. I. W. Smith, Relaxation in collisions of vibrationally excited molecules with potentially reactive atoms, Acc. Chem. Res., 1976, 9, 161–168.
67. I. W. M. Smith, Vibrational relaxation in atom exchange reactions: A classical trajectory study of H + H₂(ν) and D + D₂(ν) collisions using an accurate potential, Chem. Phys. Lett., 1977, 47, 219–224,  DOI:10.1016/0009-2614(77)80005-5.
68. I. W. M. Smith and P. M. Wood, Vibrational relaxation in atom-exchange reactions: A classical, Monte Carlo, trajectory study, Mol. Phys., 1973, 25, 441–454,  DOI:10.1080/00268977300100381.
69. F. J. Aoiz, L. Bañares and V. J. Herrero, The H + H₂ reactive system. Progress in the study of the dynamics of the simplest reaction, Int. Rev. Phys. Chem., 2005, 24, 119–190,  DOI:10.1080/01442350500195659.
70. A. L. Brunsvold, H. P. Upadhyaya, J. Zhang, R. Cooper, T. K. Minton, M. Braunstein and J. W. Duff, Dynamics of Hyperthermal Collisions of O(³P) with CO, J. Phys. Chem. A, 2008, 112, 2192–2205,  DOI:10.1021/jp710025v.
71. S. J. Greaves, E. Wrede, N. T. Goldberg, J. Zhang, D. J. Miller and R. N. Zare, Vibrational excitation through tug-of-war inelastic collisions, Nature, 2008, 454, 88–91,  DOI:10.1038/nature07079.
72. N. T. Goldberg, J. Zhang, K. Koszinowski, F. Bouakline, S. C. Althorpe and R. N. Zare, Vibrationally inelastic H + D₂ collisions are forward-scattered, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 18194–18199,  DOI:10.1073/pnas.0807942105.
73. J. Jankunas, M. Sneha, R. N. Zare, F. Bouakline, S. C. Althorpe, D. Herráez-Aguilar and F. J. Aoiz, Is the simplest chemical reaction really so simple?, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 15–20.
74. W. H. Miller, Classical S Matrix: Numerical Application to Inelastic Collisions, J. Chem. Phys., 1970, 53, 3578,  DOI:10.1063/1.1674535.
75. G. D’Ammando, M. Capitelli, F. Esposito, A. Laricchiuta, L. D. Pietanza and G. Colonna, The role of radiative reabsorption on the electron energy distribution functions in H₂/He plasma expansion through a tapered nozzle, Phys. Plasmas, 2014, 21, 093508,  DOI:10.1063/1.4895481.
76. V. Guerra, A. Tejero-del-Caz, C. D. Pintassilgo and L. L. Alves, Modelling N₂–O₂ plasmas: volume and surface kinetics, Plasma Sources Sci. Technol., 2019, 28, 073001,  DOI:10.1088/1361-6595/ab252c.
77. D. Koner, J. C. San Vicente Veliz, R. J. Bemish and M. Meuwly, Accurate reproducing kernel-based potential energy surfaces for the triplet ground states of N2O and dynamics for the N + NO ↔ O + N₂ and N₂ + O → 2N + O reactions, Phys. Chem. Chem. Phys., 2020, 22, 18488–18498,  10.1039/D0CP02509G.
78. J. C. Gray, G. A. Fraser, D. G. Truhlar and K. C. Kulander, Quasiclassical trajectory calculations and quantal wave packet calculations for vibrational energy transfer at energies above the dissociation threshold, J. Chem. Phys., 1980, 73, 5726–5733,  DOI:10.1063/1.440053.
79. E. Vervloessem, M. Aghaei, F. Jardali, N. Hafezkhiabani and A. Bogaerts, Plasma-Based N₂ Fixation into NO_x: Insights from Modeling toward Optimum Yields and Energy Costs in a Gliding Arc Plasmatron, ACS Sustainable Chem. Eng., 2020, 8, 9711–9720,  DOI:10.1021/acssuschemeng.0c01815.

---