Comparison of intermolecular energy transfer from vibrationally excited benzene in mixed nitrogen–benzene baths at 140 K and 300 K

Understanding and accurately modeling collisional intermolecular energy transfer (IET) is important for characterizing reaction rates in nonequilibrium environments, such as photoinitiated reactions, high temperature combustion systems, and after super/hypersonic shock waves.¹ This is especially important in systems where IET rates are comparable to chemical reaction rates.² Studying IET in nonequilibrium systems has often been performed by monitoring the efficiency of energy transfer of a vibrationally excited molecule of interest to relatively colder bath molecules.^3–9 These studies often involved measuring the average energy transferred per collision, ⟨∆E_c⟩.⁴ Previous studies have shown a strong dependence on the size of a bath molecule for the relative efficiency of deactivating collisions.^10,11 As a result of intermolecular interaction, large bath molecules may stay close to a high energy donor for longer periods of time, leading to multiple energy exchange events between the bath and donor molecules. These multiple interactions result in more efficient IET compared to collisions with small bath molecules.⁵ IET rates also strongly depend on the total number of available pathways for vibration–vibration (V-V) and vibration–rotation/translation (V-R/T) energy transfer for a given molecular collision pair.^9,12

Non-equilibrium between Boltzmann distributions of different types of molecular motion exists in supersonic and hypersonic gas flow fields and is known as Thermal Non-Equilibrium (TNE).^13–16 TNE is often generated by shock waves and in gas expansion since, in these examples, the gas in the flow experiences a sharp change in pressure before there are enough molecular collisions to re-thermalize different types of molecular internal energy. The rate of re-equilibration of molecular translation, rotation, vibration, and electronic excitation does not occur on the same timescale, and often IET can occur over the length scale of a flow field, such as over the body of a hypersonic vehicle. Previous work has shown that this IET can modify macroscopic flow properties, such as turbulent fluctuations in the velocity for flow fields.^17–23 

Turbulence in supersonic and hypersonic boundary layers is often caused by acoustic fluctuations in a flow field.^17,18 It has been demonstrated that if molecules in the flow can absorb the acoustic energy through IET of the correct magnitude and timescale, then the onset of transition to a turbulent flow field will be delayed.²² It is favorable in many supersonic and hypersonic flow fields to delay the transition from a laminar flow field to a turbulent flow field. Turbulence in flow fields around supersonic and hypersonic vehicles increases drag and rate of heating of the vehicle. This requires more expensive vehicle exteriors that can either withstand the increased heat or sacrificially ablate when they become too hot in order to protect the interior of the vehicle.

It is not known to what extent TNE can be utilized to modify acoustic fluctuations that cause turbulence in supersonic and hypersonic boundary layers, but several studies have begun to connect TNE and the modification of turbulence.^17–22 The first studies to examine the role of TNE in modifying turbulence focused on delaying turbulent transitions in hypersonic boundary layers by addition of varying percentages of CO₂ into the flow fields.^17–21 Some of these studies examined flow fields with CO₂ pre-mixed into the gas, while others examined the effect of injecting CO₂ into the boundary layer of a hypersonic vehicle.

Since the rate that CO₂ was able to absorb acoustic fluctuations in flow fields caused the delay of transition to turbulence, studies were extended to the injection of other types of molecules into the same flow fields.²² The injection of He, N₂, and C₄F₈ into the supersonic flow field by Schmidt and Shepherd²² demonstrated that the onset of transition to turbulence was increased in the flow field with He injection and delayed in the flow field with C₄F₈ injection as compared to a flow field with N₂ injection. Therefore, tuning molecular factors that change the rate and magnitude of IET such as molecular weight and the number of vibrational modes were shown to allow for determination of the right rate of IET to modify the transition to turbulence in this flow field. This shows that there are optimal IET rates for absorbing specific frequencies of acoustic fluctuations, and therefore, the ability to tune IET for diverse flow field conditions is desirable. Details of IET rates can also have a profound effect on fully developed turbulent flows, altering both turbulent fields and average distributions of energy across molecular modes²³ and their spectral distributions.¹³ This was found to depend, among other things, on the ratio of IET and flow timescales. TNE has also been found to affect the decay of turbulent flows in the experimental work of Fuller et al.¹⁵ and later reproduced and explained by simulations,¹⁶ which found that the interaction of TNE and turbulence depends on both the degree of TNE and the disparity of IET and timescales. This provides venues for both energy management and flow control in hypersonic applications.

Classical chemical dynamics simulation models have been developed to study IET from a vibrationally excited molecule in a bath of molecules,^{9–12,14,24–45} e.g., N₂ bath.^14,38–40 With these simulation models, IET may be studied for mixed baths to provide insights into the dynamics associated with turbulence for multiple gases in a flow field. In a recent work,⁶ IET was studied for a vibrationally excited benzene molecule, C₆H₆^*, in a mixed N₂ and C₆H₆ bath at 300 K, since unexcited C₆H₆ remains in experiments where C₆H₆^* is generated with UV light. It was found that C₆H₆–C₆H₆^* vibration-to-vibration (V-V) IET is a key process for vibrational re-thermalization following excitation. V-V pathways for C₆H₆–C₆H₆^* can result in significantly more efficient rates of IET than models that only consider N₂–C₆H₆^* IET pathways.^5,7,12 Equilibration was found to be slow for C₆H₆^* in the N₂/C₆H₆ bath. For a 1 atm bath pressure and at 10⁻⁷ s, there were four distinct temperatures, i.e., a rotational–translational (RT) temperature for C₆H₆^* and the N₂/C₆H₆ bath, a N₂ vibrational temperature, and different vibrational temperatures for C₆H₆^* and the C₆H₆ bath.

