# ChemPlasKin: a general-purpose program for unified gas and plasma kinetics
simulations

Xiao Shaoa,∗, Deanna A. Lacostea, Hong G. Ima
aCCRC, King Abdullah University of Science and Technolody (KAUST), Thuwal
23955-6900, Saudi Arabia
∗E-mail<EMAIL_ADDRESS>
###### Abstract

This work introduces ChemPlasKin, a freely accessible solver optimized for
zero-dimensional (0D) simulations of chemical kinetics of neutral gas in non-
equilibrium plasma environments. By integrating the electron Boltzmann
equation solver, CppBOLOS, with the open-source combustion library, Cantera,
at the source code level, ChemPlasKin computes time-resolved evolution of
species concentration and gas temperature in a unified gas-plasma kinetics
framework. The model allows high fidelity predictions of both chemical thermal
effects and plasma-induced heating, including fast gas heating and slower
vibrational-translational relaxation processes. Additionally, a new heat loss
model is developed for nanosecond pulsed discharges, specifically within pin-
pin electrode configurations. With its versatility, ChemPlasKin is well-suited
for a wide range of applications, from plasma-assisted combustion (PAC) to
fuel reforming. In this paper, the reliability, accuracy and efficiency of
ChemPlasKin are validated through a number of test problems, demonstrating its
utility in advancing gas-plasma kinetic studies.

Keywords: plasma-assisted combustion (PAC); ion chemistry; fuel reforming;
electron-impact reactions; reaction kinetics

Source code: https://github.com/ShaoX96/ChemPlasKin

##  1 Introduction

Non-equilibrium plasma has gained an increasing interest within both the
combustion and plasma research communities, owing to its potential to enhance
combustion characteristics and fuel reforming [1, 2, 3]. The integration of
plasma actuation into reacting flows, via kinetic, thermal, and hydrodynamic
effects, presents a complex multi-physics challenge that requires a
synergistic approach combining both experimental and computational
methodologies to gain a comprehensive understanding. Due to the formidable
computational complexity and expense associated with higher-dimensional
models, detailed chemical kinetic analyses of systems that couple neutral gas
with plasma have predominantly been confined to zero or one-dimensional
simulations. Such 0D simulations, enriched with detailed gas-plasma kinetics,
are not only foundational for the development of kinetic mechanisms but are
also pivotal in identifying important pathways of plasma energy transfer. This
insight is crucial for creating reduced-order plasma models that can be
integrated efficiently into computational fluid dynamics (CFD) solvers [4, 5].

A combined kinetic solver requires the determination of species reaction rates
for the gas-phase and plasma. For the latter, the electron energy distribution
function (EEDF) is crucial for determining reaction rate coefficients of
inelastic electron-impact collisions and calculating electron temperature
($T_{e}$). In non-equilibrium plasma, EEDF deviates from Maxwellian, and needs
to be determined by solving the electron Boltzmann equation (EBE). Table 1
summarizes various modeling approaches for computing chemical kinetics in gas-
phase coupled with non-equilibrium plasma available in the literature.
Broadly, solver development methodologies fall into three categories:
tabulation, operator-splitting, and full coupling. Solvers 1-4 employ either
the Bolsig+ EBE solver [6] or the ZDPlasKin plasma kinetics solver [7] to pre-
calculate electron-impact reaction rates to build look-up tables or polynomial
fits as functions of reduced electric field ($E/N$) or $T_{e}$, for use with
CHEMKIN [8], the standard in gas-phase chemical kinetics software. While
straightforward, tabulation may be subjected to large errors with significant
gas composition changes. Solvers 5-13, following Lefkowitz et al. [9],
integrate ZDPlasKin with chemical kinetics solvers like CHEMKIN or Cantera
[10], alternating the integration of plasma and neutral gas kinetics with
stepwise exchange of species concentration and temperature data. Solver 14 by
Cheng et al. [11] is the only approach known to date that directly couples the
BOLOS EBE solver [12] with Cantera, enabling unified gas-plasma kinetics
simulations, although the details of its code implementation were not
elaborated.

While CHEMKIN and Cantera are well-established open-source modules for
combustion community, and ZDPlasKin is for the plasma community, the
integration of the two schools for a unified version of gas-plasma kinetics
solver is not straightforward, as evidenced by the fact that most of the tools
shown in Table 1 are in-house codes and lack open availability. To address
this deficiency, we introduce ChemPlasKin, an open-source code optimized for
simulating neutral gas-phase chemical reactions in conjunction with non-
thermal plasma chemistry within a unified code module framework. Similar to
the CERFACS code [11], ChemPlasKin is developed with a particular interest in
for nanosecond repetitively pulsed (NRP) plasma applications, which exhibits a
large spectrum of time scales, offering new capabilities for plasma-assisted
ignition (PAI) and fuel reforming applications.

The main content of this paper is divided into two sections. Section 2
outlines the methodologies, including the governing equations, the heat loss
model, and code development. Section 3 presents extensive code validations
that demonstrate the consistency and efficiency of ChemPlasKin in comparison
with results from the literature. It also highlights the solver’s capability
in predicting both fast and slow gas heating processes, and assesses the
impact of the newly introduced heat loss model.

