]Theory and simulation of shock waves freely propagating through monoatomic non-Boltzmann gas

[1]Malte<EMAIL_ADDRESS>
*[1]Chair of High Pressure Gas Dynamics, Shock Wave Laboratory, RWTH Aachen University, Aachen, 52056, Germany

The effect of non-Boltzmann energy distributions on the free propagation of shock waves through a monoatomic gas is investigated via theory and simulation.
First, the non-Boltzmann heat capacity ratio $\gamma$, as a key property for describing shock waves, is derived from first principles via microcanonical integration.
Second, atomistic molecular dynamics simulations resembling a shock tube setup are used to test the theory.
The presented theory provides heat capacity ratios ranging from the well-known $\gamma = 5/3$ for Boltzmann energy-distributed gas to $\gamma \to 1$ for delta energy-distributed gas.
The molecular dynamics simulations of Boltzmann and non-Boltzmann driven gases suggest that the shock wave propagates about 9 % slower through the non-Boltzmann driven gas, while the contact wave appears to be about 4 % faster if it trails non-Boltzmann driven gas.
The observed slowdown of the shock wave through applying a non-Boltzmann energy distribution was found to be consistent with the classical shock wave equations when applying the non-Boltzmann heat capacity ratio.
These fundamental findings provide novel insights into the behavior of non-Boltzmann gases and might help to improve the understanding of gas dynamical phenomena.

§ INTRODUCTION

Non-equilibrium gas dynamics are essential for describing atmospheric entry processes [1, 2], hypersonic cruise [3], and nozzle flow [4, 5].

The mathematical formulation of such processes is well-established through the Boltzmann equation [6, 7], which allows for a continuous description of, typically, rarefied gas dynamics.
Non-Boltzmann dynamics of dense gases are commonly modeled through approximate schemes, many of which use the assumption of microcanonical equilibrium.
The present work provides a microcanonical non-equilibrium derivation of the heat capacity ratio of non-Boltzmann energy-distributed monoatomic gas.
Atomistic molecular dynamics simulations of shock waves freely propagating through Boltzmann and non-Boltzmann energy-distributed monoatomic gases are utilized to validate the theoretical findings.
Ultimately, the simulation results are compared to ideal equilibrium shock calculations using non-equilibrium heat capacity ratios.

The non-equilibrium state of a monoatomic gas is characterized by variations of the number density over the entire accessible phase space, spanned by the positions $q$ and the momenta $p$.
The Boltzmann equation provides the means for continuously describing the evolution of such a non-equilibrium state by phase space integration [7].
Solving the Boltzmann equation, however, becomes increasingly demanding with decreasing Knudsen number, hence it is typically applied only for rarefied gas flows [6].
Probably the most popular approach for averting the high computational efforts required for solving the Boltzmann equation is the use of Direct Simulation Monte Carlo (DSMC), which was pioneered by Bird [8].
For a good overview of relevant studies about and with the DSMC method, please refer to Gallies [9].

For the transition regime from large to small Knudsen numbers, usually more efficient, yet more approximate approaches are utilized [10].
For intermediate Knudsen numbers the use of moment equations has been established, for which Torrilon [10] provided a good overview of how moments are defined and how these equations can be used for modeling non-equilibrium gas flows.
For the strong non-equilibrium start-up phase of a shock tube process for instance, Au et al. [11] found that the use of moment equations allows for reliably describing the first moments of start-up, which cannot be captured by the Navier-Stokes equation.
Another example is provided in the recent special issue on "nonequilibrium thermodynamics" edited by Ván [12], in which Öttinger et al. [13] formulated moment equations for rarefied gases based on the Boltzmann equation, showing that the rotational extended thermodynamics (RET) approach [14] is a special case of the general equation for the non-equilibrium reversible-irreversible coupling (GENERIC) [15], with both approaches being used for describing non-equilibrium gas dynamics.

With increasing gas density, i.e. decreasing Knudsen number, the models used for describing non-equilibrum start using the assumption of Boltzmann energy-distributed gases.
Various approaches for describing non-equilibrium thermodynamics have been developed [16], such as superstatistics [17], stochastic thermodynamics [18], and local equilibrium thermodynamics [19].
The above methodologies all aim at formulating an unified description of non-equilibrium thermodynamics [12], yet they do not necessarily include non-Boltzmann energy distributions, i.e. local non-equilibrium in the momentum space.
The most widely used approach for modeling microcanonical non-equilibrium gas states is the use of distinct temperatures for distinguishable degrees of freedom [20].
For instance, Kosareva et al. [2] recently applied a four-temperatures kinetic model for modeling the vibrational relaxation of CO2 during atmospheric entry.
This approach can generally be used to separate excitation of translational, rotational, vibrational, and electronic degrees of freedom, yet it is still based on the concept of the Boltzmann energy distribution.

Besides the research on the Boltzmann equation and its derivatives for describing non-Boltzmann energy distributions, the role of non-Boltzmann energy distributions has been intensely investigated in the field of theoretical chemical kinetics recently [21, 22].
In 2015, Klippenstein and co-workers [21, 22] revived the topic of non-Boltzmann chemical kinetics, which had probably been initially observed in the 1970s and 1980s [23, 24, 25].
Subsequent studies by Labbe et al. [26, 27] and by Döntgen et al. [28, 29] advanced the fundamental understanding of non-Boltzmann chemical kinetics.
In particular, the prompt dissociation reactions of radical species was found to significantly affect macroscopic properties, such as the flame speed [27] and the ignition delay time [30].
The combination of non-Boltzmann gas dynamics and non-Boltzmann chemical kinetics was intensely investigated by Kustova and co-workers [31, 2, 32, 33], who investigated not only the role of vibrational exciation on CO2 chemistry during atmospheric entry [2], but also the behavior of N2 [31], O2 [32], and CO [33] in vibrationally excited gas dynamical flows.
In this context, the use of state-to-state dynamics is widely established for describing the transition between energy levels and chemical states [31, 34, 35, 2, 32, 33].

Non-equilibrium effects are known to or might be affecting various gas dynamical experiments.
For instance, Campbell et al. [36] compared calculated and experimentally derived thermodynamic variables for postshock states of vibrationally excited, i.e. non-Boltzmann, gases.
The authors highlighted that vibrational excitation can lead to deviations in temperatures of up to 8 % [36].
As a less definite example, Bedin [37] found an anomalous alteration of the velocity of a shock wave traversing a non-equilibrium plasma, which was intensely re-investigated experimentally [38, 39] and numerically [40, 41].
The controversial discussion of the underlying mechanism is not yet fully resolved [42], indicating the high complexity of the underlying mechanism and non-equilibrium gas dynamics in general.

