Disorder-induced heating in molecular atmospheric pressure plasmas

Recent work has shown that ions are strongly coupled in atmospheric pressure plasmas when the ionization fraction is sufficiently large, leading to a temperature increase from disorder-induced heating (DIH) that is not accounted for in standard modelling techniques. Here, we extend this study to molecular plasmas. A main finding is that the energy gained by ions in DIH gets spread over both translational and rotational degrees of freedom on a nanosecond timescale, causing the final ion and neutral gas temperatures to be lower in the molecular case than in the atomic case. A model is developed for the equilibrium temperature that agrees well with molecular dynamics simulations. The model and simulations are also applied to pressures up to ten atmospheres. We conclude that DIH is a significant and predictable phenomena in molecular atmospheric pressure plasmas.


I. INTRODUCTION
Cold atmospheric pressure plasmas (CAPPs) are continuously being tested and developed for a wide range of applications, including wound healing and cancer treatment 1 , water purification 2 , and food preservation 3 , among others [4][5][6] .This comes, in large part, from their unique capacity for generating large quantities of reactive oxygen and nitrogen species (RONS).Additionally, in comparison with other plasma technologies, CAPPs offer a lower cost, complexity, and carbon footprint since they do not require the use of a vacuum chamber 7 .This has driven increased interest in understanding plasma dynamics, especially with concern to the transport of RONS 7 .A significant amount of work has already been done to characterize and model these plasmas [8][9][10][11] .It is, however, also important to understand the underlying physical mechanisms that drive CAPPs to be sure they are accounted for in the models.In this paper, we examine and extend one such mechanism: disorder-induced heating (DIH), which arises from strong ion-ion coupling.
Recent work showed that ions in atmospheric pressure plasmas are strongly coupled at ionization fractions as low as 10 −5 , leading to rapid ion DIH on a picosecond timescale following an instantaneous ionization 12 .This work modelled the effect in monatomic argon plasma and found that, at high ionization fractions (≳ 10 −2 ), it can increase ion and neutral gas temperatures by thousands of Kelvin.In this paper, we extend the previous work to molecular plasmas and examine the effects of added rotational degrees of freedom.Since many applications of CAPPs rely on N 2 and O 2 [13][14][15] , it is important to understand how disorder-induced heating impacts molecular plasmas.We also extend the work to higher pressures (up to 10 atm) to account for the regimes seen in plasmaassisted combustion 16 .
Using molecular dynamics (MD) simulations of Nitrogen gas (N 2 ), we find that DIH causes the ion translational temperature to rapidly increase over approxi-mately one ion plasma period (≈ 1 ps).This is followed by translational-rotational and ion-neutral relaxation processes that occur in roughly 1 ns.The latter process leads to fast neutral gas heating.A main finding is that the translational-rotational energy exchange leads to a lower total temperature resulting from DIH than in the atomic case.This is because the kinetic energy gained by ions from DIH gets distributed amongst more degrees of freedom in the molecular case.Considering N 2 , the kinetic energy gain is eventually spread equally over 3 translational and 2 rotational degrees of freedom, for a total of 5.This contrasts with the 3 degrees of freedom in the atomic case.Furthermore, we confirm that this principle extends to molecules with 3 rotational degrees of freedom by simulating an N 4 molecular plasma.
Despite the reduction in temperature due to the presence of rotational degrees of freedom, DIH can still increase the neutral gas temperature by thousands of Kelvin.We derive a model for this change in temperature using conservation of energy arguments and find that it matches well with our MD results.This model is then generalized to account for higher pressures, where DIH becomes more significant.Finally, we consider the common case of a gradual ionization ramp over several nanoseconds.As in a previous study with atomic gases, the temperature increase of ions due to DIH and the subsequent ion-neutral energy relaxation are fast compared to a typical ionization timescale in a nanosecond discharge (∼ 10 ns).A consequence is that both the ion and neutral temperatures increase at the rate of ionization.Importantly, it is also shown that translationalrotational relaxation is fast compared to the ionization timescale -thus, in the context of fast gas heating, one must consider the effect of molecular rotation.
Fast gas heating is a particularly important implication of DIH, with applications in flow control, combustion, and CO 2 disassociation [17][18][19] .There are models for this heating based on the relaxation of excited states in atoms and molecules 20 .These models, however, do not consider DIH and as a consequence, may underestimate heating when the ion density is sufficiently high.
While many atmospheric pressure plasmas do not reach high enough levels of ionization for DIH to be significant, those generated from nanosecond discharges, arc discharges, and laser ionization often do [21][22][23][24][25][26][27] .In these instances, we expect DIH to be an important source of fast gas heating.In laser-produced plasmas, DIH is the fastest known heating mechanism present and therefore sets the ion temperature at early times 25 .In nanosecond discharges, neutral gas heating predicted by DIH has shown very good agreement with time-resolved gas measurements from experimental work 28 , giving further merit to the importance of this process.