Table 1:  Literature review of publications involving in-house 0D gas-plasma kinetics solvers $\\#$ | Origin |  Method |  Topic | Discharge | Year
---|---|---|---|---|---
1 | Ohio State |  Bolsig+ precalculation + ChemKin-Pro |  Air/$\mathrm{H_{2}}$ and air/hydrocarbon kinetics [13] | NRP | 2015
2 | UC Berkeley |  ZDPlasKin precalculation + CHEMKIN |  $\mathrm{CH_{4}}$/air ignition [14] | NRP | 2016
3 | UT Austin |  Bolsig+ precalculation + CHEMKIN |  $\mathrm{CH_{4}}$/air and $\mathrm{C_{2}H_{4}}$/air ignition [15] | NRP | 2021
4 | TU/e |  Bolsig+ precalculation + Sundials IDA |  $\mathrm{CH_{4}}$ and $\mathrm{H_{2}}$ oxidation in Ar [16] | DBD | 2023
5 | Princeton |  ZDPlasKin + CHEMKIN |  $\mathrm{CH_{4}/O_{2}/He}$ oxidization [9] | NRP | 2015
|  |  |  $\mathrm{CH_{4}}$/air ignition [17] | NRP | 2018
|  |  |  $\mathrm{H_{2}/O_{2}/He}$ ignition [18] | NRP/DC | 2019
|  |  |  $\mathrm{CH_{4}/O_{2}/He}$ ignition [19] | NRP/DC | 2019
|  |  |  $\mathrm{H_{2}}$/air ignition [20] | NRP | 2020
|  |  |  N-dodecane/$\mathrm{O_{2}/N_{2}}$ kinetics [21] | NRP | 2021
|  |  |  $\mathrm{NH_{3}/O_{2}/He}$ ignition [22] | NRP | 2021
|  |  |  $\mathrm{NH_{3}}$/air ignition [23] | NRP | 2022
|  |  |  $\mathrm{N_{2}O/NO_{x}}$ mechanism in $\mathrm{NH_{3}}$ oxidation [24] | DBD | 2023
|  |  |  $\mathrm{NH_{3}}$/air ignition [25] | NRP | 2024
|  |  |  n-pentane/air oxidation [26] | NRP | 2024
6 | USC |  ZDPlasKin + Cantera |  DME/$\mathrm{O_{2}/Ar}$ and $\mathrm{C_{3}H_{8}/O_{2}/Ar}$ ignition [27] | nanosecond | 2017
7 | Tsinghua |  ZDPlasKin + CHEMKIN |  DME oxidation kinetics [28] | DBD | 2021
8 | UMN |  ZDPlasKin + CHEMKIN |  $\mathrm{NH_{3}}$ pyrolysis and combustion [29] | NRP | 2021
|  |  |  $\mathrm{NH_{3}}$/air ignition and $\mathrm{NO_{x}}$ emission [30] | NRP | 2022
|  |  |  Plasma-based global pathway analysis [31] | NRP | 2023
9 | KAUST |  ZDPlasKin + CHEMKIN |  Lean $\mathrm{H_{2}/O_{2}}$ kinetics [32] | DBD | 2022
|  |  |  Rich $\mathrm{H_{2}/O_{2}}$ kinetics [33] | DBD | 2022
|  |  |  $\mathrm{NH_{3}}$ cracking [34] | DBD | 2023
|  |  |  $\mathrm{O_{3}}$ kinetics [35] | DBD | 2023
10 | Birmingham |  ZDPlasKin + Cantera |  $\mathrm{NH_{3}/N_{2}/O_{2}/He}$ combustion [36] | NRP | 2022
11 | MIT |  ZDPlasKin + Cantera |  1D $\mathrm{CH_{4}/air}$ flame [37] | NRP-DBD | 2023
12 | XJTU |  ZDPlasKin + in-house chemistry solver |  $\mathrm{NH_{3}/N_{2}/O_{2}}$ ignition [38] | NRP/DC | 2023
13 | SDU |  ZDPlasKin + CHEMKIN |  $\mathrm{CH_{4}/O_{2}/He}$ ignition [39] | NRP-SDBD | 2024
14 | CERFACS |  BOLOS + Cantera |  $\mathrm{CH_{4}}$/air mechanism [11] | NRP | 2022
|  |  |  $\mathrm{CH_{4}}$/air phenomenological model [5] | NRP | 2023

##  2 Methodology

###  2.1 Governing equations

The governing equations build upon the previous work by Cheng et al. [11] and
are presented in a more generalized format. ChemPlasKin integrates the mass
fraction $Y_{k}$ and gas temperature $T_{\text{gas}}$ over time $t$ in a 0D
neutral gas-plasma reactor for a total of $N$ species:

$\frac{dY_{k}}{dt}=\frac{W_{k}}{\rho}\dot{\omega}_{k},$ (1) $\displaystyle\rho
c_{v}\frac{dT_{\text{gas}}}{dt}$
$\displaystyle=-\sum_{k=1}^{N}\dot{\omega}_{k}u_{k}W_{k}+\dot{E}^{p}+\sum_{k}{(-\dot{E}_{\text{vib}}^{k}+\dot{R}_{\text{VT}}^{k})},\quad\text{(constant
volume)}$ (2a) $\displaystyle\rho c_{p}\frac{dT_{\text{gas}}}{dt}$
$\displaystyle=-\sum_{k=1}^{N}\dot{\omega}_{k}h_{k}W_{k}+\dot{E}^{p}+\sum_{k}{(-\dot{E}_{\text{vib}}^{k}+\dot{R}_{\text{VT}}^{k})}.\quad\text{(constant
pressure)}$ (2b)

Here, $W_{k}$ and $\dot{\omega}_{k}$ denote the molecular weight and molar
production rate of species $k$, respectively. $\rho$ represents the density.
The terms $c_{v}$ and $c_{p}$ signify the mass heat capacities at constant
volume and constant pressure, respectively. The term $u_{k}$ refers to the
internal energy of species $k$, and $h_{k}$ represents the enthalpy of species
$k$. The total plasma energy $\dot{E}^{p}$ is derived from non-elastic
electron collision processes $\mathcal{P}$:

$\dot{E}^{p}=\sum_{j\in\mathcal{P}}\varepsilon_{\text{th}}^{j}q_{j},$ (3)

where $\varepsilon_{\text{th}}^{j}$ and $q_{j}$ are the threshold energy and
net molar production rate for each process $j$, respectively.

A significant portion of plasma energy, not directly contributing to fast gas
heating, is stored in the vibrational states of species $k$, denoted by
$\dot{E}_{\text{vib}}^{k}$. This energy is released to facilitate slow gas
heating, primarily through vibrational-translational (V-T) relaxation, at a
rate $\dot{R}_{\text{VT}}^{k}$. Given the system’s non-equilibrium nature, a
dedicated equation for the vibrational energy $e_{\text{vib}}^{k}$ is
necessary:

$\frac{de_{\text{vib}}^{k}}{dt}=\dot{E}_{\text{vib}}^{k}-\dot{R}_{\text{VT}}^{k}.$
(4)

$\dot{R}_{\text{VT}}^{k}$ is modeled using the V-T relaxation timescale
$\tau_{\text{VT}}^{k}$:

$\dot{R}_{\text{VT}}^{k}=\frac{e_{\text{vib}^{k}}}{\tau_{\text{VT}}^{k}}.$ (5)

Extending the two-component mixture relaxation formula in [40] to multiple
species, we obtain:

$\left(\tau_{\text{VT}}^{k}\right)^{-1}=\sum_{m}\left(X_{m}/\tau_{\text{VT}}^{k,m}\right),$
(6)

where $X_{m}$ represents the mole fraction of collider $m$ and
$\tau_{\text{VT}}^{k,m}$ the relaxation time for oscillator $k$ in a high
dilution of $m$. The condition $m=k$ indicates relaxation in pure gas $k$.
Experimental data fitting lines in [40] are expressed as:

$p_{0}\tau_{\text{VT}}^{k,m}=\text{exp}\left(a(T^{-1/3}-b)-18.42\right),\
\mathrm{[atm\cdot s]},$ (7)