The present work will show that non-Boltzmann energy distributions affect the heat capacity ratio, as the most relevant gas dynamical property, even at high density conditions.
First, the heat capacity ratio of a non-Boltzmann energy-distributed monoatomic gas will be derived from first principles by applying and integrating a microcanonical formulation of the heat capacity ratio, weighted by prototypical energy distributions.
Second, the theoretical dependence of the heat capacity ratio on the energy distribution will be validated through atomistic molecular dynamics simulations of a monoatomic gas for the two presently limiting energy distributions: The Boltzmann and the delta energy distributions.
Third, the shock velocities obtained through the molecular dynamics simulations will be compared to the ideal equilibrium shock calculations, carried out with the theoretical heat capacity ratios of the Boltzmann and delta energy distributions.
The present work will provide fundamental insights into the role of non-Boltzmann energy distributions on the most fundamental gas dynamical property: The heat capacity ratio.

§ THEORY

The microcanonical heat capacity ratio, i.e. isentropic coefficient, $\gamma$ of a monoatomic ideal gas with an arbitrary energy distribution is derived here.
In the present work, the Boltzmann energy distribution (B) will be used as reference of systems in thermal equilibrium and a delta function ($\delta$) will be used to represent the upper limiting case of non-Boltzmann energy distributions.
In essence, a delta energy distribution represents a system in which all particles have exactly the same energy.
Such an energy distribution has limited physical relevance, since any collisional interaction would disrupt the delta energy distribution.
However, it allows to quantify the upper limiting effect of non-Boltzmann energy distributions on the heat capacity ratio.

The general definition of the heat capacity ratio is given as follows [43]:

\begin{equation}\label{eq:gamma_heatcapacities}
\gamma = \frac{C_{P,j}}{C_{V,j}} ,
\end{equation}

with $j$ identifying the energy distribution, $C_{P,j}$ being the the isobaric heat capacity of distribution $j$, and $C_{V,j}$ being the isochoric heat capacity of distribution $j$.
The heat capacity is typically defined as the second derivative of the Helmholtz free energy with respect to temperature.
Here, however, it is more suitable to describe the heat capacity through the energy fluctuations of the system.
In a Boltzmann energy-distributed system, the heat capacity is given as:

\begin{equation}\label{eq:heatcapacity}
C_{i,j} = k_\text{B} \cdot \frac{\langle E^2 \rangle_{i,j} - \langle E \rangle_{i,j}^2}{(k_\text{B} T)^2} ,
\end{equation}

with $i$ being $P$ or $V$, $k_\text{B}$ being Boltzmann's constant, $E$ being the energy, and $\langle ... \rangle$ being the microcanonical averaging operator.
When combining equations <ref> and <ref>, one obtains the following expression:

\begin{equation}\label{eq:gamma_energies}
\gamma = \frac{\langle E^2 \rangle_{P,j} - \langle E \rangle_{P,j}^2}{\langle E^2 \rangle_{V,j} ,
- \langle E \rangle_{V,j}^2}
\end{equation}

Since equation <ref> solely depends on the averaged squared and squared averaged energies, $\langle E^2 \rangle$ and $\langle E \rangle^2$, respectively, it is applicable for arbitrary energy distributions.
The averaged squared energy $\langle E^2 \rangle$ is given as:

\begin{equation}\label{eq:ave_E2}
\langle E^2 \rangle_{i,j} = \frac{\int\limits_0^\infty E^2 \cdot \rho_i(E) \cdot f_j(E) \cdot \text{d}E}{\int\limits_0^\infty \rho_i(E) \cdot f_j(E) \cdot \text{d}E} ,
\end{equation}

with $\rho_i(E)$ being either the isobaric density of states ($i = P$) or the isochoric density of states ($i = V$) and $f_j$ being the propabilities of energy distribution $j$.
The squared averged energy $\langle E \rangle^2$ is given as:

\begin{equation}\label{eq:ave2_E}
\langle E \rangle_{i,j}^2 = \left( \frac{\int\limits_0^\infty E \cdot \rho_i(E) \cdot f_j(E) \cdot \text{d}E}{\int\limits_0^\infty \rho_i(E) \cdot f_j(E) \cdot \text{d}E} \right)^2 .
\end{equation}

For a monoatomic gas, the density of states is obtained through inverse Laplace transformation of the translational partition function [44] and is independent of the energy distibution.
In case an isobaric state is considered, the density of states is given as:

\begin{equation}\label{eq:rho_P}
\rho_P(E) = \frac{q_\text{tr} / p}{\Gamma(5/2)} \cdot (E)^{3/2} ,
\end{equation}

with the translational partition function pre-factor $q_\text{tr} = (2 \pi m / h^2)^{3/2}$, the pressure $p$, and the gamma function $\Gamma$.
In case an isochoric state is considered, the density of states is given gas:

\begin{equation}\label{eq:rho_V}
\rho_V(E) = \frac{q_\text{tr} \cdot V}{\Gamma(3/2)} \cdot (E)^{1/2} ,
\end{equation}

with the volume $V$.

Combining equations <ref> to <ref> yields the heat capacity ratio $\gamma$ for a given energy distribution $f_j(E)$.
For the Boltzmann energy distribution $f_B(E) = \exp(-E / k_\text{B} T)$ and the delta energy distribution $f_\delta(E) = \delta(E - E_0)$, $\gamma$ can be obtained analytically and amounts to $\gamma = 5/3$ and $\gamma \to 1$, respectively.
To allow for comparability between different energy distributions, the total energy of the compared energy distributions is required to be equal.
For the delta and Boltzmann energy distributions, the following equation relates the delta peak position $E_0$ to the temperature $T$ of the Boltzmann energy distribution.

\begin{equation}
\begin{split}
\langle E_\delta \rangle &\overset{!}{=} \langle E_\text{B} \rangle \\
\Leftrightarrow \frac{\int\limits_0^\infty E \cdot \rho (E) \cdot f_\delta(E) \cdot \text{d}E}{\int\limits_0^\infty \rho (E) \cdot f_\delta(E) \cdot \text{d}E} &= \frac{\int\limits_0^\infty E \cdot \rho (E) \cdot f_\text{B}(E) \cdot \text{d}E}{\int\limits_0^\infty \rho (E) \cdot f_\text{B}(E) \cdot \text{d}E} \\
\Leftrightarrow \frac{\int\limits_0^\infty E \cdot \sqrt{E} \cdot \delta(E - E_0) \cdot \text{d}E}{\int\limits_0^\infty \sqrt{E} \cdot \delta(E - E_0) \cdot \text{d}E} &= \frac{\int\limits_0^\infty E \cdot \sqrt{E} \cdot \exp\left(-\frac{E}{k_\text{B} T}\right) \cdot \text{d}E}{\int\limits_0^\infty \sqrt{E} \cdot \exp\left(-\frac{E}{k_\text{B} T}\right) \cdot \text{d}E} \\
\Leftrightarrow E_0 &= 3/2 \cdot k_\text{B} T
\end{split}
\end{equation}