II. STRONG COUPLING AND DISORDER-INDUCED HEATING
An interaction between species s and s ′ is strongly coupled if the average potential energy is on the same order as, or greater than, the average kinetic energy.This can be quantified by the coupling parameter exceeding unity, Γ ss ′ ≥ 1, where Here, ϕ ss ′ is the interaction potential between s and s ′ , a ss ′ is the average inter-particle spacing, and T is the temperature of the two species.For particles of the same species, one can write a ss = (3/4πn s ) 1/3 , which corresponds to an average inter-particle spacing ∼ 1 nm at atmospheric pressure.Strongly coupled plasmas behave fundamentally differently than weakly coupled (ideal) plasmas.Binary collisions characteristic of weak coupling are replaced by many-body interactions that force the strongly coupled species into an ordered spatial configuration.Consequently, when a gas is ionized to the point of strong ionion coupling, the newly created ions move from a highly disordered to a more ordered state.This increased order results in a lower potential energy which, in turn, is offset by an increase in kinetic energy.This increase in kinetic energy is termed disorder-induced heating 29,30 .The process is illustrated in figure 1.It can be seen that the separation of ions leads to a spatial configuration with more apparent order than the initial configuration, implying a transfer of potential to kinetic energy.
Disorder-induced heating is not a new phenomenon.It has been studied in great detail in ultracold neutral plasmas [30][31][32][33] and more recently, fusion plasmas 34 .However, it has now been shown that DIH is also expected to occur in atmospheric pressure plasmas 12 .While this has not been confirmed experimentally, predictions from DIH have been shown to match well with experimental data 12,28 .

III. MD SIMULATIONS
Molecular dynamics simulations were run using the open-source software LAMMPS 35 .
Nitrogen (N 2 ) molecules were modelled as rigid rotors, allowing for no vibration, and all ions were modelled as N + 2 , consisting of a neutral atom rigidly bonded to an ion.A bond length of 1.09 Å36 was held constant using the "fix rigid" command in LAMMPS, which solves the forces and torques on each rigid body by summing the forces and torques on its constituent atoms.This approach ignores vibration, an assumption expected to be valid because the characteristic vibrational temperature of N 2 is 3374 K 37 and the temperature is far below this for the majority of the simulation.Additionally, DIH and the subsequent ion-neutral energy relaxation occur at the nanosecond timescale or less, while translational-vibrational energy exchange in N 2 occurs on a microsecond timescale 38 .A real ionization process will certainly create excited vibrational states, but this does not affect DIH because the ion and neutral atom temperatures are too low to couple to the vibrational modes.
Electrons were left out of the simulation because their high temperature in CAPPs make them weakly coupled.In strongly coupled plasmas, treating electrons as a non-interacting background species yields accurate ion and neutral dynamics, as commonly done in the onecomponent plasma (OCP) model 39 .Chemical effects were also left out as they are beyond the scope of this work.We aim not to simulate a CAPP as a whole, but instead examine a specific heating mechanism.Since DIH is caused solely by a change in the spatial orientation of ions, chemistry only affects the heating if it causes a significant change in the number density or temperature of ion species in the plasma on a relevant timescale.A model is provided (equation ( 7)) that can account for those changes.
Interactions were modelled on a per-atom basis, with the ion-ion, ion-neutral, and neutral-neutral interaction potentials given by the Coulomb, charge-induced dipole, and Lennard-Jones potentials as follows where ϵ = 99.8kB , σ = 0.3667 nm, and α R = 7.5 for N 2 and a 0 is the Bohr radius.The repulsive term in the charge-induced dipole potential was added to avoid bound states, as done in previous work 12 .Details on this addition can be found in the Appendix of Ref. 12.Much of this paper will focus on data from the N 2 simulations.However, for comparison's sake, simulations of N and N 4 were also run.The atomic nitrogen simulations were run in the exact same fashion as the previous work 12 .In the N 4 simulations, N 4 was modelled using an open-chain geometry as described by Glukhovtsev and Laiter 40 ; a structure which, importantly, has three rotational degrees of freedom.All N 4 ions were modelled as N + 4 , with one of the two middle atoms ionized and a modified LJ coefficient ϵ = 29.0kB was used in order to prevent significant neutral gas correlations.
To start the simulations, 10,000 molecules of neutral gas were evolved under the influence of a thermostat until equilibrium at room temperature was reached.Then, a fraction of the molecules were instantly ionized and an NVE (microcanonical) simulation was run with timestep ∆t = 5 × 10 −4 ω −1 pi , a value found to be sufficiently small for conserving energy.It has been shown in previous work 28 and is confirmed in this work (section VII) that the magnitude of heating from DIH is independent of ionization dynamics.
Molecule positions and velocities were output once every plasma period and used to calculate translational and rotational temperatures.The simulation domain was a 3-dimensional box with periodic boundary conditions.Three simulations for each of N, N 2 , and N 4 were run at every tenth ionization fraction from 0.10 to 0.70.Data was averaged across the three trials to minimize statistical deviation.