with $p_{0}$ denoting total pressure in atm, and $(a,b)$ representing the
fitting parameter pair for a $k\text{-}m$ mixture. Specifically, for air-
dominated systems where vibrational energy is primarily stored in
$\mathrm{N_{2}}(v)$ for $E/N>50$ Td, it is pertinent to consider
$\mathrm{N_{2},O_{2},}$ and $\mathrm{O}$ as major collision partners. For
$m=\mathrm{N_{2}}$, $(a,b)=(220,0.03)$ are directly obtained from [40]. For
$m=\mathrm{O_{2}}$, $(a,b)=(162,0.03)$ are inferred from experimental results
at 2500 K [41]. The most efficient V-T relaxation pathway involves reactions
with $\mathrm{O(^{3}P)}$ atoms, characterized by a rate constant of $4.5\times
10^{-15}(T/300)^{2.1}\ [\mathrm{cm^{3}/s]}$ [42], setting
$\tau_{\text{VT}}^{\mathrm{{N_{2},O}}}$ to
$488.5/(p_{0}T^{1.1}X_{\mathrm{O}})$. Note that these $\tau_{\text{VT}}^{k,m}$
parameters, are different from those in the phenomenological model by Castela
et al.[4], underscoring the need for validation and refinement for more
precise V-T relaxation predictions.

The model’s approach to slow heating via V-T relaxation, encapsulated in
equations (2a), (2b), (4) and (5), consolidates all vibrational states of
species $k$ into a singular variable, $e_{\text{vib}}^{k}$. This aggregation
significantly simplifies the species count, facilitating practicality in
multi-dimensional simulations. However, users have the option to model each
vibrational state explicitly along with their specific V-T relaxation
reactions. In such scenarios, the slow gas heating term
$\sum_{k}{(-\dot{E}_{\text{vib}}^{k}+\dot{R}_{\text{VT}}^{k})}$ in equation
(2a) and (2b) is unnecessary. Note that the operator-splitting solvers (5-13)
listed in Table 1 utilize a model proposed by Flitti & Pancheshnyi [43] to
describe gas heating. In this model, external power deposited into the system
is distributed among gas heating, electron energy, and chemical energy.
Similarly, ChemPlasKin incorporates this model option for gas temperature
calculations, analogous to the formulas used in Lefkowitz et al. [9]. For
constant volume configurations, the governing equation reads:

$\rho
c_{v}\frac{dT_{\text{gas}}}{dt}=-\sum_{k=1}^{N}\dot{\omega}_{k}u_{k}W_{k}+eN_{e}v_{e}E-\frac{3}{2}RT_{e}\dot{\omega}_{e},$
(8)

where $e$ represents the elementary charge, $N_{e}$ the electron number
density, $v_{e}$ the electron drift velocity under the electric field $E$, $R$
the gas constant, and $T_{e}$ the electron temperature. Employing this model
requires the inclusion of reactions related to all vibrational states of
species $k$, which significantly increases the number of equations that must
be solved.

Accurately depicting a weakly ionized plasma necessitates accounting for the
electron distribution function, $f_{e}$, which adheres to the Boltzmann
equation in a six-dimensional phase space [6]:

$\frac{\partial f_{e}}{\partial t}+\mathbf{v}_{e}\cdot\nabla
f_{e}-\frac{e}{m}\mathbf{E}\cdot\nabla_{v}f_{e}=\left(\frac{\partial
f_{e}}{\partial t}\right)_{\text{coll}},$ (9)

where $\mathbf{v}_{e}$ represents electron velocity, $e$ the elementary
charge, $m$ the electron mass, $\mathbf{E}$ the electric field, and
$\nabla_{v}$ the velocity gradient operator. The term on the right-hand side
quantifies $f_{e}$’s collision-induced rate of change. Employing the widely
used two-term approximation allows for the decomposition of equation (9), with
components of $f_{e}$ renormalized as probability distribution functions. The
isotropic part, denoted $F_{0}$, serves as the electron energy EEDF [6, 44].
Reaction rate coefficients for collisions $\mathcal{P}$ are derived from
$F_{0}$ and cross-section $\sigma_{j}$ by integrating over electron energy
$\epsilon$:

$k_{j}=\sqrt{2e/m}\int_{0}^{\infty}{\epsilon\sigma_{j}F_{0}}d\epsilon,\quad
j\in\mathcal{P}.$ (10)

Furthermore, $F_{0}$ facilitates the calculation of electron temperature:

$T_{e}=\frac{2}{3}\int_{0}^{\infty}{\epsilon^{3/2}F_{0}}d\epsilon,$ (11)

which is commonly featured in the rate coefficient expressions for
recombination reactions involving electrons.

###  2.2 Heat loss model for pin-pin NRP discharges

NRP discharges produced by pin-pin electrodes typically result in ultrafast
gas heating within an approximately constant volume, generating a weak shock
wave and thermal energy loss due to gas expansion [45, 46]. Such hydrodynamic
phenomena cannot be directly resolved within the limitations of an 0D
framework, leading to potential overestimations of temperature increase if
heat loss mechanisms are neglected. This discrepancy may remain negligible for
scenarios involving low plasma power; however, it becomes increasingly
significant with higher energy deposition and a greater number of nanosecond
pulses. To enhance ChemPlasKin’s ability to accurately model temperature
dynamics, particularly for applications investigating plasma-assisted ignition
delay times, we propose the following heat loss model.

During the nanosecond pulse, the gas within the discharge kernel undergoes a
rise in temperature and pressure under constant volume conditions. Once the
plasma energy input ceases, the discharge kernel experiences an isentropic
expansion from state 1 to state 2, described by:

$\frac{T_{2}}{T_{1}}=\left(\frac{\rho_{2}}{\rho_{1}}\right)^{\gamma-1}=\left(\frac{p_{2}}{p_{1}}\right)^{\frac{\gamma-1}{\gamma}},$
(12)

where $\gamma$ represents the heat capacity ratio. The propagation of the
resulting weak shock wave, occurring on the microsecond timescale, rapidly
equalizes the kernel pressure back to the ambient pressure $p_{2}$. For
practical code implementation, the temperature and density at the pulse’s
conclusion are adjusted to reflect $T_{2}$ and $\rho_{2}$ values within a
single timestep, for given initial ($p_{1}$) and final ($p_{2}$) pressures.
This model effectively modifies the gas temperature equation between constant
volume and constant pressure conditions to mimic the isentropic process, a
flexibility not typically offered by previous solvers listed in Table 1.

The heat loss model is augmented to account for radial thermal conduction from
the hot discharge kernel to the ambient environment, described by:

$\dot{q}_{loss}=-\lambda\frac{dT}{dr}\frac{A_{\text{dis}}}{V_{\text{dis}}}\approx-\lambda
C_{0}\frac{T-T_{0}}{r_{\text{dis}}^{2}},$ (13)

where $\dot{q}_{loss}$ represents the heat loss in $\mathrm{W/m^{3}}$,
$\lambda$ the thermal conductivity, and $r_{\text{dis}}$ the effective radius
of the cylindrical discharge characterized by surface area $A_{\text{dis}}$
and volume $V_{\text{dis}}$. Here, $C_{0}$ is a dimensionless empirical
coefficient, typically set to 1.0, and $T_{0}$ denotes the ambient
temperature. To simulate gas temperature dynamics during intervals without
plasma energy input, the following equation is actually solved:

$\rho
c_{p}\frac{dT_{\text{gas}}}{dt}=-\sum_{k=1}^{N}\dot{\omega}_{k}h_{k}W_{k}+\sum_{k}\dot{R}_{\text{VT}}^{k}-\lambda
C_{0}\frac{T-T_{0}}{r_{\text{dis}}^{2}}.$ (14)

###  2.3 Code development

####  2.3.1 Electron Boltzmann equation (EBE) solver

To enhance performance in simulations, ChemPlasKin incorporates an EBE solver,
named CppBOLOS, directly integrated into the Cantera library at the source
level. This integration represents a central feature of ChemPlasKin’s
functionality. Among the available EBE solvers within the plasma community,
none provides an open-source C++ implementation that employs the efficient
two-term expansion method akin to Bolsig+ [6]. For instance, while MultiBolt
[47] offers a C++ open-source version (v3.x), its multi-term model tends to be
computationally intensive and not essential for a gas-plasma kinetics solver.
Addressing this need, we developed a C++ version of BOLOS [12], an open-source
Python solver that follows the algorithms specified in [6]. Now referred to as
CppBOLOS, this solver’s development process was significantly assisted by
GPT-4 [48]. It is seamlessly integrated with Cantera, enabling dynamic
solutions of the electron Boltzmann equation and real-time updates of the EEDF
based on the temporal evolution of gas temperature, mixture composition, and
reduced electric field. CppBOLOS accepts input cross-section data in the
popular LXCat format [49], which ensures broad compatibility and simplifies
data integration. The validation of CppBOLOS is briefly presented in Section
3.1.

####  2.3.2 Code architecture

Figure 1:  Code architecture of ChemPlasKin.

Figure 1 illustrates the code architecture of ChemPlasKin, highlighting its
streamlined approach to handling both plasma and neutral gas kinetics. Key to
its design is the use of a unified YAML-type input file for all thermal
properties and reactions. This choice facilitates ease of use and integration
with existing datasets. ChemPlasKin extends Cantera’s reaction module to
register inelastic electron-impact reactions classified as type of
‘Boltzmann’, as shown in the examples of Listing LABEL:lst:yaml_example_1:

⬇

1# Example 1: Vibrational excitation

2- equation: Electron + N2 => Electron + N2

3 type: Boltzmann

4 process: N2 -> N2(v1)

5 duplicate: true

6 energy_transfer: {evib_N2: 0.291_eV}

7

8# Example 2: Electronic excitation

9- equation: Electron + O2 => Electron + O2(a1)

10 type: Boltzmann

11 process: O2 -> O2(a1)

12 energy_transfer: {e_th: 0.977_eV}

Listing 1: Examples of YAML configurations for Boltzmann type reactions.
Vibrational states of $\mathrm{N_{2}}$ are agrragated and the vibrational
energy is stored in evib_N2.

To accommodate the diverse range of rate constant expressions encountered in
plasma kinetics, ChemPlasKin employs muParser library [50] for parsing any
complex mathematical expressions not natively supported by Cantera. Listing
LABEL:lst:yaml_example_2 shows an example of $\mathrm{e^{-}+N_{4}^{+}}$
recombination reaction, where the reaction rate is calculated as a function of
electron temperature. The energy_transfer entry values enforce zero gas
heating and the corresponding energy (3.8 eV) is stored in
$e_{\text{vib}}^{\mathrm{N_{2}}}$.

⬇

1# Example 3: Recombination

2- equation: Electron + N4+ => N2 + N2(C3)

3 type: PlasmaCustomExpr

4 energy_transfer: {e_th: -3.8_eV, evib_N2: 3.8_eV} # No gas heating

5 rateExpr: {A: 2.086e+19, Expr: Te^-0.5}

Listing 2: Example of YAML configurations for PlasmaCustomExpr type reactions.

For added convenience, an optional parser tool, `parsePlasKin`, is available
to convert ZDPlaKin input mechanism files into the human-readable YAML format
automatically. This feature is particularly beneficial given the prevalence of
ZDPlaKin mechanisms published as supplementary materials in the field.
Configuration parameters are set through the text files `controlDict` and
`chemPlasProperties`, ensuring straightforward operation.

The core of ChemPlasKin’s simulation capabilities is encapsulated in
`plasmaReactor.h`, where a single ordinary differential equation (ODE)
system—encompassing equations (1, 2a, 2b 4)—is constructed and integrated
within the main time loop of the master code. The ODEs are solved using the
CVODE solver from the SUNDIALS suite [51, 52], which is efficiently utilized
by Cantera for high-precision integration.

####  2.3.3 Code efficiency

ChemPlasKin’s unified ODE system offers distinct advantages over operator-
splitting techniques, by eliminating splitting errors and enabling much larger
timesteps in the main time loop. Although it is suggested that splitting the
plasma kinetics could expedite integration due to the higher stiffness during
nanosecond pulses, such benefits become marginal during the much longer pulse
intervals where plasma stiffness diminishes. Moreover, the plasma kinetics
component of ChemPlasKin, akin to that solved by ZDPlasKin, encompasses a
broad spectrum of timescales beyond just ultrafast electron-impact reactions.
Last but not least, switching between ZDPlasKin and CHEMKIN/Cantera requires
reinitialization of their respective ODE solvers at every step, adding extra
computational overhead to the solvers in the initial transient to find
suitable integration timesteps. Consequently, ChemPlasKin’s unified approach
not only simplifies the integration process but also enhances overall
simulation efficiency compared to the operator-splitting method. To
demonstrate this, we have also built a ZDPlasKin-Cantera solver consistent
with previous counterparts and compared it with ChemPlasKin, as detailed in
Section 3.2.

##  3 Code Validation

This section begins with a brief validation of CppBOLOS against Bolsig+,
establishing the groundwork for subsequent evaluations. The subsequent five
subsections are carefully structured to assess ChemPlasKin’s performance
across diverse simulation scenarios. We present comparisons of simulation
results between literature sources, our equivalent Cantera-ZDPlasKin solver,
and ChemPlasKin in Section 3.2. The capability of ChemPlasKin to predict
ultrafast gas heating and radical species production is assessed in Section
3.3. The slow gas heating model is validated through experimental data in
Section 3.4, while a practical application case of fuel reforming is examined
in Section 3.5. Finally, the efficacy of the heat loss model is tested in
Section 3.6. For the simulations discussed in Sections 3.3, 3.4, and 3.6, we
utilize a rigorously validated detailed PAC mechanism involving 100 species
and 964 reactions for a methane-air mixture, as developed by [11].

###  3.1 CppBOLOS