The transition between the Boltzmann and delta energy distributions is modeled through the Gaussian energy distribution $f_\text{G}(E) = \frac{1}{\sqrt{2 \pi} \cdot b} \cdot \exp\left(-\frac{(E - E_{\text{G},0})^2}{2 b^2}\right)$ with the Gaussian width $b$ and the Gaussian position $E_{\text{G},0}$.
When setting $\langle E_\text{G} \rangle \overset{!}{=} \langle E_\text{B} \rangle$, the Gaussian width $b$ can be determined numerically for a given Gaussian position $E_{\text{G},0}$.
Through this procedure, the Gaussian distribution can be converged to the Boltzmann distribution for $E_{\text{G},0} \to \text{-}\infty$ and to the delta distribution for $E_{\text{G},0} \to 3/2 \cdot k_\text{B} T$.
For the Gaussian energy distributions, however, $\gamma$ cannot be determined analytically and the integrals in equations <ref> and <ref> have to be solved numerically.

§ METHODOLOGY

Pseudo three-dimensional Lennard-Jones (LJ) molecular dynamics (MD) simulations have been carried out with super-critical Argon in the driver and driven sections at an initial temperature of 300 K and ideal gas pressures of 500 and 50 bar, respectively, using the LAMMPS software package [45].
The Lennard-Jones potential used in the present MD simulations is defined as follows [46]:

\begin{equation}
V_\text{LJ} (r) = 4 \cdot \varepsilon \cdot \left ( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right) ,
\end{equation}

with distance $r$ between two LJ particles, LJ well depth $\varepsilon$, and LJ diameter $\sigma$.
The large pressures in the driver and driven sections are selected to generate a sufficiently strong shock wave for shock position post-processing and to keep the width of the shock wave rather small.
Thermodynamic real gas effects will be considered in the theoretical shock velocity calculations, as detailed below.
Argon is modeled via the LJ well depth $\varepsilon = 114\;\text{K}$ and the LJ diameter $\sigma_{3\text{D}} = 3.47\;\text{\AA}$.
To reduce the simulation size, the three-dimensional domain has been projected onto a two-dimensional domain and the LJ diameter has been adopted by requiring that the ideal gas LJ collision frequencies [47] of the three- and two-dimensional domains are equal, yielding $\sigma_{2\text{D}} = 2 \cdot \sigma_{3\text{D}}^2 / L$, with $L$ being the reference length of the reduced dimension.
Here, $L = 5\;\text{\AA}$ has been used.

The two-dimensional simulation box was $1\;\text{\textmu}\text{m}$ long and $0.2\;\text{\textmu}\text{m}$ wide.
The area defined by the first $0.1\;\text{\textmu}\text{m}$ of the long axis was filled with 11941 driver gas molecules, the remaining area was filled with 10852 driven gas molecules.
All molecules have been randomly positioned inside the respective areas.
Both ends of the short axis were modeled with periodic boundary conditions, through which the corresponding dimension of the simulation was pseudo-infinite.
The boundaries at both ends of the long axis were modeled as perfectly reflecting walls, through which collisions of gas particles with these walls were perfectly elastic.
The driver and driven sections were initially separated by five layers of particles with frozen positions.
Before removing the separating particle layers, the particle velocities were randomly distributed so that they resemble a Maxwell-Boltzmann distribution and each section was thermalized at 300 K for 0.1 ns using a Nosé-Hoover thermostat.
Then, the separating particle layers were removed and the entire system was simulated for 2 ns without thermostating.
Post-processing, however, was stopped before the shock wave reached the endwall of the driven section.
A total of five repetitions with different random number generator seeds have been carried out for the Boltzmann and the delta energy-distributed driven gases.

The propagation of the contact wave between the driver and driven gases and the shock wave inside the driven gas was tracked through the particle density over the long axis of the simulation box.
The driven gas density was fitted using two logistic functions, one for modeling the contact wave and the other for modeling the incident shock wave.
The positions of the transitions from the lower to the upper limits of the logistic functions are used as contact and shock wave positions in the post-processing.

Simulation of the Boltzmann energy-distributed driven gas was immediately possible, since the driven gas had been thermalized through the Nosé-Hoover thermostat and the subsequent constant particle (N), constant volume (V), and constant energy (E) simulations preserved the Boltzmann energy distribution in the driven gas.
A delta energy-distributed gas, however, would rapidly thermalize towards a Boltzmann energy-distributed gas at the present high-pressure conditions, which come with large collision frequencies.
Any collision within the driven gas would re-distribute the kinetic energy and relax the energy distribution towards the Boltzmann distribution.
This relaxation is especially pronounced for the delta energy distribution, as can be concluded from the H-function presented in Figure <ref>.
In order to limit the time available for relaxation of the delta energy distribution prior to arrival of the shock front at a specific point in the driven section, the delta energy distribution was actively maintained in the present simulations.
Based on a readily available thermostat routine, a novel routine has been implemented in the LAMMPS software package which rescales the velocities of the driven gas particles to a target velocity, which was defined so that the driven gas kinetic energies of the Boltzmann and delta energy distributions are equal.
To allow unbiased interactions of the driven gas with the shock front, the routine was applied to a dynamically shrinking sub-section of the driven section.
This shrinking sub-section spanned the entire short axis of the driven section, had one side fixed to the endwall of the driven section, and had the other side move towards the endwall at a constant velocity.
The shrinking velocity was the shock velocity obtained through the simulations of the Boltzmann energy-distributed gas, which was larger than the shock velocity obtained through the simulations of the non-Boltzmann energy-distributed gas.
Through this, a buffer zone between the sub-section in which the delta energy distribution was actively maintained and the shock front was implemented.