In the work presented here, the above simulation for C₆H₆^* vibrational relaxation in the N₂/C₆H₆ bath at 300 K is extended to consider a much colder bath at 140 K. Of interest is to determine how the IET dynamics may vary in a much colder bath, which is relevant to the temperatures studied in supersonic and hypersonic wind tunnels in which TNE is most often studied by aerospace engineers in ground-based blowdown facilities. The 140 K simulation results are compared with the results of an experimental study at 140 K.⁴⁶ 

N₂/N₂, C₆H₆/C₆H₆, and N₂/C₆H₆ intermolecular potentials were required for the simulations. The same potentials were used, as used in previous simulations, and are only briefly described here. They are written as sums of two-body potentials. The N₂/N₂ potential was developed from MP2 calculations extrapolated to the CBS limit.³⁸ The C₆H₆/C₆H₆ potential is represented by the Optimized Potentials for Liquid Simulations (OPLS) model.⁴⁷ The OPLS potential gives an overall good description of the benzene–benzene interaction.⁴⁸ The geometry for global potential energy minima of the benzene dimer is tilted-T, which is in excellent agreement with CCSD(T)/CBS calculations.^48,49 The OPLS Bz–Bz center-of-mass distance is 4.93 Å, whereas the CCSD(T)/CBS value is 4.96 Å. Moreover, the energy that the OPLS model gives for the potential energy minimum is −2.32 kcal/mol compared to the CCSD(T)/CBS value of −2.84 kcal/mol. The N₂/C₆H₆ two-body potentials are written as V(r) = A exp(−Br) + C/rⁿ + D/r^m. The values of the parameters A, B, C, D, n, and m for the N₂/C₆H₆ C–N and H–N interactions were assumed to be the same as those for the azulene + N₂ potential, developed from SCS-MP2/6-311++G^** calculations.⁴⁰ The N₂ intramolecular potential was represented by a Morse function, with parameters taken from the experiment.⁵⁰ The benzene intramolecular potential was represented by C–C and C–H Morse stretches, C–C–C and C–C–H harmonic bends, C–C–H harmonic wags, and torsions and is the potential used in previous studies of Na⁺ interacting with benzene.^51,52 The potential gives benzene vibrational frequencies in good agreement with the experiment.⁵³ 

Three different simulations were performed for the results presented here. One was performed to model experiments in which C₆H₆ constitutes 1.25% of the C₆H₆–N₂ system, approximately 16%–30% of the C₆H₆ molecules were vibrationally excited, and the bath was initially at 140 K. To represent these experiments, a simulation model was used in which the system consists of 396 N₂ and 3 C₆H₆ bath molecules, with one C₆H₆ molecule vibrationally excited. Thus, the simulation system was 1% C₆H₆, with 25% of the C₆H₆ molecules vibrationally excited. The initial bath temperature for this simulation was 140 K. To compare this simulation with the previous C₆H₆–N₂ bath simulation at 300 K,⁶ two additional simulations were performed. For each, the system consisted of 1 excited C₆H₆ molecule, 0 unexcited bath C₆H₆ molecules, and baths of 400 and 399 N₂ molecules at 300 K and 140 K, respectively.

To achieve the binary/single collision limit for comparison with experiments, the bath density was chosen as 40 kg/m³ or 16.2 atm, which was found to be the binary/single collision limiting density for C₆F₆ + N₂ simulations.^38,39 Due to the fact that C₆H₆ is a smaller molecule than C₆F₆, and also that there is a very small percentage of C₆H₆ in the bath, the binary/single collision limiting density is expected to be achieved at 40 kg/m³ or higher density. Performing the simulations in the binary/single collision limit allows for extrapolation of the simulation results to lower densities/pressures. However, such extrapolation may not be very accurate in the case when the possibility of complex formation among the molecular species present in the system is there. Such a possibility may reduce with the lowering bath pressure, which may lead to different IET dynamics.

The simulations were performed with the same methodology as described for previous intermolecular energy transfer bath simulations.^14,38–42 A vibrational energy of 148.1 kcal/mol was first added to one excited C₆H₆^* molecule to model the experimental 193 nm laser excitation and subsequent internal conversion.^6,46 This energy was added randomly with classical microcanonical normal mode sampling,^54,55 as implemented in a modified version³⁸ of the general chemical dynamics computer code VENUS.⁵⁶ Translational and rotational energies for 140 K or 300 K were then added randomly from a Boltzmann distribution at those corresponding temperatures to vibrationally excited C₆H₆^*. With initial conditions for C₆H₆^* chosen, the next step was to equilibrate the bath around C₆H₆^* by placing it at the center of the simulation box with its coordinates and momenta fixed. An MD simulation was then performed to thermally equilibrate the bath to the desired temperature (i.e., 140 K or 300 K), using periodic boundary conditions and nearest neighbor updating to enhance the simulation. At the end of this equilibration, the desired initial temperature for the vibration, rotation, and center-of-mass translation degrees of freedom for bath molecules was verified. To illustrate this equilibration, consider the simulation bath with 396 N₂ molecules and 3 C₆H₆ molecules at 140 K. After equilibration, the average center-of-mass translation energy for N₂ and C₆H₆ molecules was 3RT/2 = 0.41 kcal/mol, the average rotational energy for N₂ and C₆H₆ molecules was RT = 0.28 kcal/mol and 3RT/2 = 0.41 kcal/mol, respectively, and the average vibrational energy for N₂ and each mode of C₆H₆ was RT = 0.28 kcal/mol. These are the proper equilibrium average energies, which match the equipartition model.

With the above random initial conditions for C₆H₆^* and the bath, a trajectory was then calculated for 3 ns to study intermolecular energy transfer from C₆H₆^* to the bath. To obtain results that could be compared with experimental results, averaging was performed by calculating an ensemble of 39 trajectories, with random initial conditions. In a previous similar simulation for N₂ + C₆F₆ intermolecular energy transfer,³⁸ 48 trajectories gave statistically the same result as found for 96 trajectories. A simulation with only 24 trajectories gave semi-quantitative results.