Figure 2 illustrates the comparison between CppBOLOS and Bolsig+ for mean
electron energy and ionization rate in pure $\mathrm{N_{2}}$ across a wide
range of reduced electric fields ($E/N$). Additionally, relative differences
compared to Bolsig+ are plotted, demonstrating satisfactory consistency except
at very low $E/N$ values. The large relative differences in ionization rate at
low $E/N$ are considered acceptable because the infinitesimally small absolute
values have a negligible impact on the accurate prediction of the plasma
kinetics system’s evolution.

Figure 2:  Comparison between CppBOLOS and Bolsig+ for pure $\mathrm{N_{2}}$
at 300 K, using cross-section data from Phelps [49].

###  3.2 Plasma assisted $\mathbf{H_{2}/O_{2}/He}$ ignition

Mao et al. [18] conducted numerical simulations of $\mathrm{H_{2}/O_{2}/He}$
ignition assisted by hybrid NRP and DC discharges using the ZDPlasKin-CHEMKIN
solver developed by Lefkowitz et al. [9]. As no openly available PAC solver
exists for use as a benchmark, we integrated ZDPlasKin with Cantera to
evaluate the computational efficiency of the popular operator-splitting
technique, adhering to the strategy described in [9]. The Cantera-ZDPlasKin
solver includes a Python wrapper of ZDPlasKin that facilitates efficient data
communications between the two codes. The timestep in the main time loop is
dynamically adjusted to accommodate the nanosecond timescale and exponential
growth of electrons in NRP discharges, then gradually relaxed to $10^{-7}$ s
during the pulse intervals, balancing accuracy with the computational cost of
the operator-splitting method.

The mixture has an initial composition of
$\mathrm{H_{2}:O_{2}:He=0.1667:0.0883:0.75}$, and the adiabatic system is
maintained at atmospheric pressure. The reduced electric field ($E/N$) for NRP
discharges is set at 100 Td (1 Td = $\mathrm{10^{-21}Vm^{2}}$), with a
frequency of 30 kHz. DC discharges at an $E/N$ of 20 Td are applied between
the nanosecond pulses. The deposited plasma energy is fixed at 0.1
$\mathrm{mJ/cm^{3}}$ per pulse. To ensure consistency, ChemPlasKin operates
with the Flitti & Pancheshnyi [43] model for gas temperature. The plasma
kinetic mechanism used by [18], originally formatted for ZDPlasKin, has been
automatically converted into YAML format compatible with ChemPlasKin.

Figure 3 compares these three solvers for their predictions of ignition delay
times (IDT) for NRP and NRP/DC hybrid discharges assisted ignition across
various initial temperatures. IDT is defined as the point of maximum
temperature gradient during the pulse intervals in our simulations. Good
agreement is achieved among the three solvers for both discharge types.

Figure 3:  Ignition delay times (IDT) for NRP discharge assisted ignition and
NRP/DC hybrid discharge assisted ignition. Results from ChemPlasKin are
compared against those from Mao et al. [18] and our Cantera-ZDPlasKin solver.

The predictive capabilities of the solvers for electron and
$\mathrm{O_{2}(a^{1}\Delta_{g})}$ production over the first ten pulses in the
hybrid NRP/DC discharge scenario are displayed in Figure 4. These results
underscore the performance consistency of our solvers, particularly Cantera-
ZDPlasKin, with those reported by Mao et al. [18].

Figure 4:  Temporal evolution of electron and
$\mathrm{O_{2}(a^{1}\Delta_{g})}$ mole fractions during the first ten pulses
of a hybrid discharge at 400 K.

A comparison of computational efficiency between ChemPlasKin and Cantera-
ZDPlasKin was carried out for NRP discharge assisted ignition cases, as shown
in Figure 5. ChemPlasKin achieved a three-fold speed-up. Note that in Figure
5, ChemPlasKin employed the same dynamic timestep settings in the main loop to
ensure a fair comparison. In actual applications, ChemPlasKin can utilize much
larger outer timesteps, as it is not constrained by operator-splitting errors.

Figure 5:  Computational cost (clock time) comparison between ChemPlasKin and
Cantera-ZDPlasKin corresponding to the NRP assisted ignition depicted in
Figure 3. Single-thread execution at 2.3 GHz. Profiling results confirm that
there is no data communication bottleneck in the Cantera-ZDPlasKin solver.

Note that many operator-splitting solvers listed in Table 1, which operate
based on kinetic mechanisms developed in various studies, lack rigorous
validation against experimental data at micro time scales, such as radical
production and fast and slow gas heating during a single NRP discharge period.
The following two subsections will present validation of ChemPlasKin that
addresses this gap, similar to the approach taken by Cheng et al. [11].

###  3.3 Spark discharge in air

This test case aligns with Case A from [11], utilizing the experimental setup
described in [53], where NRP spark discharges are generated between pin-pin
electrodes. The gas temperature can be inferred from the rotational
temperature of $\mathrm{N_{2}(C)}$ for a duration of at least 17 ns [53]. The
initial gas mixture temperature is set at 1500 K with a composition of
$\mathrm{N_{2}:O_{2}:O=77.4:18.6:4}$, accounting for the thermal and chemical
effects of preceding pulses. Figure 6 presents a comparison of the gas
temperature evolution and $\mathrm{O(^{3}P)}$ production as predicted by
ChemPlasKin against experimental data.

Figure 6:  Comparison of ChemPlasKin predictions with experimental data on the
temporal evolution of (a) $\mathrm{O(^{3}P)}$ production during and after a
nanosecond pulse and (b) gas temperature, demonstrating good agreement.
Experimental data from [53].

###  3.4 V-T relaxation in air

This validation case aligns with Case C from [11], utilizing the experimental
setup detailed in [54]. The initial conditions of air are set at a temperature
of 300 K and a pressure of 100 Torr. Between the pin-pin electrodes, the
applied $E/N$ is approximately 100 Td, with up to $50\%$ of the plasma energy
contributing to the vibrational excitation of $\mathrm{N_{2}}$. In this
scenario, the model’s prediction of the rotational temperature, represented by
the solid curve shows good agreement with the experimental data, as shown in
Figure 7. Note that the simulation focuses solely on $\mathrm{N_{2}}$ for the
vibrational energy equation, under the assumption of isobaric conditions due
to the minimal impact of fast gas heating in this context. The exclusion of
heat loss from thermal diffusion in the model accounts for the observed
discrepancies between calculated results and experimental measurements on the
millisecond timescale.

Figure 7:  ChemPlasKin predictions of gas temperature from VT-relaxation of
$\mathrm{N_{2}}(v)$ in pure air, showcasing a comparison against experimental
measurements from [54].

###  3.5 Hydrogen oxidation with DBD discharges

This validation case compares ChemPlasKin’s performance against both
experimental and numerical results detailed in [32], which investigates the
kinetics of plasma-assisted oxidation of $\mathrm{H_{2}}$ in an undiluted lean
$\mathrm{H_{2}/O_{2}}$ system. In the referenced study, a preheated gas
mixture is subjected to treatment in a dielectric barrier discharge (DBD)
reactor under constant pressure and temperature conditions. To approximate the
filamentary nature of the discharges, it is assumed that each fluid particle
encounters a total of 288 rectangular pulses at a constant $E/N$.