For comparison, the ideal shock velocities for the above simulated conditions were calculated using the Shock- and Detonation Toolbox [48] in conjunction with the Cantera software package [49].
The effective pressures of the driver and driven sections used for shock velocity calculation were obtained through the Van-der-Waals equation of state with a molar co-volume of 0.03201 l/mol and a co-pressure of $135.5\;\text{Pa}\cdot\text{l}^2 / \text{mol}^2$ [50].
This leads to effective pressures of 821 and 48 bar for the driver and driven sections, respectively.
The pronounced real gas effects leading to the large driver pressure are acceptable in the present simulations, since the real gas interactions in the driver gas are not affecting the shock velocity in the driven section.
This shock velocity only depends on the ratio of effective pressures between driver and driven sections, the initial temperature, the molar mass, and the heat capacity ratio of the driven gas.
Since the driven gas is sufficiently well described by the ideal gas law (4 % deviation from ideal gas pressure), the presently calculated heat capacity ratios of the monoatomic ideal gases with a Boltzmann and a delta energy distribution can be used to determine the shock velocity in the driven section.
These shock velocities will be used for validating the present molecular dynamics simulations of Boltzmann energy-distributed gas and for testing if the same theories are applicable to non-Boltzmann energy-distributed gases.

§ RESULTS AND DISCUSSION

§.§ Non-Boltzmann Heat Capacity Ratio

The theoretical heat capacity ratio of non-Boltzmann energy-distributed monoatomic gases provided above is evaluated numerically for a Gaussian energy distribution with the limiting cases of the Boltzmann and delta energy distributions.
As mentioned above, the Gaussian energy distribution converges to the Boltzmann energy distribution if the Gaussian position $E_{\text{G}, 0}$ diverges to $-\infty$ and the average energies of the two distributions are set to be equal.
For the sake of compact representation, however, the ratio of the isochoric averaged squared energies (cf. equations <ref> and <ref>) of the Boltzmann and Gaussian distributions $\langle E_\text{B}^2 \rangle / \langle E_\text{G}^2 \rangle$ is used to characterize the deviation from the Boltzmann distribution.
This ratio yields unity in case the Gaussian distribution resembles the Boltzmann distribution and it yields $15/9 \approx 1.67$ in case the Gaussian distribution resembles the delta distribution.

Figure <ref> shows the heat capacity ratio of a Gaussian energy-distributed monoatomic gas $\gamma_\text{G}$ as function of the ratio of the averaged squared energies of the Boltzmann and Gaussian distributions.
The dashed lines represent the Boltzmann heat capacity ratio (cf. black dashed line, $\gamma_\text{B} = 5/3$), the Boltzmann limiting case of the Gaussian distribution (cf. blue dashed line, $\langle E_\text{G}^2 \rangle = \langle E_\text{B}^2 \rangle$), and the delta limiting case of the Gaussian distribution (cf. red dashed line, $\langle E_\text{G}^2 \rangle = \langle E_\delta^2 \rangle$).

Heat capacity ratio $\gamma_\text{G}$ of a Gaussian energy-distributed monoatomic gas as function of the ratio of the averaged squared energies $\langle E^2 \rangle$ of the Boltzmann and Gaussian energy distributions.
The horizontal dashed line represents the well-established heat capacity ratio of a Boltzmann energy-distributed monoatomic gas $\gamma_\text{B} = 5/3$.
The vertical dashed lines represent the limiting cases defined through the Boltzmann and delta energy distributions, $\langle E_\text{G}^2 \rangle = \langle E_\text{B}^2 \rangle$ and $\langle E_\text{G}^2 \rangle = \langle E_\delta^2 \rangle$, respectively.

As would be expected, the Gaussian heat capacity ratio resembles the Boltzmann and delta heat capacity ratios if the Gaussian distribution converges to either of these two limiting distributions.
However, there is an unexpected minimum for $\langle E_\text{B}^2 \rangle / \langle E_\text{G}^2 \rangle \approx 1.47$, for which the heat capacity ratio of the respective Gaussian distribution amounts to $\gamma_\text{G} \approx 0.92$.
This averaged squared energy ratio corresponds to a Gaussian position rather close to the upper limiting case of $E_0 = 3/2 \cdot k_\text{B} T$, with $E_{\text{G},0} = 0.9 \cdot E_0$.
It was tested if the observed minimum can be attributed to numerical uncertainties, yet changing the convergence threshold and the energy increment of the numerical integrations from $10^{-9}$ to $10^{-15}$ and from $1\;\text{cm}^{-1}$ to $0.1\;\text{cm}^{-1}$, respectively, had a negligible effect on $\gamma$ in the range of $10^{-5}$.
This means that the numerically calculated $\gamma_\text{G}$ curve is an accurate representation of the underlying theory and that the observed $\gamma_\text{G} < 1$ is theoretically valid.

From a classical thermodynamics perspective, the heat capacity ratio determines the volume response to adiabatic temperature changes as $\text{d}V/\text{d}T = (C_\text{p} - C_\text{V}) / p$.
In this classical context, a heated gas with $\gamma > 1$ would expand, with $\gamma = 1$ its volume would not change, and with $\gamma < 1$ it would contract.
In the present work, the effect of non-Boltzmann energy distributions on the adiabatic change of states was tested by randomly sampling velocities from Gaussian energy distributions and optimizing the Gaussian peak position $E_{\text{G}, 0}$ and width $b$ so that the pressure, the internal energy, and the entropy remained constant.
All presently tested Gaussian energy-distributed monoatomic gases expand with increasing temperature according to the well-established description of isentropic expansion within numerical uncertainties: $V_2 / V_1 = (T_1 / T_2)^{1 / (\gamma - 1)}$, with $\gamma = 5/3$.
This means that the non-Boltzmann heat capacity ratio must not be confused with the classical heat capacity ratio for thermal equilibrium.

The observed $\gamma < 1$ originates from the energy dependencies of the isochoric and isobaric densities of states.
While the isochoric density of states increases sub-proportional with energy, the isobaric density of states increases super-proportional with energy.
As a consequence, the relative population of low energy states is larger and the relative population of high energy states is smaller for isochoric conditions compared to isobaric conditions.
While this still leads to $C_\text{P} > C_\text{V}$ for probability distributions which monotonically decline with energy, such as the Boltzmann distribution or the Gaussian distribution with $\langle E_\text{B}^2 \rangle / \langle E_\text{G}^2 \rangle < 1.22$, it results in $C_\text{P} \leq C_\text{V}$ for a range of non-monotone Gaussian distributions with $\langle E_\text{B}^2 \rangle / \langle E_\text{G}^2 \rangle \gtrapprox 1.3$.
For these Gaussian energy distributions, the contributions of the densely populated low energy states to the isochoric heat capacity exceed the contributions of the sparsely populated high energy states to the isobaric heat capacity.