In the simulation with C₆H₆^* in the N₂/C₆H₆ bath at 140 K, the center-of-mass distances between the benzene molecules were monitored to study the possibilities of dimer, trimer, and tetramer long-lived collision complex formation between the benzene molecules. There were four benzene molecules, one excited and three unexcited. In Figs. 1 and 2, distances between the benzene molecules are plotted vs time for two trajectories. The excited benzene molecule is Bz1, with the other three Bz2, Bz3, and Bz4 unexcited, and there are six center-of-mass distances between these molecules. Two simultaneous criteria were used to identify complex formation between Bz molecules, i.e., the distance between the centers-of-mass of two Bz molecules is less than 10 Å (the equilibrium distance between Bz molecules is 4.93 Å for the Bz dimer)⁴⁸ and the two Bz molecules retain a distance less than 10 Å for at least 10 ps. The orientation averaged potential energy curve as done in a previous work⁴⁸ showed a potential energy of −0.039 kcal/mol when the center-of-mass separation between two Bz molecules is 10 Å (see Fig. S2 with a related discussion in the supplementary material). The period for the Bz–Bz intermolecular harmonic stretching vibration is 0.5 ps for the Bz dimer. However, the lowest normal mode frequency of the Bz dimer is 6.4 cm⁻¹, corresponding to a vibrational period of ∼5 ps. For dimer formation, a reasonable time criterion could be two Bz–Bz vibrational periods. Moreover, fitting the orientation averaged Bz–Bz intermolecular potential energy curve to a Morse function of the form D_e[1 − exp(β(r − r₀))]² provides D_e = 0.62 kcal/mol, β = 2.88 Å⁻¹, and r₀ = 6.58 Å. With the help of these parameters, the Bz–Bz vibrational period is obtained as ∼1 ps at the orientation averaged potential energy of −0.05 kcal/mol. Therefore, the selection of 10 ps as a criterion for dimer formation is at least two times larger than both the Bz dimer vibrational frequency and any of the Bz vibrational frequencies. Figure 1 is a plot of the six Bz–Bz distances for the trajectory with a minimum amount of complex formation, while Fig. 2 is the plot for the trajectory with the maximum amount of complex formation.

Averaged over the 39 trajectories and with the criteria for benzene complex formation given above, the average time a 3000 ps trajectory spent as a complex was 612 ps. Thus, 20.4% of the time, the trajectory consisted of a benzene complex. The percentage of time the complex was a dimer was 17.8%, while the percentages as a trimer and tetramer were much smaller and 2.25% and 0.34%, respectively. When only complexes involving the excited benzene (Bz1) are considered, these percentages became 7.90% as total benzene complexes, 6.72% as a dimer, 1.01% as a trimer, and 0.17% as a tetramer. For the sample trajectory depicted in Fig. 1, the total percentage of complex formation was 7.06% and all are dimers, while in Fig. 2, the total percentage of complex formation was 37.1%. Out of this, 13.8% was a dimer, 11.7% was a trimer, and 11.6% was a tetramer.

Average energies of C₆H₆^* and the N₂ and C₆H₆ bath molecules, for the 140 K N₂/C₆H₆ bath simulation, are plotted vs time in Figs. 3 and 4. Figure 3 gives the C₆H₆^* total energy, while individual translational, rotational, and vibrational energies of the bath molecules are given in Fig. 4. For a C₆H₆ bath molecule, the average vibrational energy for one mode of the molecule is given, which equals the average vibrational energy of a C₆H₆ bath molecule divided by 30, the number of vibrational modes. While IET efficiency from different C₆H₆ modes has been previously shown to not be equal,⁵⁷ this method of calculating the per mode vibrational energy allows for comparison with individual energies of rotational and translational modes of the bath. In the following, average energies are given for the 39 trajectories, with uncertainties as standard deviations of the mean. For a N₂ bath molecule, the initial average translation, rotation, and vibration energies, i.e., ⟨E_trans⟩, ⟨E_rot⟩, and ⟨E_vib⟩, are 0.416 kcal/mol ± 0.003 kcal/mol, 0.278 kcal/mol ± 0.002 kcal/mol, and 0.276 kcal/mol ± 0.002 kcal/mol, respectively, and at the 3 ns conclusion of the trajectories, these average energies are 0.470 kcal/mol ± 0.003 kcal/mol, 0.313 kcal/mol ± 0.003 kcal/mol, and 0.277 kcal/mol ± 0.002 kcal/mol. Within statistical uncertainty, the average initial temperatures for translation, rotation, and vibration are each 140 K, and at the conclusion of the trajectories, the average temperatures for these respective degrees of freedom are 158 K ± 1 K, 157 K ± 2 K, and 140 K ± 1 K. The heating of translation and rotation is similar, while there is no energy transfer to N₂ vibration, which is consistent with the energy partitioning in similar previous calculations that had an initial bath temperature of 298 K.³⁸ 

For the benzene (Bz) bath molecule, in the N₂/C₆H₆ bath, the initial ⟨E_trans⟩ and ⟨E_rot⟩ are 0.445 kcal/mol ± 0.03 kcal/mol and 0.433 kcal/mol ± 0.04 kcal/mol, respectively, and the average energy in one vibration mode of Bz, ⟨E_vib⟩, is 0.270 kcal/mol ± 0.01 kcal/mol. At the 3 ns conclusion of the trajectories, these average energies are 0.450 kcal/mol ± 0.04 kcal/mol, 0.560 kcal/mol ± 0.04 kcal/mol, and 1.08 kcal/mol ± 0.04 kcal/mol. Within statistical uncertainty, the average initial temperatures for translation, rotation, and vibration are each 140 K, and at the conclusion of the trajectories, the average temperatures for these respective degrees of freedom are 150 K ± 13 K, 188 K ± 13 K, and 543 K ± 20 K. There are two interesting features for the energy transfer to Bz. Energy transfer to Bz vibration is substantially more efficient than transfer to Bz translation and rotation, dynamics observed in the previous simulation of C₆H₆^* relaxation in a N₂/C₆H₆ bath initially at 300 K.⁶ In addition, energy transfer to Bz rotation is slightly more efficient than transfer to Bz translation. The fluctuations in the average Bz rotational energy are much greater than those for the translational energy. Such dynamics were not observed in the previous N₂/C₆H₆ bath simulation at 300 K.