IV. MD RESULTS
The averaged time-evolution of the component temperatures from a particular set of simulations, N 2 at a 0.30 ionization fraction, can be seen in figure 2. This plot shows results of the NVE stage simulation and t = 0 refers to the time immediately following ionization.A sharp spike can be seen in the translational temperature of N + 2 , reaching its peak after 2.22 ps (1.52 ω −1 pi ).This is disorder-induced heating.Molecules, upon being ionized, structurally rearrange themselves until they reach their state of minimal potential energy, increasing their translational kinetic energy in the process.This heating is followed by collisional relaxation, during which the N + 2 translational temperature equilibrates with the rest of the degrees of freedom in the simulation over the course of one nanosecond, ultimately reaching an equilibrated temperature of 732 K.For a clear physical picture of this process, recall figure 1.
Simulations of each gas type, N, N 2 , and N 4 , followed a similar time-evolution to figure 2. Ion translational temperatures always peaked after 1.5 plasma periods, consistent with previous work 12 , and equilibration always occurred in roughly one nanosecond.However, they all produced differing equilibrium temperatures.We found the equilibrium temperature scales with the ionization fraction x 4/3 i , again consistent with the previous work 12 .This is caused by an increased ion density and thus an increased chance of two ions being generated in close enough proximity to gain an energy exceeding the initial thermal energy as they repel.
It is also observed that the equilibrated temperature is inversely proportional to the number of degrees of freedom present, implying T eq,N > T eq,N2 > T eq,N4 , as shown in figure 3.This comes as a direct consequence of the equipartition theorem.Figure 4 shows that, for a given ion density, N, N 2 , and N 4 plasmas convert the same amount of potential energy to kinetic energy.However, in N 4 , this energy must be spread across 6 degrees of freedom (3 translational and 3 rotational), while it is spread across just 5 in N 2 and 3 in N. It is clear, then, that DIH models must consider the number of active degrees of freedom.