Figure 8 presents a comparison of the $\mathrm{H_{2}}$ oxidation rates derived
from ChemPlasKin against the experimental and numerical results under varying
gas temperatures. The agreement between ChemPlasKin’s predictions and the
reference results is generally good, with a larger discrepancy at higher
temperatures, potentially attributed to the omission of thermal cracking
residence time in our model. ChemPlasKin completes each data point presented
in Figure 8 within approximately 45 seconds of single-threaded execution on a
2.3 GHz core.

Figure 8:  Comparison of ChemPlasKin simulations with reference data on
plasma-assisted oxidation rates of $\mathrm{H_{2}}$ over a range of gas
temperatures ($T_{\text{gas}}$) for an $\mathrm{H_{2}/O_{2}}$ mixture
($\phi=0.01$) at 10 W of plasma power. Reference data from [32].

###  3.6 Heat loss model

The proposed heat loss model, detailed in Section 3.6, is validated focusing
on the fast gas-heating phenomenon, thereby omitting the need for detailed
plasma kinetics and slow gas heating. Figure 9 illustrates a typical scenario
for applying the heat loss model, aligned with the test case described in
Section 3.3. Here, NRP discharges at 10 kHz are generated between pin-pin
electrodes, with each pulse delivering an energy input of $670\ \mathrm{\mu
J}$. According to the phenomenological plasma model by [4], $20\%$ of the
pulse energy is allocated to fast gas-heating within a 10 ns duration. A
horizontal 2D computational domain is defined, featuring the pulsed heating
source at its center, modeled with an energy density spatial function
$\mathcal{F}=\text{erfc}\left(r^{2}/{6.9\times 10^{-8}}\right)^{2.5}$, where
$r$ denotes the radial distance from the discharge center and an effective
discharge radius $r_{\text{dis}}$ of $225\ \mathrm{\mu m}$ is specified. The
domain employs a non-reflective boundary condition, ensuring the pressure
field near the discharge zone recovers the ambient level within microseconds
[55].

The averaged temperature within $r_{\text{dis}}$ in the 2D domain, solved
using OpenFOAM [56], represents the discharge kernel temperature. In parallel,
a 0D model, simply written as an independent Python code, utilizes the Cantera
library to update thermal conductivity based on temperature changes. After
applying ten heating pulses, Figure 10 displays the temperature comparison
between the 2D and 0D models. The models demonstrate good overall agreement,
with the 0D model’s assumption of instantaneous isentropic expansion capturing
the sharp temperature drop at each pulse’s conclusion, as elaborated in the
zoom-in subplot of Figure 10. Without incorporating the heat loss model, a
conventional 0D chemistry-plasma model would likely predict a step-wise,
monotonic temperature increase for NRP discharges, contrasting with the quasi-
steady state indicated by the 2D model.

Figure 9:  Schematic of a 2D computational domain for simulating ultrafast
heating effects from pin-pin generated NRP discharges using phenomenological
model. Figure 10:  Comparison of temperature evolution within the discharge
kernel across 10 NRP discharges between 2D and 0D phenomenological models. V-T
relaxation heating not included. The red bars highlight the highest
temperature for each nanosecond pulse in the 2D model. $C_{0}=1.0$ is employed
in the 0D model.

Following the standalone examination discussed previously, the impact of the
heat loss model on gas temperature in ChemPlasKin, which incorporates detailed
gas-plasma kinetics of a methane/air mixture [11], is presented in Figure 11.
The discrepancy in temperature evolution with and without the heat loss model
during the first ten pulses clearly demonstrates its potential influence on
ignition delay time calculations. We acknowledge that further validation and
improvement of this model are necessary, ideally through direct comparison
with experimental data or high-fidelity multi-dimensional simulations.

Figure 11:  Temperature evolution in a
$\mathrm{0.0499CH_{4}/0.7505N_{2}/0.1996O_{2}}$ mixture with NRP discharge
assisted ignition over the first ten pulses. Deposited plasma energy is
$\mathrm{10\ mJ/cm^{3}}$ per pulse, with $C_{0}=1.0$.

##  4 Conclusions

This study describes ChemPlasKin, a freeware tool developed for simulating
gas-plasma kinetic processes. ChemPlasKin integrates the electron Boltzmann
equation solver CppBOLOS with the Cantera library, enabling the solution of
neutral gas and plasma kinetics within a unified ODE system. Additionally, a
supplementary heat loss model is proposed to enhance the accuracy of
temperature predictions for nanosecond repetitively pulsed (NRP) discharges in
configurations utilizing pin-pin electrodes.

To evaluate the computational efficiency of ChemPlasKin, we constructed a
Cantera-ZDPlasKin PAC solver using the widely used operator-splitting method.
In the test cases, this configuration achieved a threefold speedup, and we
anticipate even faster performance with a relaxed outer timestep.

The C++ solver has been validated against experimental results across various
aspects and timescales, including ultrafast gas heating and radical
production, slow gas heating from V-T relaxation and fuel reforming involving
hundreds of pulses. ChemPlasKin shows its capability as a versatile tool,
particularly useful for PAC and fuel reforming. It can also be used purely as
a plasma kinetics solver, akin to ZDPlasKin. Its development provides a
resource for researchers seeking to explore gas-plasma kinetics without the
need to invest in custom coding efforts. Continued validation and enhancement
by the user community are encouraged. We also anticipate ChemPlasKin to be
incorporated into CFD codes to enable high-fidelity simulations of fully-
coupled plasma-assisted reacting flows, acknowledging that computational costs
remain a significant consideration.

## Acknowledgments

This research was funded by King Abdullah University of Science and Technology
(KAUST). The authors would like to thank Seunghwan Bang and Dr. Ramses Snoeckx
for useful discussion. The first author also thanks friends Renston, Vijay,
and Alessandro for the great tunes and good times that helped lighten the
load.