⟨E_trans⟩ and ⟨E_rot⟩ and per mode vibrational energy ⟨E_vib⟩ of C₆H₆^* vs time, for the N₂/C₆H₆ bath simulation, are given in Fig. 5. The initial ⟨E_trans⟩ and ⟨E_rot⟩ are 0.45 kcal/mol ± 0.05 kcal/mol and 0.48 kcal/mol ± 0.07 kcal/mol, respectively, and ⟨E_vib⟩ is 4.93 kcal/mol ± 0.40 kcal/mol. At the 3 ns conclusion of the trajectories, these average energies are 0.49 kcal/mol ± 0.06 kcal/mol, 0.67 kcal/mol ± 0.14 kcal/mol, and 1.70 kcal/mol ± 0.14 kcal/mol. ⟨E_trans⟩ does not change with time, a result similar to that for a C₆H₆ bath molecule (Fig. 4 and above discussion). ⟨E_rot⟩ rapidly increases but then decreases to a value slightly higher than that for ⟨E_trans⟩. For a C₆H₆ bath molecule, ⟨E_rot⟩ increases with time (Fig. 4). The large fluctuations in ⟨E_rot⟩, for C₆H₆^* in Fig. 5, are partly due to the small number of C₆H₆^* molecules in the analysis, i.e., 39, with one for each of the 39 trajectories. At the end of the 3 ns simulation, the C₆H₆^* translational, rotational, and vibrational temperatures are 164 K ± 20 K, 225 K ± 47 K, and 855 K ± 70 K.

The 3 ns translational, rotational, and vibrational temperatures for the N₂ and C₆H₆ bath molecules, and C₆H₆^*, are summarized in Table I. It should be noted that a complete re-equilibration of the C₆H₆^* vibration has not been achieved at 3 ns, i.e., at the termination of the trajectories. The final C₆H₆^* per mode vibrational energy at 3 ns is 1.70 kcal/mol ± 0.14 kcal/mol, which is about 1/3 of its initial energy. The 3 ns energy of C₆H₆ bath molecules for the same mode is 1.08 kcal/mol ± 0.04 kcal/mol.

The average energies of C₆H₆^* vs time, for the simulations with the N₂ bath at 140 K and 300 K, are plotted in Fig. 6. These results are obtained from the simulations where there is no unexcited C₆H₆ molecule in the bath. One can see from this figure that the energy transfer at 300 K is more efficient than that at 140 K. The translational, rotational, and vibrational temperatures of these simulations at both 140 K and 300 K are obtained and compared with the ones from mixed bath simulations mentioned above. At the 3 ns conclusion of the 140 K simulations, T_trans and T_rot are 155 K ± 1 K and 157 K ± 1 K, respectively, for a N₂ bath molecule, and 137 K ± 16 K and 319 K ± 74 K for C₆H₆^*. T_vib for C₆H₆^* is 1974 K ± 25 K. These temperatures, except T_vib, are statistically the same as those given above for the N₂/C₆H₆ bath simulation at 140 K. For the 300 K simulation of the N₂ bath, T_trans and T_rot are 316 K ± 2 K and 317 K ± 2 K for N₂ and 310 K ± 40 K and 362 K ± 52 K for C₆H₆^*. T_vib for C₆H₆^* is 1749 K ± 52 K. Within statistical uncertainties, these temperatures, except T_vib, are the same, as found for the previous N₂/C₆H₆ bath simulation at 300 K.⁶ However, for both the 140 K and 300 K N₂ bath simulations, there is some indication that T_rot for C₆H₆^* is a bit higher than T_trans and T_rot for N₂ and T_trans for C₆H₆^*. However, better statistics are required to determine if this is indeed the case.

An important analysis from the simulations is determination of the average energy transferred per collision, ⟨ΔE_c⟩, from vibrationally excited C₆H₆^*. ⟨ΔE_c⟩ is found from the simulation ⟨E(t)⟩ and given by

where d⟨E(t)⟩/dt is the energy transferred per unit time and ω is the collision frequency in s⁻¹. To facilitate the analysis, the simulation ⟨E(t)⟩ are fit analytically. As found for previous simulations,^6,48–42 the ⟨E(t)⟩ are well fit by the bi-exponential

where f₁ + f₂ = 1, E(0) is the initial energy of C₆H₆^*, k₁ and k₂ are rate constants, and E(∞) is the corresponding energy value of C₆H₆^* at complete re-equilibration of the baths. The fits are shown in Figs. 3 and 6, with the fitting parameters listed in Table II, along with properties of the baths. Also included are the fitting parameters for the previous N₂/C₆H₆ bath simulation at 300 K.⁶ ⟨∆E_c⟩ in Eq. (1) includes all collisions, both of those that transfer the energy from and to C₆H₆^*.

For the simulation with the N₂/C₆H₆ bath, the collision frequency is the sum of the collision frequencies for C₆H₆^* colliding with the N₂ and C₆H₆ bath molecules and is ω = ω(C₆H₆^* − N₂) + ω(C₆H₆^* − C₆H₆). For the N₂ bath simulations, the only collision frequency is ω(C₆H₆^*–N₂). The collision frequency for each bath component may be expressed as ω = ω_P × P, where P is the partial pressure of the bath gas. The ω_P used here is the same as the values to interpret the experiments,³² with which the simulations are compared. These values for N₂ and C₆H₆ along with the resulting ω(C₆H₆^*–C₆H₆) and ω(C₆H₆^*–N₂) are summarized in Table III for both 140 K and 300 K mixed bath simulations. The pressures used to determine the collision frequencies were those for the bath when the simulation was initiated. The partial pressures of N₂ and C₆H₆ for the 140 K and 300 K mixed bath simulations are also presented in Table III. The pressure for the pure N₂ bath is 35.0 atm and 16.3 atm at 300 K and 140 K, respectively. It is interesting to note here that numbers of C₆H₆^*–N₂ and C₆H₆^*–C₆H₆ collisions in 2.4 ns of the trajectory calculation at 300 K are ∼847 and ∼62, respectively, whereas the numbers for the 3 ns trajectory integration at 140 K are 1161 and ∼13, respectively, yielding a lower percentage of C₆H₆^*–C₆H₆ collisions at 140 K.