While the energy of the presently investigated probability distributions is set to be equal, the entropy is certainly not.
Here, Boltzmann's $H$-theorem [51] is utilized to obtain insights into the stability / entropy of the Boltzmann, Gaussian, and delta energy distributions.
Figure <ref> shows the $H$-function of the Boltzmann energy distribution (left-most limit), of the delta energy distribution (right-most limit), and of the Gaussian energy distributions interpolating between these two limiting cases.
The H-function is calculated as follows:

\begin{equation}
H = \sum\limits_k p_k \cdot \ln\left( p_k \right) ,
\end{equation}

with the probability $p_k$ of energy $E_k$, defined as $p_k = \rho_i(E_k) \cdot f_j(E_k) / \int\limits_0^\infty \rho_i(E) \cdot f_j(E) \cdot \text{d}E$.

$H$-function of the Boltzmann, Gaussian, and delta energy distributions as function of the ratio of the averaged squared energies $\langle E^2 \rangle$ of the Boltzmann and Gaussian energy distributions.
The vertical dashed lines represent the limiting cases defined through the Boltzmann and delta energy distributions, $\langle E_\text{G}^2 \rangle = \langle E_\text{B}^2 \rangle$ and $\langle E_\text{G}^2 \rangle = \langle E_\delta^2 \rangle$, respectively.

As would be expected, the Boltzmann energy distribution comes with the smallest value of the $H$-function.
The delta energy distribution comes with a $H$-function value of zero, which indicates that the entropy of the delta energy distribution is also zero.
This underlines that the Boltzmann and delta energy distributions are suitable limits for the investigation of non-Boltzmann energy distributions.
Starting from the delta energy distribution, the H-function of Gaussian energy distributions rapidly drops from zero to a monotonous convergence towards the $H$-function value of the Boltzmann energy distribution.
The rapid drop of the $H$-function close to the delta limiting case indicates that the delta energy distribution is entropically extremely unstable.

§.§ Atomistic Shock Wave Simulations

The present atomistic molecular dynamics simulations are used to test and validate the above theoretical findings on the non-Boltzmann heat capacity ratio.
The simulations resemble an ideal shock tube process, as the heat capacity ratio directly affects the velocity at which shock waves propagate through the driven gas.
First, potential non-idealities of the simulations will be discussed, followed by discussion of the simulation results and their comparison to theory.

Fixed particles have been used as diaphragm between the driver and driven sections of the simulation to allow independent thermalization of the two gases.
The abrupt removal of the diaphragm particles, however, leaves a 12.32 Å wide void between the two sections.
Based on the rather large effective pressure ratio between the driver and driven sections of about 17 it is assumed that the void between the driver and driven sections is solely penetrated by the driver gas.
This penetration through the driver gas into the void generates expansion and compression waves bouncing forth and back within this void during the first 12 ps of the simulation, which have been traced through the driver gas count.
Figure <ref> shows an exemplary driver gas count during the inception phase of the main simulation directly after abruptly removing the diaphragm particles.

Driver gas count during the inception phase directly after abrupt removal of the diaphragm particles which separated the driver and driven sections during thermalization.
The theoretical vacuum expansion velocity [52] for the simulated condition is provided as reference.

When loosely attributing kinks and jumps in the driver gas count to expansion and compression waves, one can find three characteristic velocities: 899-926 m/s, 616-680 m/s, and 513 m/s.
Each characteristic velocity comes in pairs, indicating that each wave which traverses the void towards the driven section is being reflected at the interface between the void and the driven section.
With each expansion and compression process, the void is filled with more driver gas particles, leading to a reduction in the characteristic velocities of the succeeding waves.
While the very first forth and back propagating wave has a characteristic velocity comparable to the theoretical vacuum expansion velocity $v_\text{vac} = 963\;\text{m}/\text{s}$ based on the formula provided by Greenspan and Butler [52], the characteristic velocities of the two succeeding forth and back propagating waves are already rather close to the shock velocity of the main shock wave propagating through the bulk of the driven gas after the inception phase ($v_\text{S} = 545.56\;\text{m}/\text{s}$ for the particular example in Figure <ref>, cf. Table <ref> for all resutls).

After the third forth and back propagating wave has traversed the void between the driver and driven sections, the driver gas flow is established and the driver gas count remains constant at the time scale relevant for the inception phase.
The large  volume ratio of the driven section to the driver section of 9, however, causes a significant drop in driver gas concentration at later times, as visible in the sequence of frames provided in Figure <ref>.
Neither the dropping pressure in the driver section nor the observed bouncing expansion and compression waves during the inception phase are expected to affect the shock velocity of the main shock wave freely propagating through the bulk of the driven gas.

Each simulation is post-processed based on the driver and driven gas particle counts along the long axis of the simulation, which ranges from -1000 Å to 9000 Å.
In order to improve statistics, ten subsequent simulation steps are always merged into a single frame, yielding a total of 119410 driver gas particles and 108520 driven gas particles for evaluating the particle distribution along the x-axis.
Figure <ref> shows two sequences of such particle distributions together with snapshots of the simulation boxes for a Boltzmann energy-distributed driven gas (<ref>) and for a non-Boltzmann energy-distributed driven gas (<ref>).
Note that both simulations have an identical random number generator seed, thus are identical except for the energy distribution of the driven gas preceding the shock wave.
The particle counts are discretized along the long axis in 100 Å increments with the driver gas particles shown in red and the driven gas particles shown in blue for both the simulation snapshots and the particle count plots.
The red left pointing arrow and the blue right pointing arrow in the particle count plots assign the left and right y-axes of the plots to the driver and driven gas particle counts, respectively.
The upward pointing red and blue arrows directly below the simulation snapshots in turn indicate the contact and shock wave positions, respectively.
These positions were determined through the steps of the double sigmoid fits to the driven gas particle counts (black lines).
The full video of the presently provided simulation frames can be found in the supplementary information.

Sequences of frames taken from the post-processing of the simulations of Boltzmann and non-Boltzmann energy-distributed driven gases with the same random number generator seed.