V. MODELLING DIH A. Precise Formulation
Disorder-induced heating results from a change in the spatial ordering of ions.The ion-ion radial distribution function, which outputs the ratio of the local ion density FIG. 3. Simulated and predicted temperatures of N, N2, and N4 plasma using the precise formulation from equation (8)  with gii(r) computed from MD, and from the approximate model of equation (10).
around each ion to the bulk ion density 41 can capture this change.Figure 5 shows an example of g ii (r) computed from MD simulation data using equation (5) at various times throughout DIH.Immediately following ionization, g ii (r) ≈ 1, indicating a disordered spatial configuration where each ion's position is effectively random.Over time, a gap forms at small radial distances, indicating the ions have spread apart from one another.This spreading apart, also depicted in figure 1, affects the potential energy of ions.
The radial distribution function provides a convenient means to quantify the potential energy.Integrating radially and multiplying by the ion density, ρ i , gives the density of ions at each distance r + dr away from each reference ion.Therefore, the potential energy of ions in terms of the ion-ion distribution function is given by 42 where N i is the total number of ions.For modelling the effect of disorder-induced heating on the equilibrium temperature, we are only concerned with how the potential energy of ions changes from immediately following ionization to equilibrium.Additionally, to avoid simulation artifacts like the number of ions, N i , we are interested only in an energy density.Therefore, we convert the previous equation to a change in energy density as follows where x i is the ionization fraction, ϕ ii (r) is the Coulomb potential, and g ii (r, t eq ) and g ii (r, t 0 ) are the ion-ion radial distribution functions at equilibrium and immediately following ionization respectively.A change in the spatial orientation of ions also affects the ion-neutral interaction, which can be quantified through the ion-neutral radial distribution function (g in (r)).However, we found the contribution to the potential energy change from ion-neutral interactions to be much less than from ion-ion interactions, and thus deemed the ion-neutral interaction to have a negligible effect on the equilibrium temperature.It is worth noting that at sufficiently high pressures, this assumption is unlikely to be valid because the ion-neutral interaction becomes strongly coupled, becoming a source of additional heating.However, at the pressures presented in this paper (up to 10 atm), we found this to be a valid approximation.
Hence, we use equation (7) to model the entirety of the change in potential energy from DIH and convert this to a temperature change using traditional thermodynamics arguments where f is the number of degrees of freedom.
Figure 3 shows close agreement between the model and simulated temperatures.Here, the radial distribution functions used in the model were computed from the MD simulation data itself.We observe a slight discrepancy from the simulated temperatures, which can be attributed to noise in the radial distribution function.This noise can be seen in figure 5, and comes as an artifact of the finite number of ions and finite number of timesteps over which we can compute this function.
Equation ( 8) is, by nature, a flexible formulation for the equilibrium temperature.One can consider the effects of recombination by splitting equation ( 7) into having an "initial" and "equilibrium" ionization fraction and ion density.Additionally, it is additive, such that if multiple species of ions exist in a plasma, one need only sum the contributions of each.To use the model, there exist several methods for estimating radial distribution functions including the Percus-Yevick 43 and hypernetted-chain 44 approximations.Future work will use these methods to model DIH and be implemented into a fluid model.This approach relies on numerically evaluating a model for the radial distribution functions.Next, we consider a simpler more approximate approach to estimate the potential energy change in DIH directly.

B. Approximate Model
A simpler estimate for DIH can be constructed by analyzing the change in average ion-ion spacing.In the limiting case of transitioning from a completely disordered to perfectly ordered state, we expect the average ion spacing to change from approximately a ii /2 to a ii , leading to a change in potential energy ∆u max = x i e 2 4πϵ 0 We found that in the ion densities relevant to strongly coupled CAPPs, a linear correction factor of κ = 2/3 is required to fit our data, such that ∆u = κ∆u max .Converting equation ( 9) to a change in temperature, we find Equation ( 10) is a slight modification of the model used in Ref. 12 to account for the added degrees of freedom in molecular motion.
A comparison of the prediction of equation ( 10) and MD simulations is shown in figure 3. The close agreement with the simple model is in part due to the use of the κ = 2/3 fit coefficient, but the model captures well the scaling with ionization fraction.It also highlights the importance of the degrees of freedom in distributing the energy gained by DIH.The change in temperature of N (f = 3), N 2 (f = 5), and N 4 (f = 6) differ only by a factor of 1/f .While it lacks the physical precision associated with the previous model, it can still be useful as a fairly accurate means to estimate the extent to which DIH will affect a given plasma.

VI. HIGHER PRESSURES
Simulations were also run with initial gas pressures of 2.5, 5.0, 7.5, and 10 atm at constant ionization fractions to explore the pressure dependence of DIH.This is particularly applicable to plasma-assisted combustion, where initial gas pressures can range from 0.1-10 atm 16 and fast gas heating is known to play a role in accelerating and intensifying the combustion process 20 .
A higher initial gas pressure leads to a larger ion density and thus, we can expect more heating.MD results, shown in figure 6, confirm this.Additionally, it is shown that the precise formulation from equation ( 8) remains accurate, while the approximate model from equation ( 10) requires an increased value for the correction factor κ to fit MD data.Scaling κ as P 1/20 , where P is the initial gas pressure, gives good agreement.