Of interest from the analysis of ⟨E(t)⟩, as described by Eq. (1), are the values of ⟨ΔE_c⟩ for C₆H₆^* energy transfer to the N₂ and C₆H₆ bath molecules, i.e., ⟨ΔE_c⟩_N2 and ⟨ΔE_c⟩_Bz. Values of ⟨ΔE_c⟩_N2 for collisions of N₂ with C₆H₆^*, at 140 K and 300 K, were determined from the simulation plots of ⟨E(t)⟩ in Fig. 6 and the fitting parameters in Table II. The resulting values of ⟨ΔE_c⟩_N2 vs ⟨E(t)⟩ are plotted in Fig. 7, where it is seen that ⟨ΔE_c⟩_N2 is significantly smaller at 140 K than 300 K. The same analysis was applied to ⟨E(t)⟩ vs time in Fig. 3, using the fitting parameters in Table II, to determine ⟨ΔE_c⟩ vs ⟨E(t)⟩ for collisions of C₆H₆^* with the N₂/C₆H₆ bath at 140 K. This information was given previously for the simulations at 300 K. The plots of ⟨ΔE_c⟩ vs ⟨E(t)⟩, for the 140 K and 300 K simulations, are given in the supplementary material. Using Eq. (A6) in the Appendix, the values of ⟨ΔE_c⟩ for collisions of C₆H₆^* with the N₂/C₆H₆ bath, and values of ⟨ΔE_c⟩_N2, values of ⟨ΔE_c⟩_Bz were determined for collisions of C₆H₆^* with the C₆H₆ bath molecules. The values are plotted in Fig. 7, and it is seen that C₆H₆^* + C₆H₆ energy transfer is much more efficient at 140 K than 300 K. The ⟨ΔE_c⟩_Bz at 140 K is ∼11 kcal/mol, whereas the same at 300 K is ∼5 kcal/mol at the excitation energy of 148.1 kcal/mol. Note here that in a previously done classical trajectory study⁵⁷ on Bz^* + Bz single-collision energy transfer, the ⟨ΔE_c⟩_Bz was obtained as ∼2.6 kcal/mol at a Bz^* excitation energy of 116.4 kcal/mol and Bz temperature of 300 K. However, it was also shown that when the rotationally frozen condition was implemented in the classical trajectory study, the chance of formation of the collisional complexes with longer lifetimes was larger, and consequently, the ⟨ΔE_c⟩_Bz value increased. At the same excitation energy, ⟨E(t)⟩_Bz^* = 116.4 kcal/mol, and the ⟨ΔE_c⟩_Bz from our simulations in Fig. 7 is ∼3.7 kcal/mol, which could be due to the formation of more collisional complexes in the present simulation model. However, at 140 K, ⟨ΔE_c⟩_Bz becomes almost two times larger at ⟨E(t)⟩_Bz^* = 116.4 kcal/mol than at 300 K. In the previous classical trajectory study⁵⁸ mentioned above, the lowest Bz temperature simulated was 200 K, where ⟨ΔE_c⟩_Bz was found to be ∼4% more than the one at 300 K and was ∼2.7 kcal/mol for the down Bz^*–Bz collisions. Therefore, it can be assumed that the value of ⟨ΔE_c⟩_Bz increases sharply when the temperature decreases from 200 K to 140 K. The possible reason for such a behavior is discussed in Sec. V.

For comparisons with the experiments, described below, it is important to have simulation values for the average energy transferred from vibrationally excited Bz^* per collision, ⟨ΔE_c⟩_Bz^*, vs the bath temperature. The atomistic simulation studies for the model presented here provide the details of IET dynamics, namely, the effect of intermolecular potential energy parameters, energy distribution of the excited molecule, and energy partitioning among center-of-mass translational, rotational, and vibrational modes for the collision partners, along with the underlying dynamics. It is also possible to identify the gateway modes of energy transfer. Overall, the simulations provide many more IET parameters than experiments, which only measured the bath N₂ transient rotational/translational temperature by proxy with 1% seeded NO. Previously,⁶ a simulation was performed for Bz^* relaxation in a N₂/Bz bath at 300 K, and for the current work, such a simulation was performed at 140 K. However, the ⟨ΔE_c⟩_Bz^* vs ⟨E(t)⟩_Bz^* simulation results do not provide a direct comparison since the composition of the N₂/Bz bath was different for the two simulations. The model for the 300 K simulation consisted of 1 excited benzene molecule Bz^* and 9 Bz and 190 N₂ bath molecules, to represent experiments with a N₂/Bz bath of 5% Bz, with 10% of Bz vibrationally excited.⁶ In contrast, the model for the current 140 K simulation consisted of a Bz^* molecule and 3 Bz and 396 N₂ bath molecules, constituting a model that is 1% Bz, with 25% of Bz vibrationally excited. This model represents experiments at 140 K in which Bz constitutes 1.25% of the N₂/Bz system and 16%–30% of Bz is vibrationally excited.