The first observation from Figure <ref> is that the shock positions differ between the Boltzmann and non-Boltzmann energy-distributed driven gases.
In the Boltzmann case, the shock velocity extracted from the provided example is 545.56 m/s, while the shock velocity in the corresponding non-Boltzmann case is only 477.18 m/s.
Interestingly, the shock front appears to be thicker in the Boltzmann case compared to the non-Boltzmann case, with a sigmoid width of 70.12 Å for the Boltzmann case and 28.64 Å for the non-Boltzmann case, averaged over the entire simulations.
The differences in shock velocities and in shock front widths can be qualitatively explained by the differences in the energy distributions.
The presently applied non-Boltzmann energy distribution closely resembles a delta energy distribution, thus the driven gas particles have almost exclusively one randomly oriented velocity.
This means that the driven gas particles do not have high or low velocity particles, as it is the case for a Boltzmann energy-distributed driven gas.
The high velocity gas particles facilitate a faster shock propagation, leading to the observed larger shock velocity in the Boltzmann case.
At the same time, the low velocity particles are trailing the shock front, leading to a more stretched shock front in the Boltzmann case.

For each replica simulation of the Boltzmann and the non-Boltzmann cases, the contact and shock wave positions are obtained as described above and traced over time.
Figure <ref> shows the wave position profiles for the two exemplary simulations provided in Figure <ref>.
The reported wave velocities and uncertainties are obtained via linear fitting, with the uncertainties being explicitly for the slope of the linear fit.

Contact and shock wave positions over simulation time for two exemplary simulations of a comparable Boltzmann and non-Boltzmann case.

This linear fitting procedure is carried out for all replica simulations and the resulting contact and shock wave velocities with their respective fitting uncertainties are tabulated in Table <ref>.
The wave velocities are sorted by the energy distribution case of the driven gas, either Boltzmann or non-Boltzmann.

Shock wave velocities $v_\text{S}$ and contact wave velocities $v_\text{C}$ of the Boltzmann and non-Boltzmann driven gases.

Replica     $v_{\text{S}, \text{sim}}$     $v_{\text{S}, \text{sim}}$ err.     $v_{\text{C}, \text{sim}}$     $v_{\text{C}, \text{sim}}$ err.

m/s     m/s     m/s     m/s

1     545.56     1.99     185.77     1.25

2     555.43     2.58     182.73     1.10

3     558.66     1.63     180.83     1.16

4     549.24     1.75     181.44     1.06

5     536.65     1.52     183.45     1.33

1     477.18     2.02     197.87     0.80

2     502.92     2.33     188.64     1.05

3     516.63     3.62     191.27     1.08

4     515.82     3.13     188.12     1.05

5     493.42     2.17     182.76     1.53

The differences observed in Figure <ref> for a single example of a Boltzmann and a non-Boltzmann case prevail in the other replica simulations, as well as in the averaged velocities.
The five simulations with the Boltzmann energy-distributed driven gas give an average shock velocity of $v_{\text{S}, \text{B}, \text{sim}} = 549.11\pm6.40\;\text{m}/\text{s}$ and the five simulations with the non-Boltzmann energy-distributed driven gas give an average shock velocity of $v_{\text{S}, \delta, \text{sim}} = 501.19\pm12.72\;\text{m}/\text{s}$.
This means that the presently employed non-Boltzmann energy distribution, which closely resembles the delta distribution, hinders the free propagation of shock waves through the monoatomic gas by about 9 % with respect to the shock velocity.
Note that convergence of the presented results has been tested by simulating with a doubled temporal resolution and post-processing with doubled spatial resolution, resulting in -0.08 % and -0.1 % deviations in the shock velocity, respectively.
In addition, the total energy of the present NVE simulations should be constant and was found to deviate by 0.13 % over the entire simulation time at most.

Therefore, the presented results are assumed to be sufficiently converged and the observed scatter is solely attributed to the statistical uncertainties inherent to molecular dynamics simulations [53].
For the non-Boltzmann molecular dynamics simulation, the speed of the dynamically shrinking sub-section which is used to maintain the delta energy distribution has been modified as well, resulting in a shock velocity within the aforementioned methodological uncertainty.

Interestingly, the contact wave appears to accelerate from an average of $v_{\text{C}, \text{B}, \text{sim}} = 183.04\pm1.25\;\text{m}/\text{s}$ to $v_{\text{C}, \delta, \text{sim}} = 189.73\pm3.87\;\text{m}/\text{s}$ when applying a delta energy distribution to the driven gas.
This effect amounts to about 4 % and is larger than the observed statistical uncertainties of the simulations.
Noteworthy, this effect is trailing the incident shock wave, yet the gas behind the shock wave would typically be expected to be unaffected by the energy distribution in front of the shock wave.

The averaged shock velocities of the Boltzmann and non-Boltzmann cases agree well with the ideal shock calculations using the incident shock prediction of the Shock- and Detonation Toolbox [48], which yields $v_{\text{S}, \text{B}, \text{calc}} = 548\;\text{m}/\text{s}$ for $\gamma = 5/3$ and $v_{\text{S}, \gamma \to 1, \text{calc}} = 494\;\text{m}/\text{s}$ when slowly converging $\gamma$ to unity.
While the shock velocity is ill-defined for $\gamma = 1$, it converges with $\gamma \to 1$.
Although the value for $\gamma \to 1$ is slightly below the average simulation result, it is clearly within the statistical uncertainty of the simulations.
The comparison of the ideal shock calculation with the present simulations show that a non-Boltzmann heat capacity ratio can be combined with the standard description of ideal shock waves, which has originally been formulated for Boltzmann energy-distributed gases.
As mentioned above, however, the non-Boltzmann heat capacity ratio should not be confused with the classical heat capacity ratio for thermal equilibrium.

§ CONCLUSIONS

In the present study, a formulation of the heat capacity ratio of non-Boltzmann energy-distributed monoatomic gas has been derived from first principles.
This theory is applied to non-Boltzmann energy distributions which all exhibit the same total energy and are limited by the Boltzmann energy distribution as the equilibrium case and the delta energy distribution as the most extreme non-equilibrium case.
A continuous description of distributions between the two limiting cases is facilitated using a Gaussian distribution with variable peak position.
Applying the presented theory in the Boltzmann case yields the well-established heat capacity ratio of $\gamma_\text{B} = 5/3$, while in the delta case it yields $\gamma_\delta \to 1$.
Interestingly, applying the presented theory to the Gaussian case with a peak position close to the delta peak position, the heat capacity ratio appears to drop below unity.
Investigation of the isentropic thermal expansion behavior of non-Boltzmann energy-distributed monoatomic gas, however, revealed that the energy distribution does not affect thermal expansion and that a monoatomic gas always expands with a $\gamma = 5/3$ in the isentropic case.
To stress it once more, the non-Boltzmann heat capacity ratio should not be confused with the classical heat capacity ratio for thermal equilibrium.