## References

  * [1] Yiguang Ju and Wenting Sun.  Plasma assisted combustion: Dynamics and chemistry.  Progress in Energy and Combustion Science, 48:21–83, June 2015.
  * [2] Deanna A. Lacoste.  Flames with plasmas.  Proceedings of the Combustion Institute, 39(4):5405–5428, January 2023.
  * [3] Yaolin Wang, Yanzhen Chen, Jonathan Harding, Hongyuan He, Annemie Bogaerts, and Xin Tu.  Catalyst-free single-step plasma reforming of CH4 and CO2 to higher value oxygenates under ambient conditions.  Chemical Engineering Journal, 450:137860, December 2022.
  * [4] Maria Luis Gracio Bilro Castela.  Direct Numerical Simulations of plasma-assisted ignition in quiescent and turbulent flow conditions.  phdthesis, Université Paris-Saclay, May 2016.
  * [5] Nicolas Barléon, Lionel Cheng, Bénédicte Cuenot, and Olivier Vermorel.  A phenomenological model for plasma-assisted combustion with NRP discharges in methane-air mixtures: PACMIND.  Combustion and Flame, 253:112794, July 2023.
  * [6] G. J. M. Hagelaar and L. C. Pitchford.  Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models.  Plasma Sources Science and Technology, 14(4):722, October 2005.
  * [7] S Pancheshnyi, B Eismann, GJM Hagelaar, and LC Pitchford.  Computer code ZDPlasKin.  University of Toulouse, LAPLACE, CNRS-UPS-INP, Toulouse, France, 2008.
  * [8] R. J. Kee, F. M. Rupley, and J. A. Miller.  Chemkin-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics.  Technical Report SAND-89-8009, Sandia National Lab. (SNL-CA), Livermore, CA (United States), September 1989.
  * [9] Joseph K Lefkowitz, Peng Guo, Aric Rousso, and Yiguang Ju.  Species and temperature measurements of methane oxidation in a nanosecond repetitively pulsed discharge.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 373(2048):20140333, August 2015.  Publisher: Royal Society.
  * [10] David G. Goodwin, Harry K. Moffat, Ingmar Schoegl, Raymond L. Speth, and Bryan W. Weber.  Cantera: An Object-oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes, 2022.
  * [11] L. Cheng, N. Barleon, B. Cuenot, O. Vermorel, and A. Bourdon.  Plasma assisted combustion of methane-air mixtures: Validation and reduction.  Combustion and Flame, 240:111990, June 2022.
  * [12] aluque.  aluque/bolos, July 2023.  original-date: 2014-06-10T08:17:00Z.
  * [13] Igor V. Adamovich, Ting Li, and Walter R. Lempert.  Kinetic mechanism of molecular energy transfer and chemical reactions in low-temperature air-fuel plasmas.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 373(2048):20140336, August 2015.  Publisher: Royal Society.
  * [14] Anthony C. DeFilippo and Jyh-Yuan Chen.  Modeling plasma-assisted methane–air ignition using pre-calculated electron impact reaction rates.  Combustion and Flame, 172:38–48, October 2016.
  * [15] Nicholas Deak, Aurélie Bellemans, and Fabrizio Bisetti.  Plasma-assisted ignition of methane/air and ethylene/air mixtures: Efficiency at low and high pressures.  Proceedings of the Combustion Institute, 38(4):6551–6558, January 2021.
  * [16] T. Hazenberg, J. van Dijk, and J. A. van Oijen.  Chemical flux analysis of low-temperature plasma-enhanced oxidation of methane and hydrogen in argon.  Combustion and Flame, 257:113037, November 2023.
  * [17] Xingqian Mao and Qi Chen.  Effects of Vibrational Excitation on Nanosecond Discharge Enhanced Methane–Air Ignition.  AIAA Journal, 56(11):4312–4320, 2018.  Publisher: American Institute of Aeronautics and Astronautics _eprint: https://doi.org/10.2514/1.J057304.
  * [18] Xingqian Mao, Qi Chen, Aric C. Rousso, Timothy Y. Chen, and Yiguang Ju.  Effects of controlled non-equilibrium excitation on H2/O2/He ignition using a hybrid repetitive nanosecond and DC discharge.  Combustion and Flame, 206:522–535, August 2019.
  * [19] Xingqian Mao, Aric Rousso, Qi Chen, and Yiguang Ju.  Numerical modeling of ignition enhancement of CH4/O2/He mixtures using a hybrid repetitive nanosecond and DC discharge.  Proceedings of the Combustion Institute, 37(4):5545–5552, January 2019.
  * [20] Yuan Wang, Peng Guo, Haitao Chen, and Zheng Chen.  Numerical modeling of ignition enhancement by repetitive nanosecond discharge in a hydrogen/air mixture II: forced ignition.  Journal of Physics D: Applied Physics, 54(6):065502, November 2020.  Publisher: IOP Publishing.
  * [21] Hongtao Zhong, Xingqian Mao, Aric C Rousso, Charles L Patrick, Chao Yan, Wenbin Xu, Qi Chen, Gerard Wysocki, and Yiguang Ju.  Kinetic study of plasma-assisted n-dodecane/O2/N2 pyrolysis and oxidation in a nanosecond-pulsed discharge.  Proceedings of the Combustion Institute, 38(4):6521–6531, January 2021.
  * [22] Galia Faingold and Joseph K. Lefkowitz.  A numerical investigation of NH3/O2/He ignition limits in a non-thermal plasma.  Proceedings of the Combustion Institute, 38(4):6661–6669, January 2021.
  * [23] Galia Faingold, Omer Kalitzky, and Joseph K. Lefkowitz.  Plasma reforming for enhanced ammonia-air ignition: A numerical study.  Fuel Communications, 12:100070, September 2022.
  * [24] Hongtao Zhong, Xingqian Mao, Ning Liu, Ziyu Wang, Timothy Ombrello, and Yiguang Ju.  Understanding non-equilibrium N2O/NOx chemistry in plasma-assisted low-temperature NH3 oxidation.  Combustion and Flame, 256:112948, October 2023.
  * [25] Xingqian Mao, Hongtao Zhong, Ning Liu, Ziyu Wang, and Yiguang Ju.  Ignition enhancement and NOx formation of NH3/air mixtures by non-equilibrium plasma discharge.  Combustion and Flame, 259:113140, January 2024.
  * [26] Nan Liu, Qi Chen, Xianwu Jiang, Jie Chen, and Lidong Zhang.  Kinetic roles of excited state species in the oxidation of n-pentane/air assisted by nanosecond-pulsed discharge.  Combustion and Flame, 259:113142, January 2024.
  * [27] Vyaas Gururajan and Fokion N. Egolfopoulos.  Transient plasma effects on the autoignition of DME/O2/Ar and C3H8/O2/Ar mixtures.  Proceedings of the Combustion Institute, 36(3):4165–4174, January 2017.
  * [28] Ruzheng Zhang, Handong Liao, Jiuzhong Yang, and Bin Yang.  Exploring chemical kinetics of plasma assisted oxidation of dimethyl ether (DME).  Combustion and Flame, 225:388–394, March 2021.
  * [29] Taaresh Sanjeev Taneja and Suo Yang.  Comparing Low-Mach and Fully-Compressible CFD Solvers for Phenomenological Modeling of Nanosecond Pulsed Plasma Discharges with and without Turbulence.  In AIAA SCITECH 2022 Forum. American Institute of Aeronautics and Astronautics, 2021.  _eprint: https://arc.aiaa.org/doi/pdf/10.2514/6.2022-0976.
  * [30] Taaresh Sanjeev Taneja, Praise Noah Johnson, and Suo Yang.  Nanosecond pulsed plasma assisted combustion of ammonia-air mixtures: Effects on ignition delays and NOx emission.  Combustion and Flame, 245:112327, November 2022.
  * [31] Praise N. Johnson, Taaresh S. Taneja, and Suo Yang.  Plasma-based global pathway analysis to understand the chemical kinetics of plasma-assisted combustion and fuel reforming.  Combustion and Flame, 255:112927, September 2023.
  * [32] Ramses Snoeckx, Daeyoung Jun, Bok Jik Lee, and Min Suk Cha.  Kinetic study of plasma assisted oxidation of H2 for an undiluted lean mixture.  Combustion and Flame, 242:112205, August 2022.
  * [33] Ramses Snoeckx and Min Suk Cha.  Kinetic study of plasma assisted oxidation of H2 for an undiluted rich mixture.  Combustion and Flame, 250:112638, April 2023.
  * [34] Seunghwan Bang, Ramses Snoeckx, and Min Suk Cha.  Kinetic Study for Plasma Assisted Cracking of NH3: Approaches and Challenges.  The Journal of Physical Chemistry A, 127(5):1271–1282, February 2023.  Publisher: American Chemical Society.
  * [35] Seunghwan Bang, Ramses Snoeckx, and Min Suk Cha.  Temperature-Dependent Kinetics of Ozone Production in Oxygen Discharges.  Plasma Chemistry and Plasma Processing, August 2023.
  * [36] Mohammad Shahsavari, Alexander A. Konnov, Agustin Valera-Medina, and Mehdi Jangi.  On nanosecond plasma-assisted ammonia combustion: Effects of pulse and mixture properties.  Combustion and Flame, 245:112368, November 2022.
  * [37] Colin A. Pavan and Carmen Guerra-Garcia.  Modeling Flame Speed Modification by Nanosecond Pulsed Discharges to Inform Experimental Design.  In AIAA SCITECH 2023 Forum, AIAA SciTech Forum. American Institute of Aeronautics and Astronautics, January 2023.
  * [38] Yifan Qiu, Yifei Zhu, Yun Wu, Ningqiu Zhao, Zhenyang Li, Mai Hao, Boya Zhang, and Di Pan.  Numerical investigation of the hybrid pulse–DC plasma assisted ignition and NOx emission of NH3/N2/O2 mixture.  Combustion and Flame, 258:113078, December 2023.
  * [39] Ziying Xin, Zhencao Zheng, Yong Hu, Ao Sun, Feiyang Zhao, and Wenbin Yu.  Numerical modeling of plasma assisted ignition of CH4/O2/He mixture by the nanosecond repetitive pulsed surface dielectric barrier discharge.  Fuel, 357:129975, February 2024.
  * [40] Roger C. Millikan and Donald R. White.  Systematics of Vibrational Relaxation.  The Journal of Chemical Physics, 39(12):3209–3213, December 1963.
  * [41] S. J. Colgan and B. P. Levitt.  Vibrational relaxation of nitrogen by various collision partners at 2500°K.  Transactions of the Faraday Society, 63(0):2898–2905, 1967.  Publisher: Royal Society of Chemistry.
  * [42] N. A. Popov.  Fast gas heating in a nitrogen–oxygen discharge plasma: I. Kinetic mechanism.  Journal of Physics D: Applied Physics, 44(28):285201, June 2011.
  * [43] A. Flitti and S. Pancheshnyi.  Gas heating in fast pulsed discharges in N2–O2 mixtures.  The European Physical Journal Applied Physics, 45(2):21001, February 2009.  Number: 2 Publisher: EDP Sciences.
  * [44] L. L. Alves, A. Bogaerts, V. Guerra, and M. M. Turner.  Foundations of modelling of nonequilibrium low-temperature plasmas.  Plasma Sources Science and Technology, 27(2):023002, February 2018.  Publisher: IOP Publishing.
  * [45] Da A. Xu, Deanna A. Lacoste, and Christophe O. Laux.  Schlieren Imaging of Shock-Wave Formation Induced by Ultrafast Heating of a Nanosecond Repetitively Pulsed Discharge in Air.  IEEE Transactions on Plasma Science, 42(10):2350–2351, October 2014.  Conference Name: IEEE Transactions on Plasma Science.
  * [46] Fabien Tholin and Anne Bourdon.  Simulation of the hydrodynamic expansion following a nanosecond pulsed spark discharge in air at atmospheric pressure.  Journal of Physics D: Applied Physics, 46(36):365205, August 2013.  Publisher: IOP Publishing.
  * [47] J. Stephens.  A multi-term Boltzmann equation benchmark of electron-argon cross-sections for use in low temperature plasma models.  Journal of Physics D: Applied Physics, 51(12):125203, March 2018.  Publisher: IOP Publishing.
  * [48] GPT-4.
  * [49] S. Pancheshnyi, S. Biagi, M. Bordage, G. J. M. Hagelaar, W. Morgan, Arthur Phelps, and L. Pitchford.  The LXCat project: Electron scattering cross sections and swarm parameters for low temperature plasma modeling.  Chemical Physics, 398:148–153, April 2012.
  * [50] muparser - fast math parser library.
  * [51] Alan C. Hindmarsh, Peter N. Brown, Keith E. Grant, Steven L. Lee, Radu Serban, Dan E. Shumaker, and Carol S. Woodward.  SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers.  ACM Transactions on Mathematical Software, 31(3):363–396, September 2005.
  * [52] David J. Gardner, Daniel R. Reynolds, Carol S. Woodward, and Cody J. Balos.  Enabling New Flexibility in the SUNDIALS Suite of Nonlinear and Differential/Algebraic Equation Solvers.  ACM Transactions on Mathematical Software, 48(3):31:1–31:24, September 2022.
  * [53] D. L. Rusterholtz, D. A. Lacoste, G. D. Stancu, D. Z. Pai, and C. O. Laux.  Ultrafast heating and oxygen dissociation in atmospheric pressure air by nanosecond repetitively pulsed discharges.  Journal of Physics D: Applied Physics, 46(46):464010, October 2013.  Publisher: IOP Publishing.
  * [54] A. Montello, Z. Yin, D. Burnette, I. V. Adamovich, and W. R. Lempert.  Picosecond CARS measurements of nitrogen vibrational loading and rotational/translational temperature in non-equilibrium discharges.  Journal of Physics D: Applied Physics, 46(46):464002, October 2013.  Publisher: IOP Publishing.
  * [55] Xiao Shao, Narjisse Kabbaj, Deanna A. Lacoste, and Hong G. Im.  A computational study of a laminar methane–air flame assisted by nanosecond repetitively pulsed discharges.  Journal of Physics D: Applied Physics, 57(20):205201, February 2024.  Publisher: IOP Publishing.
  * [56] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  Computer in Physics, 12(6):620–631, November 1998.