Although the total energy transfer from Bz^*, ⟨ΔE_c⟩_Bz^*, may not be directly compared for the 140 K and 300 K simulations, the values for energy transfer vs temperature to the N₂ and Bz bath molecules, ⟨ΔE_c⟩_N2 and ⟨ΔE_c⟩_Bz, provide a means to determine how ⟨ΔE_c⟩_Bz^* depends on temperature. As published previously,⁶ ⟨ΔE_c⟩_Bz^* = 0.44 kcal/mol for Bz^* with an average energy ⟨E(t)⟩_Bz^* = 150 kcal/mol in the 300 K simulation. In order to compare simulated ⟨ΔE_c⟩_Bz^* values at 140 K and 300 K for Bz^* with ⟨E(t)⟩_Bz^* = 150 kcal/mol, the bath composition at 140 K (i.e., 396 N₂ + 3 C₆H₆) was converted to the bath composition at 300 K (i.e., 190 N₂ + 9 C₆H₆) by utilizing the values of ⟨ΔE_c⟩_N2 and ⟨ΔE_c⟩_Bz in Fig. 7 for 140 K in Eq. (A6). The resulting value for ⟨ΔE_c⟩_Bz^* is 0.83 kcal/mol, about two times larger than the value at 300 K. The same comparison may be made for the simulation at 140 K, for which ⟨ΔE_c⟩_Bz^* = 0.18 kcal/mol at ⟨E(t)⟩_Bz^* = 150 kcal/mol. Using the bath composition for 140 K in Eq. (A6), with values of ⟨ΔE_c⟩_N2 and ⟨ΔE_c⟩_Bz in Fig. 7 for 300 K, ⟨ΔE_c⟩_Bz^* is 0.12 kcal/mol, which is almost 1.5 times smaller than the value at 140 K. Values of ⟨ΔE_c⟩_Bz^* vs ⟨E(t)⟩_Bz^*, for the 140 K and 300 K simulations, are given in Fig. S1, and as described above, values of ⟨ΔE_c⟩_Bz^* may be determined for the 140 K and 300 K bath compositions at 300 K and 140 K, respectively. The resulting plots of ⟨ΔE_c⟩_Bz^* at 140 K and 300 K, for the 140 K and 300 K bath compositions, are given in Fig. 8 vs ⟨E(t)⟩_Bz^*. ⟨ΔE_c⟩_Bz^* is approximately 1.5 to 2 times larger at 140 K at all ⟨E(t)⟩_Bz^* for both baths.

The experimental procedures were described in detail previously,^7,46 and the experiments at 300 K were summarized in a comparison of simulation and experimental studies of Bz^* relaxation in a N₂/Bz bath initially at 300 K.⁶ In the experiments, Bz^* was generated by 193 nm photoexcitation and subsequent rapid internal conversion, and then Bz^* was collisionally relaxed in N₂/Bz baths of 300 K and 140 K initial temperatures. The transient rotational–translational temperature rise of (97.75%) N₂ in the bath was monitored by proxy via LIF rotational temperature measurements of seeded NO, which comprised 1% of the bath. Seeded NO was chosen due to its quick RT–RT equilibration with N₂, a negligibly different heat capacity to N₂ (to not disturb low-temperature Laval flow fields), and the relatively high S/N obtained from LIF measurements. Temperature-dependent Bz^*–N₂ and Bz^*–Bz collision frequencies were calculated with Lennard-Jones parameters for Bz and N₂. For unambiguous comparison with the experiment, the same collision frequencies were used to analyze the simulations. The previous comparisons of the initially 300 K experiments and simulations gave values for ⟨ΔE_c⟩_Bz^* in agreement, and the simulations revealed the importance of Bz^*–Bz vibrational energy transfer in the collisional relaxation of Bz^*.

For the 140 K experiment considered here, the N₂/Bz bath consisted of 1.25% Bz, and ∼16%–30% of Bz was estimated to have been vibrationally excited. The temperature rise of the bath reached 95% of its maximum value, ∼225 K, by 10 µs, which was calculated to have occurred over ∼400 Bz^*–N₂ and ∼8 Bz^*–Bz collisions. After ∼10 µs, the lack of any further significant observable rise in NO rotational temperature indicated that either Bz^* had finished collisional relaxation or any further relaxation occurred on a timescale that was too slow to measure in the Laval apparatus. In Fig. 8, a plot of the experimental average energy transferred per collision vs the internal energy of Bz^* is given and compared with the above simulation results with similar bath composition and percentage of Bz molecules vibrationally excited. As seen from Fig. 8, the experimental ⟨ΔE_c⟩_Bz^* is about 14 times larger than the simulation ⟨ΔE_c⟩_Bz^* at the ⟨E(t)⟩_Bz^* of 150 kcal/mol. There is a discrepancy in final N₂ rotational temperature between the experiment and simulation. The ΔT from the experiment is about five times larger than that from the simulation. The possible sources of this mismatch between experimental and simulation ⟨ΔE_c⟩_Bz^* are addressed in Sec. V.

By comparing the current experimental results with those from previous studies,³² the average energy transferred per collision was found to increase with the decreasing temperature below 300 K, which is opposite of the behavior observed above 300 K. The IET efficiency for the 140 K curve, where the temperature rises from 140 K to 225 K, is approximately two times larger than the IET efficiency for the temperature rise from 300 K to 610 K. Therefore, the IET rate increases more quickly per kelvin below 300 K than it does above 300 K. A few possible origins of the more efficient energy transfer at 140 K are discussed in Sec. V.

While the percentage of Bz molecules was up to 5.6 times lower in the 140 K experiment than in the 300 K experiments, it is important to consider whether both the extent to which Bz^*–Bz and Bz–Bz collision lifetimes increase at lower temperatures due to longer lived “chattering” collisions and the possible formation of collisionally stabilized dimers at low temperatures would allow for more efficient Bz^*–Bz IET. It may be useful to remember that based on an orientation-averaged calculation,⁴⁸ the Bz–Bz potential energy parameters used in the simulation provide a potential energy minimum of −0.54 kcal/mol for a Bz–Bz center-of-mass separation of 6.7 Å. However, the global minimum was calculated as −2.32 kcal/mol for a tilted-T geometry. The formation of long-lived complexes in the simulations of this work is consistent with the observations from simulations of Bernshtein and Oref in which ∼7 ps “chattering” collision complexes formed in classical trajectory calculations at 200 K (as opposed to ∼30 fs “hit-and-run” collisions).⁵⁸ Therefore, these “chattering” collision complexes would have longer collision lifetimes at 140 K due to the lower average molecular velocities, and the increased IET efficiency can be explained by multiple energy transfer events during each “chattering” collision complex.