The presented atomistic molecular dynamics simulations resemble a shock tube process and have been carried out for Boltzmann and non-Boltzmann energy-distributed driven gases.
The contact wave and shock wave velocities extracted from these simulations reveal a consistent difference between the Boltzmann and non-Boltzmann cases.
While the shock wave propagates about 9 % slower through the non-Boltzmann driven gas, the contact wave appears to be faster by 4 % when trailing a shock wave through a non-Boltzmann driven gas.
The difference between the Boltzmann and non-Boltzmann shock wave velocities are reproduced through the ideal shock equations when utilizing the heat capacity ratios of $\gamma = 5/3$ and $\gamma \to 1$ for the Boltzmann and non-Boltzmann cases, respectively, as provided by the present theory.
Firstly, this validates the present theory for non-Boltzmann heat capacity ratios in the context of shock propagation and secondly, it appears that non-Boltzmann heat capacity ratios can be combined with the classical shock wave formulations for describing non-Boltzmann gas dynamics.

Building on the presented theory, the non-Boltzmann heat capacity ratio should be formulated for polyatomic compounds by replacing the monoatomic density of states with a polyatomic formulation in future studies.
Further simulations of gas dynamical processes, both atomistic and continuous, would be necessary to explore the limits of the present theory.
The non-Boltzmann energy distributions in between the Boltzmann and the delta energy distributions are of particular interest, but molecular dynamics simulations will require rather sophisticated thermostating routines for maintaining such non-Boltzmann energy distributions.
Ultimately, experiments need to be designed which would allow to test the presented theoretical and computational findings.

The present work provides novel insights into gas dynamics of non-Boltzmann gases and proposes interesting and partly counter-intuitive aspects of the heat capacity ratio.
Since the heat capacity ratio is a key property for describing gas dynamical processes, the present findings might help to improve the understanding of gas dynamical phenomena and could potentially allow discovering novel phenomena.

Supplementary information
Animated sequences of exemplary Boltzmann and non-Boltzmann molecular dynamics simulations, corresponding to the frames presented in Figure <ref>.

§ DECLARATIONS

Conflict of interest The author has no financial or proprietary interests in any material discussed in this article.

[1]
R. R. Boyce, D. R. N. Pulford, A. F. P. Houwing, and Ch. Mundt.
Rotational and vibrational temperature measurements using CARS in a
hypervelocity shock layer flow and comparison with CFD calculations.
Shock Waves, 6:41–51, 1996.

[2]
A. Kosareva, O. Kunova, E. Kustova, and E. Nagnibeda.
Four-temperature kinetic model for CO2 vibrational relaxation.
Phys. Fluids, 33:016103, 2021.

[3]
V. M. Doroshenko, N. N. Koudriavtsev, and V. V. Smetanin.
Nonequilibrium excitation of internal molecular degrees of freedom in
the shock layer during hypersonic flight.
Shock Waves, 2:139–146, 1992.

[4]
M. Louviot, N. Suas-David, V. Boudon, R. Georges, M. Rey, and S. Kassi.
Strong thermal nonequilibrium in hypersonic CO and CH4
probed by CRDS.
J. Chem. Phys., 142:214305, 2015.

[5]
E. Dudás, N. Suas-David, S. Brahmachary, Y. Kukarni, A. Benidar, S. Kassi,
C. Charles, and R. Georges.
High-temperature hypersonic laval nozzle for non-LTE cavity
ringdown spectroscopy.
J. Chem. Phys., 152:134201, 2020.

[6]
C. Cercignani.
Boltzmann Equation and Rarefied Gas Dynamics, pages 75–89.
Springer, Boston, United States, 1985.

[7]
K. Huang.
Statistical Mechanics, page 53.
Wiley, New York, United States, 1987.

[8]
G. A. Bird.
Molecular Gas Dynamics and the Direct Simulation of Gas Flows.
Oxford University Press, Oxford, 1994.

[9]
M. A. Gallies.
Graeme a. bird.
<https://www.osti.gov/servlets/purl/1595027>, 2019.

[10]
M. Torrilhon.
Modeling nonequilibrium gas flow based on moment equations.
Annu. Rev. Fluid Mech., 48:429–458, 2016.

[11]
J. D. Au, M. Torrilhon, and W. Weiss.
The shock tube study in extended thermodynamics.
Phys. Fluids, 13:2423, 2001.

[12]
P. Ván.
Nonequilibrium thermodynamics: emergent and fundamental.
Phil. Trans. R. Soc., 378:20200066, 2020.

[13]
H. C. Öttinger, H. Struchtrup, and M. Torrilhon.
Formulation of moment equations for rarefied gases within two
frameworks of non-equilibrium thermodynamics: Ret and generic.
Phil. Trans. R. Soc., 378:20190174, 2020.

[14]
I. Müller and T. Ruggeri.
Rotational extended thermodynamics.
Springer, New York, United States, 1998.

[15]
M. Grmela and H. C. Öttinger.
Dynamics and thermodynamics of complex fluids. I. development of a
general formalism.
Phys. Rev. E, 56:6620, 1997.

[16]
G. Lebon, D. Jou, and J. Casas-Vázquez.
Understanding Non-equilibrium Thermodynamics.
Springer, Heidelberg, Germany, 2008.

[17]
C. Beck and E. G. D. Cohen.
Physica A Stat. Mech. Appl., 322:267–275, 2003.

[18]
U. Seifert.
Stochastic thermodynamics: principles and perspectives.
Eur. Phys. J. B, 64:423–431, 2008.

[19]
B. Hafskjold and S. K. Ratkje.
Criteria for local equilibrium in a system with transport of heat and
J. Stat. Phys., 78:463–494, 1995.

[20]
R. N. Schwartz, Z. I. Slawsky, and K. F. Herzfeld.
Calculation of vibrational relaxation times in gases.
J. Chem. Phys., 20:1591–1599, 1952.

[21]
M. P. Burke, C. F. Goldsmith, Y. Georgievskii, and S. J. Klippenstein.
Towards a quantitative understanding of the role of non-boltzmann
reactant distributions in low-temperature oxidation.
Proc. Combust. Inst., 35:205–213, 2015.

[22]
C. F. Goldsmith, M. P. Burke, Y. Georgievskii, and S. J. Klippenstein.
Effect of non-thermal product energy distribution on
ketohydroperoxide decomposition kinetics.
Proc. Combust. Inst., 35:283–290, 2015.