VII. GRADUAL IONIZATION
All simulations discussed thus far have involved an instantaneous ionization process.However, in atmospheric pressure plasmas, ionization processes often occur over the course of several nanoseconds 21,45 .It has already been shown that, in atomic plasma, the final temperature following disorder-induced heating and subsequent ion-neutral temperature equilibration are unaffected by the rate of ionization 28 .Considering that translationalrotational relaxation occurs on the same timescale as ionneutral relaxation in these plasmas, we expect this to be the case for molecular plasmas as well.
To confirm this, an MD simulation was run with a gradual ionization process.The ion density was increased linearly over the course of 2.2 ns until an ionization fraction of 0.30 was reached.The results, shown in figure 7, confirm expectations.The simulation reaches an equilibrium temperature of 752 K, matching the simulations from an instantaneous ionization process to a 2.7% error.It can also be observed that all degrees of freedom present in the simulation heat at, roughly, the rate of ionization.This is because ion-neutral and translationalrotational relaxation in these plasmas occur on a faster timescale than the ionization process.Finally, it should be noted that early on in the simulation, when the ionization fraction is low, there are not enough ions to have a well-defined temperature.Consequently, the early temperature data has a significant amount of noise.The relaxation dynamics present in both the instantaneous and gradual ionization process are complex.Strongly coupled ion-ion collisions play a large role in the system's translational-rotational energy exchange, so attempts to model the relaxation with the Boltzmann equation 46 and Parker's model 47 (a model for RT relaxation in neutral gas) fail, as they do not account for strong coupling.Thus, this work motivates the need for a theory of translational-rotational relaxation in a strongly coupled plasma.

VIII. CONCLUSIONS
This work examined how plasma and gas heating due to DIH changes with the introduction of molecules and higher pressures.Using MD, it was shown that energy from DIH is spread to rotational degrees of freedom on a nanosecond timescale, leading to a lower equilibrium temperature than in the atomic case.It was also shown that a higher initial gas pressure leads to a larger temperature increase from DIH.Two models were derived that can predict this behavior to a high degree of accuracy.Finally, by simulating a gradual ionization over the course of a few nanoseconds, it was shown that the temperature change due to DIH in molecular plasmas is unaffected by ionization dynamics.
This work shows that the effects of strong coupling remain significant in molecular plasmas and establish that one must consider all active degrees of freedom in order to accurately model the influence of DIH on the overall heating process.It further motivates the need to develop new modelling techniques that can account for strong coupling in CAPPs, as current models based on the Boltzmann equation fail to do this.These models do not consider DIH and as such, they risk vastly underestimating the temperature.

FIG. 1 .
FIG.1.Illustration of the disorder-induced heating process.Ionization transforms neutral gas into a nonequilibrium plasma.Coulomb repulsion causes the ions to separate on a picosecond timescale, causing significant translational heating.Collisions then distribute this energy across neutral and rotational degrees of freedom until they reach an equilibrium.The collision relaxation timescale near atmospheric pressure is often characteristic of nanoseconds.

FIG. 2 .
FIG. 2. Simulated time-evolution of the temperature components of an N2 plasma at a 0.30 ionization fraction.Lines correspond to the ion translational temperature (blue), ion rotational temperature (red), neutral molecule translational temperature (green), and neutral molecule rotational temperature (yellow).

FIG. 4 .
FIG.4.Change in total kinetic energy density (with respect to t = 0) for N, N2, and N4.This data comes from simulations at an ionization fraction of 0.30.

FIG. 5 .
FIG. 5. Snapshots of the ion-ion radial distribution function over the course of DIH from a simulation of N2 at a 0.30 ionization fraction.Over time, the ion-ion radial distribution function shifts to the right as ions spread apart from each other, reducing their potential energy in the process.

FIG. 6 .
FIG.6.Simulated and predicted temperatures of N2 as a function of the initial gas pressure.The approximate model from equation (10) must be scaled by P 1/20 to fit MD data.

FIG. 7 .
FIG.7.Simulated time-evolution of the temperature components of an N2 plasma subject to a gradual ionization process (2.2 ns) until a 0.30 ionization fraction is reached.Initial ion temperatures are made a lighter color because early in the simulation, there are not enough ions to have a well-defined temperature.