Along with the Bz^*/Bz–Bz complex formations, one could also expect the possible formation of Bz^*/Bz–N₂ complexes during the simulation at low T. If such complexes were formed, the IET could be enhanced to some extent. Many trajectories were analyzed with a similar method to evaluate the presence of Bz dimers, i.e., by monitoring the center-of-mass distance between Bz and N₂. Several random N₂ molecules were selected from the bath to calculate the center-of-mass distance with all four Bz molecules. However, there is no evidence of any Bz^*/Bz–N₂ complex formation in the simulation. An orientation average potential energy vs center-of-mass distance as utilized for Bz–Bz complexes was also calculated for Bz–N₂ interactions using the potential energy parameters used in the simulation (see Fig. S3 in the supplementary material). The minimum energy was obtained as −0.24 kcal/mol, which is more than two times larger than that of Bz–Bz. Moreover, the global minimum for the Bz–N₂ system was only −0.75 kcal/mol and is much smaller than −2.32 kcal/mol of the Bz dimer. Thus, the Bz–N₂ interaction is much weaker than the Bz–Bz interaction, which is perhaps the reason for the Bz–N₂ complexes not being formed in the simulation.

In order to avoid the formation of collisionally stabilized benzene dimers in experiments, experimental conditions were chosen in a regime in which no significant formation of collisionally stabilized benzene dimers would have been able to occur on the timescale of the 140 K experiment according to experimental rate coefficients for the formation of collisionally stabilized benzene dimers measured by Hamon et al.⁵⁹ (i.e., the reason that a lower [Bz] was utilized for the 140 K experiment than for the 300 K experiments). The experimental conditions were also similar to those of previous NO quenching experiments by Bz.⁶⁰ Since benzene has a melting point of 278.68 K and is a solid at 140 K given a long enough timescale, experiments were carried out such that the Bz/N₂/NO gas mixture was cooled and measured quickly enough at a low enough pressure that negligible formation of collisionally stabilized benzene dimers would have been seen (where collisions with other bath species occurred on the order of every 100 ns during a 10 µs experiment). However, if this assumption is untrue, then it is one possible explanation for the dramatically increased IET observed in the 140 K experimental results. A recent study on the simulated IET dynamics of the C₆H₆–C₆F₆^* complex in the N₂ bath showed that the IET efficiency increased significantly for aromatic complexes.⁶¹ At the excitation energy of 67 kcal/mol and 125 kcal/mol, the ⟨ΔE_c⟩ was obtained as ∼2.7 kcal/mol and ∼4.8 kcal/mol, respectively, for the N₂ bath density of 20 kg/m³. By comparison, for an excitation energy of 104 kcal/mol, C₆F₆^* gave a ⟨ΔE_c⟩ of ∼2.7 kcal/mol at the same bath density.^38,39 If collisionally stabilized benzene dimers (Bz–Bz) were present during the experiments, then some Bz^*–Bz IET would not have needed to wait for the collision of a Bz^* monomer and a Bz monomer if dimers were excited by the 193 nm light. In this case, excited Bz dimers as well as Bz oligomers would have occurred more often than originally anticipated, both of which would have led to increased IET. Based on Hamon’s dimerization rate coefficient that was measured at a much lower pressure,⁵⁹ the percent of complex at 123 K would be much less than 1% during the current 140 K simulation time of 3 ns. There is a negligible percentage of the Bz–Bz complex formed in the simulation with vibrational temperature below 300 K and no complex with vibrational temperature similar to the overall bath temperature. Alternatively, if the dimers with lifetime 50 ps or more are taken, the percentage of those dimers formed in the simulation is about 2% with respect to the total simulation time. Thus, considering collisionally stabilized Bz–Bz dimers and not just Bz–Bz collisional complexes, those numbers are more-or-less consistent with that of Hamon’s.

Although the IET dynamics are non-statistical within the simulation timescale, one may consider that at infinite time, the final bath temperature becomes 186 K with an overall ΔT of 46 K. The rotational ΔT of N₂ at the termination of the trajectory at 3 ns in the current 140 K simulation is 17 K. Therefore, at 3 ns, the simulation is about 37% of the way to a hypothetical complete equilibration. On the other hand, the N₂ rotational ΔT is obtained as 85 K from the experiment. If it is assumed that the experiment is completely re-equilibrated with respect to the translational/rotational modes of N₂, then the experimental ΔT is twice as large as in the simulation. The observed differences in the total magnitude of the N₂ temperature rise and IET efficiency from the simulations and the experiments could have arisen from several sources. There are several potential sources of experimental error (discussed in more detail in the supplementary material), which cannot individually account for all of the differences including the error in the initial [Bz], the [Bz^*] generated at 140 K, 193 nm multiphoton absorption of Bz, and possible enhanced absorption of 193 nm light by Bz dimers. More likely experimental sources of error that may explain the difference between the experiment and simulations are whether NO behaves sufficiently like the bath N₂ molecules at low temperatures such that NO can still act as a proxy for determining N₂ IET properties or if it is possible for NO to quench electronically excited Bz more efficiently at low temperatures. The possible sources of error in the simulation (discussed in more detail in the supplementary material) could be intermolecular potential energy parameters and trajectory initial conditions. However, any errors from the simulation method utilized are not likely to have caused the magnitude of difference between simulated and experimental results.