[23]
F. L. Dryer, D. Naegeli, and I. Glassman.
Temperature dependence of the reaction CO + OH = CO2 + H.
Combust. Flame, 17:270–272, 1971.

[24]
A. M. Dean.
Predictions of pressure and temperature effects upon radical addition
and recombination reactions.
J. Phys. Chem., 89:4600–4608, 1985.

[25]
P. R. Westmoreland, J. B. Howard, J. P. Longwell, and A. M. Dean.
Prediction of rate constants for combustion and pyrolysis reactions
by bimolecular QRRK.
AIChE J., 32:1971–1979, 1986.

[26]
N. J. Labbe, R. Sivaramakrishnan, C. F. Goldsmith, Y. Georgievskii, J. A.
Miller, and S. J. Klippenstein.
Weakly bound free radicals in combustion: "prompt" dissociation of
formyl radicals and its effect on laminar flame speeds.
J. Phys. Chem. Lett., 7:85–89, 2016.

[27]
N. J. Labbe, R. Sivaramakrishnan, C. F. Goldsmith, Y. Georgievskii, J. A.
Miller, and S. J. Klippenstein.
Ramifications of including non-equilibrium effects for HCO in
flame chemistry.
Proc. Combust. Inst., 36:525–532, 2017.

[28]
M. Döntgen, L. C. Kröger, and K. Leonhard.
Hot $\beta$-scission of radicals formed via hydrogen abstraction.
Proc. Combust. Inst., 36:135–142, 2017.

[29]
M. Döntgen and K. Leonhard.
Discussion of the separation of chemical and relaxational kinetics of
chemically activated intermediates in master equation simulations.
J. Phys. Chem. A, 121:1593–1570, 2017.

[30]
A. Wildenberg, M. Döntgen, I. S. Roy, C. Huang, B. Lefort, L. Le Moyne,
A. Kéromnès, K. Leonhard, and K. A. Heufer.
Solving the riddle of the high-temperature chemistry of
Proc. Combust. Inst., 39, 2023.

[31]
O. Kunova, E. Kustova, M. Mekkonoshina, and E. Nagnibeda.
Non-equilibrium kinetics, diffusion and heat transfer in shock heated
flows of N2/N and O2/O mixtures.
Chem. Phys., 463:70–81, 2015.

[32]
D. S. Kravchenko, E. V. Kustova, and M. Y. Melnik.
Modeling of state-to-state oxygen kinetics behind reflected shock
Vestn. St. Petersbg., 55:281–289, 2022.

[33]
M. Y. Melnik and E. V. Kustova.
Impact of electronic excitation on the state-to-state
vibrational-chemical CO kinetics.
J. Phys. Conf. Ser., 2308:012014, 2022.

[34]
Y. Sakamura.
A master equation study of vibration-dissociation coupling in
shock-heated oxygen molecules.
Shock Waves, 12:361–373, 2003.

[35]
I. N. Kadochnikov and I. V. Arsentiev.
Modelling of vibrational nonequilibrium effects on the H2–air
mixture ignition under shock wave conditions in the state-to-state and mode
Shock Waves, 30:491–504, 2020.

[36]
M. F. Campbell, K. G. Owen, D. F. Davidson, and R. K. Hanson.
Dependence of calculated postshock thermodynamic variables on
vibrational equilibrium and input uncertainty.
J. Thermophys. Heat Transf., 31:586–608, 2017.

[37]
A. P. Bedin.
Gasdynamic phenomena accompanying the motion of shock waves and
objects in a low-temperature nonequilibrium plasma.
Tech. Phys. Lett., 24:725–727, 1998.

[38]
B. N. Ganguly, P. Bletzinger, and A. Garscadden.
Shock wave damping an dispersion in nonequilibrium low pressure argon
Phys. Lett. A, 230:218–222, 1997.

[39]
V. Bityurin, A. Klimov, S. Leonov, Y. Brovkin, Y. Kolesnichenko, N. Popov, and
D. M. Van Wie.
Shock wave structure and velocity at propagation through
non-homogeneous plasma.
In 31${}^\text{st}$ AIAA Plasmadynamics and Lasers Conference,

[40]
J. Poggie.
Modeling the effects of weak ionization on supersonic flow and shock
In 37${}^\text{th}$ AIAA Aerospace Sciences Meeting and
Exhibit, 1999.

[41]
N. E. Molevich.
Acoustical properties of nonequilibrium media.
In 42${}^\text{nd}$ AIAA Aerospace Sciences Meeting and
Exhibit, 2004.

[42]
Q. Zhou, Z. Dong, and W. Yang.
Revisiting the thermal effect on shock wave propagation in weakly
ionized plasmas.
Phys. Plasmas, 23:073508, 2016.

[43]
P. Atkins and J. De Paula.
Physical Chemistry, page 64.
Oxford University Press, Oxford, 2010.

[44]
M. Döntgen.
Analytic energy-level densities of separable harmonic oscillators
including approximate hindered rotor corrections.
AIP Adv., 6:095318, 2016.

[45]
S. J. Plimpton.
Fast parallel algorithms for short-range molecular dynamics.
J. Comput. Phys., 117:1–19, 1995.

[46]
J. E. Lennard-Jones.
Proc. Phys. Soc., 43:461, 1931.

[47]
H.-H. Carstensen and A. M. Dean.
Chapter 4 the kinetics of pressure-dependent reactions.
In Robert W. Carr, editor, Comprehensive Chemical Kinetics,
volume 42, pages 101–187. Elsevier, 2007.

[48]
J. E. Shepherd.
Shock and detonation toolbox.
https://shepherd.caltech.edu/EDL/publicresources.html, 2014.
accessed 2021/02/11.

[49]
D. G. Goodwin, H. K. Moffat, and R. L. Speth.
Cantera: An object-oriented software toolkit for chemical kinetics,
thermodynamics, and transport processes.
<http://www.cantera.org>, 2016.
Version 2.2.1.

[50]
Engineering ToolBox.
Non-ideal gas – van der waal's equation and constants.
accessed 2022/12/14.

[51]
L. Boltzmann.
Weitere Studien über das Wärmegleichgewicht unter
Sitzungsberichte Akad. Wiss., 66:275–370, 1872.

[52]
H. P. Greenspan and D. S. Butler.
On the expansion of a gas into vacuum.
J. Fluid Mech., 13:101–119, 1962.

[53]
L. C. Kröger, W. A. Kopp, and M. Döntgen.
Assessing statistical uncertainties of rare events in reactive
molecular dynamics simulations.
J. Chem. Theory Comput., 13:3955–3960, 2017.