Classical trajectory calculations at 140 K were performed with one vibrationally excited benzene molecule having an excitation energy of 148.1 kcal/mol undergoing collisional re-equilibration in a benzene/N₂ bath in order to determine the time-dependent internal energy partitioning between internal modes of the Bz^* molecule, the Bz bath molecules, and the N₂ bath molecules. Two different baths, one with 190 N₂ + 9 Bz and another with 396 N₂ + 3 Bz, were considered in these simulation studies. These simulations were performed since previously only two (RT) modes of the bath N₂ molecules had been examined experimentally by proxy via NO LIF temperature measurements, and the experiment did not directly give information about the other molecular modes during the re-equilibration of Bz^*.⁴⁶ Additionally, the 140 K bath experiment and simulations were compared to 300 K experiments and simulations in order to examine the temperature dependence of the IET processes.^6,7

At 140 K, it was found that Bz^*–Bz V-V IET was highly efficient, which is consistent with previous temperature-dependent Bz^*–Bz classical trajectory calculations.⁵⁸ It was previously found that longer lived, multiple-IET-event, “chattering” Bz^*–Bz complexes were formed, which increase the IET efficiency as the lifetime of the complexes increased at lower temperatures. For the bath with 396 N₂ and 3 Bz molecules, the average energy transferred per collision ⟨ΔE_c⟩ was found to be ∼12 kcal/mol for the Bz^*–Bz collisions at the ⟨E(t)⟩_Bz^* of 150 kcal/mol and the bath temperature of 140 K. This value is only about 5 kcal/mol at 300 K. In the present simulations, benzene dimer and trimer formation was quantified at 140 K in order to estimate the effect of chattering collisions on IET efficiency. The partitioning of the energy of R/T modes in these simulations was also consistent with the work of Oref et al., which showed minimal V-R/T IET in Bz^*–Bz collisions.

In Bz^*–N₂ IET, V-R/T IET dominated with negligible V-V IET observed during the duration of the simulations. The Bz^*–N₂ V-R/T IET per collision was lower in magnitude at 140 K than at 300 K, consistent with an impulsive collision model. For the bath with 396 N₂ and 3 Bz molecules, the ⟨ΔE_c⟩ was found to be ∼0.050 kcal/mol for the Bz^*–N₂ collisions at the ⟨E(t)⟩_Bz^* of 150 kcal/mol and the bath temperature of 140 K. This value is about 0.064 kcal/mol at 300 K. The overall ⟨ΔE_c⟩ resulted from all the collisions was found to be about two times higher at 140 K than at 300 K. However, the experimental value is about ten times larger than the result obtained from the simulation. A few of the most likely possible reasons for that could be a higher Bz^* concentration in the experiment than in the simulation and/or the dipole–quadrupole interaction of Bz^*–NO at low temperatures. The latter may lead to a remarkable enhancement of energy transfer to the bath. Therefore, the overall increased rethermalization of Bz^* observed at 140 K in the simulations was due to the increase in Bz^*–Bz V-V IET and not due to other IET processes for these percentages of excited and unexcited Bz molecules in N₂. Therefore, the timescale of collisional relaxation of Bz^* can be tuned in gas systems by varying the molecular composition of the bath. The ability to tune IET in gas systems could allow for the tuning of the onset of turbulence in the gas around hypersonic vehicles allowing for more efficient hypersonic vehicle design; in addition, the ability to tune IET in gas systems could allow for the control of collisional stabilization vs reaction in combustion systems, which could allow for more efficient and/or cleaner combustion.

See the supplementary material for ⟨ΔE_c⟩_Bz^* vs ⟨E(t)⟩_Bz^*, for the 140 K and 300 K simulations, orientation averaged Bz–Bz and Bz–N₂ intermolecular potential energies vs the respective center-of-mass distances, and an additional detailed discussion of potential reasons for differences between the experiment and simulations at 140 K. The factors controlling the IET dynamics at low temperatures as observed in the simulation are also discussed.

The authors would like to thank Bill Hase for all of his hard work and great insights into this project, and they regret that he could not see the final published paper since he passed away during the writing of this paper. Funding for this project was provided by the United States Air Force Office of Research (No. C13-0027), the Air Force Office of Scientific Research under AFOSR Award Nos. FA9550-12-1-0443, FA9550-16-0133, and FA9550-17-1-0107, and the Robert A. Welch Foundation under Grant No. D-0005. Support was also provided by the High-Performance Computing Center (HPCC) at Texas Tech University, under the direction of Philip W. Smith. Parts of the computations were also performed on Robinson, a general computer cluster of the Department of Chemistry and Biochemistry, Texas Tech University, purchased by NSFCRIF-MU under Grant No. CHE-0840493. A.K.P. also acknowledges SERB-DST under File No. ECR/2017/001434 and CSIR under File No. 01 (2998)/19/EMR-II. S.S.A. thanks NIT Meghalaya and MHRD for his fellowship. Parts of this paper were adapted from the following thesis: West, N. Characterization of Benzene Laser-Induced Nonthermal Equilibrium via Nitric Oxide Laser Induced Fluorescence Temperature Measurements, Texas A&M University, College Station, Texas, 2018, with permission from N.A.W.

Some of the data that support the findings of this study are available within this article and its supplementary material in terms of plots. However, all the numerical data are available from the corresponding author upon reasonable request.

The negative of the decrease in the energy of the vibrationally excited benzene molecule, Bz^*, vs time is equal to the sum of the increase in the average energies of the benzene and nitrogen molecules in the bath vs time, i.e.,

The total collision frequency for Bz^* is the sum of its collision frequencies with the Bz and N₂ molecules in the bath, i.e.,

The average energy transfer from Bz^* per collision, ⟨ΔE_c⟩_Bz^*, is given by the left-hand side of Eq. (A3) and may be written as the properly weighted sum of the average energy transfers per collision to Bz and N₂ molecules in the bath as given by Eq. (A5), i.e.,

The average energy transfer per collision to a Bz molecule in the bath is then