# Coupled Kinetic and Electrical Simulations of Methane Plasmalysis by Nanosecond Discharges

Pierre-Antoine Goutier, Sara Delahaie, Sean Mcguire, Nicolas Minesi,
Christophe Laux, Erwan Pannier

► To cite this version:

Pierre-Antoine Goutier, Sara Delahaie, Sean Mcguire, Nicolas Minesi, Christophe Laux, et al. Coupled Kinetic and Electrical Simulations of Methane Plasmalysis by Nanosecond Discharges. AIAA SCITECH 2026 Forum, Jan 2026, Orlando, United States. 10.2514/6.2026-1424. hal-05483196

HAL Id: hal-05483196

https://hal.science/hal-05483196v1

Submitted on 29 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Coupled kinetic and electrical simulations of methane plasmalysis by nanosecond discharges

P.-A. Goutier<fnref n="1" />, S. Delahaie<fnref n="2" />, S. D. McGuire<fnref n="3" />, N. Q. Minesi<fnref n="4" />, C. O. Laux<fnref n="5" />

Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, Gif-sur-Yvette, 91190, France

E. Pannier<fnref n="6" />

Spark Cleantech, Gif-sur-Yvette, 91190, France

We model the conversion of methane (CH₄) towards hydrogen (H₂) and acetylene (C₂H₂) at atmospheric pressure via Nanosecond Repetitively Pulsed (NRP) discharges. The chemical kinetic mechanism includes hydrocarbon species, CₓHᵧ, up to ethane (C₂H₆), along with their corresponding ions. The reaction rate constants available in the literature are reviewed and compared to build a kinetic model valid in the temperature range of interest for plasmalysis of methane. A zero-dimensional kinetic model is then used in a closed two-temperature system to solve the temporal evolution of species concentrations, the gas (heavy particle) energy, and the electron energy. This model also includes the interactions between the pulser, the transmission line (coaxial cable), and the load (plasma). Furthermore, pressure relaxation effects are considered. The simulated electron number density reproduces the trends observed in experimental data within one order of magnitude. The simulation provides the plasma resistance, which can then be used to design optimal pulsing strategies to increase the energy efficiency.

## I. Introduction

The transition from fossil fuels to hydrogen could reduce the CO₂ emissions from combustion processes used in high-temperature industries (such as cement plants, steel mills, and glass factories). However, as of today, the H₂ molecule is mainly synthesized by steam methane reforming (CH₄ + 2H₂O → CO₂ + 4H₂), which emits approximately 10 tons of CO₂ per ton of hydrogen produced [1]. Water electrolysis (2H₂O → O₂ + 2H₂) is an alternative process that requires high energy input to break the molecular bonds of water. There is a third pathway, the so-called plasmalysis of methane, wherein methane is converted to hydrogen (CH₄ → C(s) + 2H₂) using an electric discharge to provide the necessary energy input. Plasmalysis is less electricity-intensive than water electrolysis (thanks to a reaction enthalpy of 5 kWh/kgH₂ for plasmalysis vs 40 kWh/kgH₂ for water electrolysis), produces a valuable co-product (solid carbon), and emits no greenhouse gases if the electrical power is produced in a CO₂-free manner.

Plasmalysis can be performed using a DC arc discharge to increase the gas temperature and achieve thermal conversion [2]. Warm non-equilibrium Nanosecond Repetitively Pulsed (NRP) discharges can also be employed to initiate plasmalysis kinetics at reduced temperature, thereby improving the thermal efficiency and/or the properties of the solid carbon product [3]. Multiple regimes exist in NRP discharges (corona, glow, non-equilibrium spark [4], and

<footnote n="1">¹ PhD student, EM2C Laboratory and Spark Cleantech</footnote>

<footnote n="2">² PhD student, EM2C Laboratory and Spark Cleantech</footnote>

<footnote n="3">³ Associate Professor, CentraleSupélec, EM2C Laboratory</footnote>

<footnote n="4">⁴ Associate Professor, CentraleSupélec, EM2C Laboratory</footnote>

<footnote n="5">⁵ Professor, CentraleSupélec, EM2C Laboratory, AIAA Fellow</footnote>

<footnote n="6">⁶ Chief Technical Officer, Spark Cleantech</footnote>
thermal spark [5]), and their chemical, thermal, and flow dynamic effects must be understood to design efficient
reactors and pulsers.

This work focuses on the time-varying electrical resistance of the plasma generated by nanosecond discharges, in the (non-equilibrium and thermal) spark regime. The ability to predict this resistance is crucial for optimizing electrical power injection. Section II describes our working hypothesis, the model, and the derivation of the balance equations for the two-temperature homogeneous discharge. We then present the chemical species considered in the model, their thermodynamic properties, and the kinetic mechanism. In Section III, the simulation results - particularly the electron density - are compared to experimental data. Appendix A describes the method to calculate the electrical conductivity, which requires momentum-transfer collision frequencies shown in Appendix B.

II. Model description

In Section A, the principal working hypotheses of the model are listed and discussed. In Section B, the species mass balance, the electron energy balance, and the gas energy balance are derived. The numerical implementation is also presented. Finally, Section C describes the thermodynamic properties and kinetic reactions used in the model.

A. System definition and hypothesis

The system studied consists of the charged volume created by the electric discharge, assumed to be a cylinder of radius $r_d$ and length $L_{gap}$ (the interelectrode gap distance), as represented in Fig. 1. Whereas $L_{gap}$ is constant, $r_d$ is allowed to vary with time during the simulations. Since the electrical resistance is the parameter of interest, it is sufficient to consider only the first hundreds of nanoseconds because most electrical energy is deposited during this interval (as seen later in the companion presentation of experimental data by Delahaie *et al.* [6]). To simplify computations, the system is assumed to be spatially uniform and closed (no mass entering or leaving). The electric field $E_p(t)$ is assumed to be uniform and is calculated by dividing the applied voltage across the plasma $V_p(t)$ (neglecting cathode and anode falls) by the interelectrode gap $L_{gap}$. Moreover, we only consider gas-phase chemical reactions (surface reactions are neglected). These simplifications allow a 0-D kinetic model to be used.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Schematic representation of the system studied. (left) Geometry of the discharge at $t = 0$ ns. (right) Discharge volume has grown due to pressure increase during the pulse, followed by pressure relaxation after the pulse. The mass is conserved within the volume throughout the entire process.</figcaption></figure>

The temporal evolution of a nanosecond spark discharge at atmospheric pressure can be decomposed into several phases that are simplified to correspond to the 0-D nature of our simulation. These phases include the streamer propagation, the pressure relaxation (and the resultant shock wave), the cooling by radiation and diffusion, and finally the recirculation (gas ejection, in torus shape for pin-pin electrodes [7]). In this work, we assume that the streamer phase does not occur. This assumption holds if the initial electron density is above a certain critical density, i.e. $n_{e,0} > 10^8$ cm$^{-3}$, as was observed by Pai *et al.* [8] for atmospheric pressure discharges in preheated air repeated at high frequency. In their case, the electron number density increased uniformly via avalanche breakdown in the volume. We will therefore assume an initial electron number density of $10^{13}$ cm$^{-3}$.

During the discharge, the temperature and electron number density increase, resulting in a corresponding pressure increase. For a non-equilibrium spark, the pressure and temperature increase up to a few bars and ~1,000 K [8]. For a thermal spark, the temperature can increase by up to 50,000 K [5], which corresponds to a pressure of several hundred bars. Table 1 summarizes the characteristic times of post-discharge phenomena for various discharges. For the case of a non-equilibrium spark (electron density below $10^{17}$ cm$^{-3}$), pressure relaxation can be neglected up to hundreds
of nanoseconds. However, in the case of a thermal spark, pressure relaxation occurs at much shorter times - on the order of tens of nanoseconds. Thus, for the simulations presented in this work, the isochoric assumption may not be used. Pressure relaxation will be accounted for using an expansion term. Radiation, diffusion, and recirculation are neglected, as their characteristic times are much longer than the time window studied.

<table><caption>Table 1. Characteristic times reported for the main post-discharge phenomena.</caption><thead><tr><th>Ref.</th><th>Conditions</th><th>τ<sub>rarefaction</sub></th><th>τ<sub>radiation</sub></th><th>τ<sub>diffusion</sub></th><th>τ<sub>recirculation</sub></th></tr></thead><tbody><tr><td rowspan="2">Maillard [9]</td><td>Non-equilibrium Spark, CO<sub>2</sub></td><td>0.5–1 μs</td><td>&gt; 500 ns</td><td>10 μs – 10 ms</td><td>-</td></tr><tr><td>Thermal Spark, CO<sub>2</sub></td><td>10 ns</td><td>250 ns (full spectra)</td><td>1 – 10 μs</td><td>-</td></tr><tr><td>Minesi et al. [5]</td><td>Thermal Spark, Air</td><td>14 ns</td><td>&gt; 60 ns (black body)</td><td>12 μs</td><td>&gt; 1 μs</td></tr><tr><td>Roger et al. [7]</td><td>Thermal Spark, Air</td><td>10 ns – 1 μs</td><td>-</td><td>10 – 100 μs</td><td>10 – 100 μs</td></tr></tbody></table>

Electrons and heavy species are assumed to follow a Maxwell-Boltzmann distribution at their respective temperatures ($T_e$ and $T_g$). This approximation will be relaxed in future work, because the electron energy distribution function $f_e$ is usually not Maxwellian in plasma discharges. The temperatures $T_e$ and $T_g$ correspond to the translational temperature of electrons and of heavy species (atoms, molecules, ions), respectively. Finally, all chemical species are assumed to be in their ground states - vibrational and electronic excitation are neglected. Excited electronic states of molecular and atomic species will be considered in future work, since they are known to accelerate the transition from non-equilibrium to thermal spark [10], [11].

## B. Conservation equations

Based on the previous assumptions, this section derives balance equations for mass, volume, mass density, total energy, electron energy, and gas energy [12]. The incident and reflected voltages along the transmission line connecting the power supply to the electrodes are computed by analytically solving the telegraph equations [13], with a time-varying resistance for the plasma as a boundary condition.

### 1. Total mass conservation

The total mass $m$ is conserved, since we are studying a closed system:

$$ \frac{dm}{dt} = 0 \qquad (1) $$

With a constant mass $m$, the evolution of the mass density $\rho = m/V$ is:

$$ \frac{d\rho}{dt} = -\frac{\rho}{V} \frac{dV}{dt} \qquad (2) $$

where $V$ is the discharge volume. It is assumed that pressure relaxes to atmospheric conditions with a timescale given by the rarefaction characteristic timescale $\tau_{rarefaction}$, such that:

$$ \frac{dp}{dt} = - \frac{p - p_{\text{atm}}}{\tau_{\text{rarefaction}}} \qquad (3) $$

where the rarefaction (or hydrodynamic) timescale is given by the ratio of the discharge radius over the sound velocity of the hot gases inside the plasma:

$$ \tau_{\text{rarefaction}} = \frac{r_d}{c_{\text{sound}}} = \frac{r_d}{\sqrt{\gamma \frac{R}{M} \bar{T}}} \qquad (4) $$

$\gamma$ is the adiabatic index, $R = 8.314 \text{ kJ} \cdot \text{mol}^{-1} \cdot \text{K}^{-1}$ is the molar gas constant, $\bar{M}$ is the mean molar mass of the mixture, and $\bar{T}$ is the (average) temperature in the gap. The volume is then determined by assuming a polytropic process:

$$ PV^k = cst \Leftrightarrow \frac{dV}{dt} = - \frac{V}{k} \frac{1}{p} \frac{dp}{dt} \qquad (5) $$

where $k$ is the polytropic index, assumed constant. If $k \to +\infty$, we retrieve an isochoric process. In the other cases, the volume will vary, and so will the mass density $\rho$. The special case $k = \gamma$ corresponds to an isentropic process.
## 2. Species mass fraction balance equation

The mass of species k is noted $m_k$. The rate of change of $m_k$ is only due to chemical reactions inside the system (no mass in nor out, since the system is closed), hence:

$$ \frac{dm_k}{dt} = \dot{m}_k^{gen} \qquad (6) $$

The generated mass can be written as $\dot{m}_k^{gen} = V\dot{\omega}_k M_k$, with $\dot{\omega}_k$ the molar rate of production of species k and $M_k$ the molar mass of species k. For a set $\mathcal{R}$ elementary reversible reactions, involving K chemical species, represented by the following equations:

$$ \forall r \in \mathcal{R}, \sum_{k=1}^{K} \nu_{kr} \chi_k \Leftrightarrow \sum_{k=1}^{K} \nu_{kr}^{\prime\prime} \chi_k \qquad (7) $$

where $\chi_k$ is the chemical symbol of the k-th species, and $\nu_{kr}^{\prime\prime}(\nu_{kr}^{\prime\prime})$ the stoichiometric coefficient of the k-th reactant (product) species in the r-th reaction, the molar rate of production is expressed as:

$$ \dot{\omega}_k = \sum_{r \in \mathcal{R}} (\nu_{kr}^{\prime\prime} - \nu_{kr}^{\prime}) q_r \qquad (8) $$

where the rate of progress variable of the r-th reaction is given by:

$$ q_r = k_{fr} \prod_{k=1}^{K} \left( \frac{n_k^{mol}}{V} \right)^{\nu_{kr}'} - k_{rr} \prod_{k=1}^{K} \left( \frac{n_k^{mol}}{V} \right)^{\nu_{kr}''} \qquad (9) $$

where $k_{fr}$ and $k_{rr}$ are the forward and reverse reaction rate constants of reaction r, and $n_k^{mol}$ is the molar quantity of species k. The mass of species k can be expressed using the mass fraction $Y_k: m_k = mY_k = \rho VY_k$. Since the mass is constant:

$$ \frac{dY_k}{dt} = \frac{\dot{\omega}_k M_k}{\rho} \qquad (10) $$

## 3. Total energy conservation

Applying the first principle of thermodynamics to the closed system gives:

$$ dU + dE_k + dE_p = \delta Q + \delta W \qquad (11) $$

where the change of the sum of internal energy U, kinetic energy $E_k$, and potential energy $E_p$ is related to the change of work W and heat Q. The variation of kinetic and potential energy is assumed to be negligible compared to the internal energy. For the time interval studied, thermal losses are neglected, thus $\delta Q = 0$. The electrical energy applied via the discharge is contained within $\delta W$, which also accounts for any mechanical work due to the change in volume. Taking the derivative with respect to time of Eq. (11) gives:

$$ \frac{dU}{dt} = \dot{W}_{\text{Joule}} + \dot{W}_{\text{p}} = \mathcal{P}_{\text{Joule}}V - p \frac{dV}{dt} \qquad (12) $$

where $\dot{W}_p = -p \frac{dV}{dt}$ is the power related to pressure forces. $\dot{W}_{\text{Joule}}$ is the injected power, related to the power density brought by the Joule effect $\mathcal{P}_{\text{Joule}}$, which can be expressed using the local form of Ohm's law:

$$ \mathcal{P}_{\text{Joule}} = j_{\text{p}}(t)E_{\text{p}}(t) = \sigma_{\text{p}}(t)E_{\text{p}}^{2}(t) \qquad (13) $$

where $\sigma_p$ is the electrical conductivity of the plasma (computed in Appendix A) and $E_p$ the electric field across the discharge. With the assumption of a homogeneous discharge, the electric field is written as the ratio of the plasma voltage $V_p(t)$ over the interelectrode gap $L_{gap}$, yielding:

$$ \mathcal{P}_{\text{Joule}}(t) = \sigma_{\text{p}}(t) \left( \frac{V_{\text{p}}(t)}{L_{\text{gap}}} \right)^2 \qquad (14) $$

In separate work to be presented in an upcoming paper, we demonstrated how to consider the reflections of the voltage wave sent by the generator from the plasma discharge. Under the assumptions of a time-varying resistance for the plasma and a lossless transmission cable, an analytical formula for the plasma voltage $V_p(t)$ can be derived:
$$ 
\begin{aligned}
V_p(t) &= (1 + \Gamma_p(t)) \sum_{n=0}^{N(t-L/c)} \alpha_g V_g \left( t - \frac{(2n + 1)L}{c} \right) \prod_{k=1}^{n} \left[ \Gamma_g \Gamma_p \left( t - \frac{2kL}{c} \right) \right] \\
&= \frac{2R_p(t)}{R_p(t) + Z_c} \left[ \alpha_g V_g \left( t - \frac{L}{c} \right) + \alpha_g V_g \left( t - \frac{3L}{c} \right) \Gamma_g \Gamma_p \left( t - \frac{2L}{c} \right) + \dots \right]
\end{aligned} 
\quad (15) $$

where $V_g$ is the generator voltage, $Z_c$ the cable impedance, $L$ the cable length, $c$ the wave propagation velocity in the cable, and $\alpha_g = Z_c/(Z_c + R_g)$ an attenuation coefficient (with $R_g$ the internal resistance of the generator). The reflection coefficient (for voltage) at the generator is given by:

$$ \Gamma_g = \frac{R_g - Z_c}{R_g + Z_c} \quad (16) $$

and the reflection coefficient (for voltage) at the plasma is:

$$ \Gamma_p(t) = \frac{R_p(t) - Z_c}{R_p(t) + Z_c} \quad (17) $$

Under the homogeneous discharge hypothesis, the plasma resistance $R_p(t)$ is calculated via:

$$ R_p(t) = \frac{L_{\text{gap}}}{S_d(t) \sigma_p(t)} \quad (18) $$

where $S_d(t) = \pi r_d^2(t)$ is the cross-sectional area of the plasma cylinder, with $r_d(t)$ the plasma radius. Each new incident wave is assumed to be applied across a homogeneous time-varying volume. Furthermore, the total internal energy $U$ in Eq. (12) is the sum of the mass internal energy of each species $u_k$ evaluated at their corresponding temperature $T_k$, weighted by their mass fraction $Y_k$:

$$ U = m u = m \sum_k u_k(T_k) Y_k \quad (19) $$

Note that the previous expression requires a temperature to be defined for each species, which is achieved here by assuming that species follow a Maxwellian distribution. The time derivative of total internal energy is then equal to:

$$ \frac{dU}{dt} = m \sum_k \left( u_k(T_k) \frac{dY_k}{dt} + Y_k \frac{du_k(T_k)}{dt} \right) + \frac{dm}{dt} \cdot \sum_k u_k(T_k) Y_k \quad (20) $$

Using conservation of mass and species (Eq. (1) and Eq. (10)), one then gets:

$$ \frac{dU}{dt} = m \sum_k \left( u_k(T_k) \frac{\dot{\omega}_k M_k}{\rho} + Y_k \frac{du_k(T_k)}{dt} \right) \quad (21) $$

Using the definition of the specific heat capacity at constant volume, $du_k(T_k) = c_{v,k}(T_k)dT_k$, Eq. (21) can be rearranged:

$$ \frac{dU}{dt} = V \sum_k \left( u_k(T_k) \dot{\omega}_k M_k + \rho Y_k c_{v,k}(T_k) \frac{dT_k}{dt} \right) \quad (22) $$

Injecting Eq. (22) into Eq. (12) and dividing by the volume, we obtain:

$$ \rho \sum_k \left( Y_k c_{v,k}(T_k) \frac{dT_k}{dt} \right) = - \sum_k (u_k(T_k) \dot{\omega}_k M_k) + P_{\text{Joule}} - \frac{p}{V} \frac{dV}{dt} \quad (23) $$

Finally, for a two-temperature system, with the heavy (neutral and ions) temperature equal to $T_g$ and the electron temperature to $T_e$, Eq. (23) yields:

$$ \rho \left\{ \sum_{k \neq e} \left( Y_k c_{v,k}(T_g) \frac{dT_g}{dt} + Y_e c_{v,e} \frac{dT_e}{dt} \right) \right\} = - \sum_{k \neq e} \left( u_k(T_g) \dot{\omega}_k M_k - c_{v,e} T_e \dot{\omega}_e M_e + P_{\text{Joule}} - \frac{p}{V} \frac{dV}{dt} \right) \quad (24) $$

The specific heat capacity at constant volume of an electron is given by:

$$ c_{v,e} = \frac{3 R}{2 M_e} \quad (25) $$

The thermodynamic properties of other gas species are detailed in Section C and depend on the gas temperature.

### 4. Electron energy balance equation

The electron energy distribution function $f_e$ is obtained by solving the Boltzmann equation. When $f_e$ is assumed to be Maxwellian, as given by Eq. IV-5.2 of Ref. [14], the electron energy conservation is given by:
$$ \frac{D}{Dt} \left( n_e \frac{3}{2} k_B T_e \right) = - \left( n_e \frac{5}{2} k_B T_e \right) \vec{\nabla} \cdot \vec{u} - \vec{\nabla} \cdot \vec{q}_e + \mathcal{P}_{\text{Joule}} - \dot{Q}_{\text{loss, vol}}^{\text{elas}} - \dot{Q}_{\text{loss, vol}}^{\text{inel}} \quad (26) $$

The left-hand side expresses the time rate of change of the electron thermal energy per unit volume, with $D/Dt = \partial/\partial t + \vec{u} \cdot \vec{\nabla}$ being the total derivative. The electron pressure is defined as:

$$ p_e = n_e k_B T_e \quad (27) $$

where $n_e$ is the number density of electrons, $k_B$ the Boltzmann constant and $T_e$ the electron temperature. The first term on the right-hand side of Eq. (26) represents the power required to change the electron volume. This can be rewritten using the electron pressure definition, the conservation of total mass (Eq. IV–4.1, [14]), and the relationship between volume and mass density (Eq. (2)):

$$ - \left( n_e \frac{5}{2} k_B T_e \right) \vec{\nabla} \cdot \vec{u} = - \left( \frac{5}{2} p_e \right) \left( -\frac{1}{\rho} \frac{D\rho}{Dt} \right) = - \left( \frac{5}{2} p_e \right) \left( \frac{1}{V} \frac{DV}{Dt} \right) = - \frac{5 p_e}{2V} \frac{DV}{Dt} \quad (28) $$

The second term on the right-hand side of Eq. (26) represents the electron heat flux, which is neglected here, given the short timescale studied. The third term on the right-hand side of Eq. (26) is the volumetric Joule power, as given by Eq. (13). The electric field accelerates the electrons, and not the ions, due to the mass difference between these two species. The fourth term on the right-hand side of Eq. (26) is the rate of electron energy loss per unit volume resulting from elastic collisions with heavy particles:

$$ \dot{Q}_{\text{loss, vol}}^{\text{elas}} = \frac{3}{2} n_e k_B (T_e - T_g) \sum_{k=1; k \neq e}^{K} \frac{2m_e}{m_k} \overline{\overline{v_{ek}}} \quad (29) $$

where $\overline{v_{ek}}$ denotes the average momentum transfer collision frequency between an electron and heavy particles of species k (see Appendix B). The last term on the right-hand side of Eq. (26) is the kinetic energy lost by electrons due to inelastic collisional processes:

$$ \dot{Q}_{\text{loss, vol}}^{\text{inel}} = \sum_{r \in \mathcal{R}_{ei}} \Delta H_r q_r \quad (30) $$

where the set of electron-impact (i.e., electron-temperature dependent) reactions is $\mathcal{R}_{ei}$, $q_r$ is the net rate of progress of the r-th reaction as detailed in Eq. (9), and $\Delta H_r$ is the net enthalpy change of the corresponding reaction. Due to the uniform discharge hypothesis, all terms in Eq. (26) involving a divergence operator are zero. Thus, the total derivative $D/Dt$ becomes a time-derivative $d/dt$ in Eq. (26) and in Eq. (28). Furthermore, using the equality $\frac{3}{2} k_B n_e = \rho Y_e c_{v,e}$, the mass species conservation (Eq. (10)), and the relationship between volume and mass density (Eq. (2)), one finally obtains the electron temperature balance equation:

$$ \rho Y_e c_{v,e} \frac{dT_e}{dt} = - \dot{\omega}_e M_e c_{v,e} T_e + \mathcal{P}_{\text{Joule}} - \dot{Q}_{\text{loss, vol}}^{\text{elas}} - \dot{Q}_{\text{loss, vol}}^{\text{inel}} - \frac{p_e}{V} \frac{dV}{dt} \quad (31) $$

### 5. Gas energy balance equation

The gas energy equation is obtained by subtracting the electron energy (Eq. (31)) from the total gas energy (Eq. (24)):

$$ \rho \sum_{k \neq e} \left( Y_k c_{v,k} (T_g) \right) \frac{dT_g}{dt} = - \sum_{k \neq e} \left( \dot{\omega}_k M_k u_k (T_g) \right) + \dot{Q}_{\text{loss, vol}}^{\text{elas}} + \dot{Q}_{\text{loss, vol}}^{\text{inel}} - \frac{p - p_e}{V} \frac{dV}{dt} \quad (32) $$

### 6. State equation

The state equation considered here is the multi-temperature perfect gas law, which relates the total pressure $p$ to the temperature $T_k$ and number density $n_k$ of each species:

$$ p = \sum_{k} n_{k} k_{B} T_{k} = \left( \sum_{k \neq e} n_{k} \right) k_{B} T_{g} + n_{e} k_{B} T_{e} = (n_{\text{tot}} - n_{e}) k_{B} T_{g} + n_{e} k_{B} T_{e} = n_{\text{tot}} k_{B} \bar{T} \quad (33) $$

where $\bar{T} = (1 - x_e)T_g + x_e T_e$ is defined as a mean temperature of the plasma ($x_e$ being the mole fraction of electrons).

### 7. Summary and numerical implementation

Fig. 2 summarizes the equations to solve and the required inputs. The chemical mechanism and thermodynamic data are detailed in Section C. The required electrical circuit inputs are the voltage waveform $V_g(t)$ and internal resistance $R_g$ of the generator, as well as the coaxial cable properties: length $L$, impedance $Z_c$, and wave propagation
velocity $c$. The interelectrode gap distance $L_{gap}$ is fixed. Initial values are required for the electron temperature $T_e$, gas temperature $T_g$, number density $n_k$, discharge radius $r_d$, and polytropic index $k$. The volume evolution (Eq. (5)), combined with pressure evolution Eq. (3)), species mass fraction balance (Eq. (10)), electron energy balance (Eq. (31)), and gas energy balance (Eq. (32)) are solved numerically. The chemical terms $\dot{\omega}_k$ are computed using Cantera 3.1 [15], around which an interface is built. The species mass balance and the gas and electron energy equations are solved using an implicit method based on backward-differentiation formulas (BDF) solver implemented in the Python package SciPy [16].

<figure><img src="image_3.png" /><figcaption>Fig. 2 Summary of the equations solved. The arrows represent the main dependencies between the blocks of the model.</figcaption></figure>

### C. Chemical mechanism from methane to acetylene

Experiments with NRP spark discharges at atmospheric pressure reveal that CH$_4$ is primarily converted to C$_2$H$_2$ and H$_2$ [3], [17], [18], [19], [20], [21]. Thus, we first focus on modeling the neutral and electron reactions, considering hydrocarbons up to C$_2$H$_6$ as well as their corresponding ions. A literature survey was conducted to obtain the reaction rates and thermodynamic data, taking into account their respective temperature ranges of validity. The rate constants of electron-impact reactions are computed assuming a Maxwellian distribution of electrons. Reverse reaction rates are calculated using detailed balance.

#### 1. Thermodynamic data and equilibrium

From a thermodynamic viewpoint only, at atmospheric pressure and in equilibrium, methane converts itself to hydrogen (H$_2$) and solid carbon (C(s)) between 1,000 K and 3,000 K. Above these temperatures, molecular dissociation of H$_2$ into atomic hydrogen (H) occurs, and the solid carbon vaporizes into gaseous carbon (C(g)). A further increase in temperature results in ionization, yielding hydrogen ions (H$^+$), carbon ions (C$^+$), and electrons (e$^-$). This process is illustrated in Fig. 2, where thermodynamic data obtained from the Glenn Research Center [27] are compared with the results of Fincke et al. [22]. The thermodynamic data used in Ref. [22] are not explicitly stated, but the main difference arises from the solid carbon properties employed.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Equilibrium composition of gaseous methane at atmospheric pressure. The mole fractions computed from thermodynamic data [23] (full line) are compared to the values of Fincke et al. [22] (dotted line).</figcaption></figure>
From a kinetic perspective, the conversion of methane can be divided into two sequential parts: a rapid endothermic reaction from methane to acetylene, ethylene, and benzene, followed by a slower exothermic reaction associated with the growth of solid carbon [22], [24]. If the products are rapidly cooled, the composition is frozen - soot cannot nucleate and grow, and the equilibrium condition is never reached. Thus, if solid carbon is removed from the list of possible products, the composition plotted in Fig. 4 is obtained.

<figure><img src="image_5.png" /><figcaption>Fig. 4 (Metastable) equilibrium composition of gaseous methane at constant atmospheric pressure, without including solid carbon C(s). The mole fractions are computed from thermodynamic data [23] (full line) and compared to the values of Fincke et al. [22] (dotted line).</figcaption></figure>

In nanosecond thermal spark discharges, electron temperatures of up to 50,000 K have been reported [5]. In the case of thermalization, the gas temperature will also achieve such high temperatures. To correctly compute the energy conservation equations (Eq. (31)–(32)), thermodynamic data up to 50,000 K are required. Valid properties are also required to calculate reverse reactions. Table 2 reports the species considered, along with the sources of their thermodynamic properties. As justified later in Section C-2, only hydrocarbons up to C₂H₆ are considered. Note that negative ions are not present, and exotic species, such as H₃⁺ or CH₅⁺ [25], are ignored at this stage. It is also assumed that all species are in their fundamental states - no excited electronic or vibrational states are considered. Excited states of molecular and atomic species will be considered in future work, since they are known to accelerate the transition from non-equilibrium to thermal spark [10], [11].

<table><caption>Table 2. List of species considered, source of their thermodynamic properties, and temperature range of validity. Refer to the text for details on how the extension is made, as well as the computed data.</caption><thead><tr><th>Ref.</th><th>Temperature range</th><th>Extended to</th><th>Species</th></tr></thead><tbody><tr><td rowspan="3">GRC [23]</td><td>≥ 300 K</td><td>-</td><td>e<sup>-</sup>, H<sup>+</sup></td></tr><tr><td>300–20,000 K</td><td>50,000 K</td><td>CH, CH<sup>+</sup>, C, e<sup>+</sup>, C<sub>2</sub>, C<sub>2</sub><sup>+</sup>, H<sub>2</sub>, H<sub>2</sub><sup>+</sup>, H</td></tr><tr><td>300–6,000 K</td><td>50,000 K</td><td>CH<sub>4</sub>, CH<sub>3</sub>, CH<sub>2</sub>, C<sub>2</sub>H<sub>6</sub>, C<sub>2</sub>H<sub>5</sub>, C<sub>2</sub>H<sub>4</sub>, C<sub>2</sub>H<sub>3</sub>, C<sub>2</sub>H<sub>2</sub>, C<sub>2</sub>H</td></tr><tr><td>Burcat [26]</td><td>300–6,000 K</td><td>50,000 K</td><td>CH<sub>4</sub><sup>+</sup>, CH<sub>3</sub><sup>+</sup>, CH<sub>2</sub><sup>+</sup>, C<sub>2</sub>H<sub>6</sub><sup>+</sup>, C<sub>2</sub>H<sub>5</sub><sup>+</sup>, C<sub>2</sub>H<sub>4</sub><sup>+</sup>, C<sub>2</sub>H<sub>3</sub><sup>+</sup>, C<sub>2</sub>H<sub>2</sub><sup>+</sup>, C<sub>2</sub>H<sup>+</sup></td></tr><tr><td>Computed</td><td>300–50,000 K</td><td>-</td><td>C<sup>+</sup>, C<sup>++</sup>, C<sup>+++</sup></td></tr></tbody></table>

Species standard-state thermodynamic properties are reported in the NASA9 format [15], [23], where nine coefficients represent the temperature dependence of the heat capacity at constant pressure c<sub>p</sub>, the enthalpy h and the entropy s:

$$ 
\begin{aligned}
\frac{\hat{c}_p^\circ(T)}{R} &= a_0 T^{-2} + a_1 T^{-1} + a_2 + a_3 T + a_4 T^2 + a_5 T^3 + a_6 T^4 \\
\frac{\hat{h}^\circ(T)}{RT} &= -a_0 T^{-2} + a_1 \frac{\ln T}{T} + a_2 + \frac{a_3}{2} T + \frac{a_4}{3} T^2 + \frac{a_5}{4} T^3 + \frac{a_6}{5} T^4 + \frac{a_7}{T} \\
\frac{\hat{s}^\circ(T)}{R} &= -\frac{a_0}{2} T^{-2} - a_1 T^{-1} + a_2 \ln T + a_3 T + \frac{a_4}{2} T^2 + \frac{a_5}{3} T^3 + \frac{a_6}{4} T^4 + a_8
\end{aligned}
\qquad (34) $$

For atoms, diatomic molecules, and their corresponding ions, data are obtained from the GRC database [23] and are valid up to 20,000 K, except for electrons and H⁺, whose values are analytically determined for all relevant temperatures. For neutral species with more than two atoms, coefficients are also given in [23], and are valid up to 6,000 K. Their corresponding ion properties come from the Burcat database [26], also valid up to 6,000 K. To avoid
numerical error in the extrapolation performed by Cantera [15] outside these ranges of validity[^7], the thermodynamic values are extended up to 50,000 K using the following procedure. The heat capacity $c_p$ at the highest temperature available is assumed to be constant up to 50,000 K. The enthalpy (entropy) is then computed by integrating this constant $c_p$ ($c_p/T$). Finally, coefficients of Eq. (34) are obtained via a fitting algorithm detailed in the Python package pMuTT [27]. We justify this approach with the following considerations. For polyatomic species, the translational, rotational, and vibrational modes are fully excited above a few thousand kelvins. Thus, at 6,000 K, most species have reached a plateau in their heat capacity. This approach neglects electronic excitation, but its impact on heat capacity is important only for $T \gtrsim 10,000$ K, at which point molecular species (with more than two atoms) are dissociated. For atomic and diatomic species (and their ions), electronic excitation is already included, as their properties are calculated up to 20,000 K. Above 20,000 K, these species will also dissociate and ionize before reaching this temperature. To summarize, our approach neglects the electronic excitation of species by assuming a constant $c_p$ above 6,000 K to achieve numerical convergence and avoid non-physical behaviour. The results of this procedure are shown in Fig. 4 for acetylene, with both the simple extrapolation and the corrected values displayed. In the simple approach, at 40,000 K, the divergence of the free Gibbs energy leads to the presence of $C_2H_2$ in equilibrium, which is non-physical.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Thermodynamic properties of acetylene against temperature, in the validity range of [23] (green full line), and numerical divergence arising when evaluating at high temperature (red dashed line), as computed with Cantera 3.1 [15]. The proposed extrapolation method is also illustrated (black full line). Subplot (a) represents the molar heat capacity at constant pressure, (b) the molar enthalpy, (c) the molar entropy, and (d) the molar Gibbs free energy.</figcaption></figure>

The above extrapolation method does not work for the carbon ion ($C^+$), which ionizes into $C^{++}$ at temperatures above 20,000 K. Therefore, for completeness, the thermodynamic properties of $C^+$ are recomputed from the partition function, using the degeneracy and associated energy given by the NIST Atomic Spectra Database [28]. Then, coefficients of Eq. (34) are again fitted using pMuTT [27], and compared with the available data of $C^+$ in its validity range of the GRC database [23], as illustrated in Fig. 6. This methodology is also applied to compute the properties of $C^{++}$ and $C^{+++}$.

<footnote>[^7]: Even though temperature ranges can be provided in Cantera 3.1 [15] for thermodynamic properties, there are not hard limits and no error/warning is displayed if the temperature goes outside of this interval.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6 Computation of thermodynamic properties of C⁺ using the partition function with energies and degeneracies from NIST Atomic Spectra Database [28] (black full line), compared with data in their validity range available in [23] (green dashed line). Subplot (a) represents the molar heat capacity at constant pressure, (b) the molar enthalpy, (c) the molar entropy, and (d) the molar Gibbs free energy.</figcaption></figure>

Once we have all the properties for the given species, an equilibrium computation is done by minimizing the Gibbs free energy in Cantera. The equilibrium composition at atmospheric pressure is shown in Fig. 7, starting from an initial composition of pure CH₄. Methane is decomposed into hydrogen and acetylene (when there is no solid carbon) up to 3,000 K, after which dissociation occurs – around 10,000 K, ionization starts, yielding electrons, ionized hydrogen, and ionized carbon. C⁺⁺ appears above 20,000 K, while C⁺⁺⁺ becomes important around 35,000 K. For verification, the total enthalpy is compared to the computation of Wu *et al.* [29], as shown in Fig. 8.

<figure><img src="image_8.png" /><figcaption>Fig. 7 (Metastable) equilibrium composition of methane at constant atmospheric pressure (solid carbon C(s) is not included). The label “others” is the sum of the other species listed in Table 2.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8 Enthalpy of a mixture of methane in equilibrium state, at atmospheric pressure and temperatures up to 50,000 K. Black line corresponds to calculated values, and red crosses are the computations of Wu et al. [29].</figcaption></figure>

## 2. Kinetic mechanism

As stated above, the conversion of methane can be divided into two sequential parts: a rapid endothermic conversion from methane to acetylene, followed by slower exothermic reactions describing the growth of solid carbon [22], [24]. In NRP discharges, experiments have shown that acetylene is one of the primary products at atmospheric pressure [3], [19], [21], [30]. Other major products include hydrogen, ethylene, and ethane. Minor proportions of higher hydrocarbons are found in gas chromatography, and traces of solid carbon are found in these experiments. Based on these observations, a chemical mechanism is constructed, which includes the following species: $H_y$, $CH_y$, and $C_2H_y$. All corresponding ions are also added, as well as electrons, as reported in Table 2.

A literature survey of some chemical mechanisms proposed for methane plasmalysis was conducted. We only considered the reactions involving methane, neglecting all reactions with oxygen and nitrogen used to model combustion. The reaction rates and cross-sections proposed by Heijkers et al. [25], Yuan et al. [31], and Fincke et al. [22] were reviewed. The final kinetic mechanism is not yet reported, as ongoing work is being carried out to determine the relevant reaction rates and assess their sensitivity. However, the methodology used is presented here. All forward reaction rates $k_f$ are reported in the extended Arrhenius format given by:

$$k_f(T) = AT^b \exp\left(-\frac{E_a}{RT}\right), T = T_g \text{ or } T_e \qquad (35)$$

where A is the pre-exponential constant, b the power exponent, and $E_a$ the activation energy. $k_f$ depends on the electron temperature whenever the electron is present as a reactant, whereas the gas temperature is used for all other reaction rates. Reverse reaction rates are computed using the detailed balance principle:

$$k_r(T) = \frac{k_f(T)}{K(T)} \qquad (36)$$

where K is the equilibrium constant computed through the thermodynamic properties defined in Eq. (34). While this computation is automatically done in Cantera 3.1 [15] for reaction rates depending on $T_g$, it is not the case for electron reactions. Therefore, reverse rates for reactions involving electrons are computed using an in-house routine that first uses Eq. (36) to calculate $k_r(T)$ for a range of temperatures. $k_r(T)$ is then fit to an Arrhenius law defined by Eq. (35) using a least-square minimization.

Our mechanism includes the 46 neutral-neutral reactions considered by Fincke et al. [22], with rate constants updated from the NIST Kinetics database [32]. In cases where multiple rate constants were available for a reaction, preference was given to second-order reactions with the widest temperature ranges. For ion-neutral reactions, the UMIST 2022 Database [33] was used. For reactions with electrons, rate constants are computed assuming a Maxwellian distribution at $T_e$, yielding:

$$k_f(T_e) = \int_0^{+\infty} \sigma(E) v_{T_e}(E) f_{T_e}(E) dE \qquad (37)$$

In Eq. (37), $k_f$ is the forward reaction rate constant, $T_e$ the temperature of electrons, $\sigma$ the cross section of the electron-impact reaction, $v_{T_e}$ the thermal velocity of electrons, $f_{T_e}$ the electron energy distribution function (assumed Maxwellian), and the integral is carried out over all energies E. Various cross-section databases were explored, including LXCAT [34], [35], [36] and Quantemol [37]. For consistency, all electron impact cross-sections of
hydrocarbons C$_{x}$≥1H$_{y}$≥1 are taken from the Janev databases [38], [39]. Those for H$_{y}$ species are taken from the Biagi [40], ISTLisbon [41], [42], and Morgan [43] databases. Those of C$_{y}$ species come from BSR [44], [45], [46], Tawara [47], and Quantemol [48]. For recombination reactions of C$^{+}$ and H$^{+}$ with electrons, the UMIST 2022 Database [33] was used, and the rate constant of the H$_{2}^{+}$ recombination reaction was taken from the KIDA Database [49]. In total, more than 300 reactions were included in this mechanism, along with specified temperature ranges and references, summarized in Table 3.

<table><caption>**Table 3. List of rate constants used, and their references. Reverse rate constants are automatically calculated by Cantera 3.1 [15] for heavy-heavy reactions (*). Arrhenius fits (Eq. (35)) of the rate constants determined by detailed balance (Eq. (36)) are used for electronic reactions (**).**</caption><thead><tr><th>Ref.</th><th>Reaction type</th><th># of forward reactions</th><th># of reverse reactions</th></tr></thead><tbody><tr><td>NIST Kinetics database [32]</td><td>Neutral-Neutral</td><td>46</td><td>46 (*)</td></tr><tr><td>UMIST 2022 Database [33]</td><td>Neutral-Ion</td><td>103</td><td>103 (*)</td></tr><tr><td>Janev databases [38], [39]<br>(C<sub>x≥1</sub>H<sub>y≥1</sub>) + others (see text)</td><td>Electron-Heavy</td><td>228</td><td>209 (**)</td></tr></tbody></table>

The correct implementation of reverse reactions is verified by comparing the thermodynamic equilibrium composition with the steady-state chemical composition obtained from simulations at constant atmospheric pressure and constant temperature ($T_e = T_g = T$). We will compare the chemical evolution of the species, including or neglecting reverse reactions, starting from a mixture of pure CH$_4$. Fig. 9–11 show the evolution of chemical species at $T = 5,000, 15,000$ and $30,000$ K. At $5,000$ K, the plasma reaches the equilibrium composition whether reverse electron reactions are included or not, since there are no electrons and ions present at this temperature. At $15,000$ K, however, the equilibrium composition can only be reached if reverse electron reactions are included (Fig. 10). At $30,000$ K, numerical simulations diverge if reverse reactions are not considered, as illustrated in Fig. 11. However, the simulations converge towards equilibrium if reverse electron reactions are included. This verification was performed up to $50,000$ K. We also observe in these three figures that methane converts into CH$_3$ and H$_2$ for the “low” temperature cases, and into C and H at higher temperatures. For the “high” temperature cases, ionization starts around $1–10$ µs at $15,000$ K, and around $10–100$ ns at $30,000$ K. At $50,000$ K (not shown here), quasi-full ionization occurs in around $10$ ns.

<figure><img src="image_10.png" /><figcaption>**Fig. 9 Evolution of chemical species at constant atmospheric pressure and temperature, with $T_e = T_g = 5,000$ K. Thermodynamic equilibrium is represented by dashed lines. Reverse electron impact reactions have no influence on the results at this temperature.**</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10 Evolution of chemical species at constant atmospheric pressure and temperature, with $T_e = T_g = 15,000$ K. Thermodynamic equilibrium is represented by dashed lines. Without the reverse reaction (long dashes), the simulation diverges from equilibrium around 100 µs.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Evolution of chemical species at constant atmospheric pressure and temperature, with $T_e = T_g = 30,000$ K. Thermodynamic equilibrium is represented by dashed lines.</figcaption></figure>

### III. Results and discussion

Numerical simulations are performed to reproduce and compare with the experimental data of Delahaie et al. [6]. In that reference, NRP discharges are applied between two electrodes inside a reactor filled with methane at atmospheric pressure. The inter-electrode gap distance is $L_{gap} = 3.8$ mm, and the pulses are applied at a repetition frequency of $f = 10$ kHz. Methane is injected into the reactor with a volumetric flow rate of 2 splm. Nitrogen is introduced in trace amounts (less than 1% in volume) as a temperature tracer; however, the impact of $N_2$ on the kinetics will be neglected here, and we will assume that the initial plasma volume is composed of only $CH_4$ and a trace amount of seed electrons (resulting from previous pulses). The electrical setup is briefly summarized in Fig. 12, where a generator imposes its high voltage at the boundary ($x = 0$) of a cable. The incident wave travels along the transmission line to the load at the other end ($x = L = 6.2$ m), where it is absorbed and/or reflected depending on the impedance-matching conditions. The cable is characterized by an impedance $Z_c = 75$ Ω and a wave propagation velocity $c = 1.9 \ 10^8$ m·s⁻¹. Voltage and current probes are placed in the middle of the cable at $x_{meas} = 3.1$ m. Average electrical measurements are presented in Fig. 13, where we observe a first incident voltage pulse between 0 and ~25 ns (labelled as '$1 \to $'), corresponding to the generator voltage pulse<fnref n="8" /> traveling toward the load. At the plasma end, the voltage pulse is partially absorbed but mainly reflected – this can be seen from the fact that the reflected voltage $V^{1\to}$ has a smaller amplitude relative to the incident $V^{1\to}$. A second set of incident/reflected waves occurs, with a stronger attenuation of $V^{2\to}$ (second reflected pulse) compared to $V^{2\to}$ (second incident pulse), meaning that the plasma resistance is dropping, but still above the cable impedance. Finally, in the third set of reflections, the voltage passes through zero for $V^{3\to}$, indicating that the plasma resistance falls below the cable impedance. These successive reflections continue until all the energy is dissipated. Note that these are real reflections of electrical waves, and not another pulse being emitted by the generator (since $\tau_{PRF} = 1/f = 100$ µs ≫ 100 ns). Delahaie et al. [6] also

<footnote n="8">⁸ Strictly speaking, the propagating voltage pulse has an amplitude $\alpha_g V_g(t)$, with $\alpha_g = \frac{Z_c}{R_g+Z_c}$ the attenuation factor, close to 1 when $R_g \ll Z_c$, as is the case here.</footnote>
performed optical measurements of the discharge, obtaining time-resolved electron number density from optical emission spectroscopy of $H_\alpha$, as well as the evolution of the discharge diameter.

<figure><img src="image_13.png" /><figcaption>Fig. 12 Schematic representation of the electrical setup. The high-voltage pulse sent by the generator travels along the coaxial cable and is absorbed and reflected at the plasma.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13 Voltage and current measurements at the middle of the cable (averaged over 1000 pulses). The 2-standard deviation is indicated by the shaded area. The time origin is taken when signals appear at probes.</figcaption></figure>

For the simulations, numerous inputs are required. These are specified in Fig. 1 above. The simulation is run for 200 ns, at which point thermal losses (by radiation) cease to be negligible. Four cases will be considered, where most of the initial conditions are summarized in Table 4. Case (A) corresponds to an isochoric simulation, while Cases (B), (C), and (D) undergo an isentropic expansion. The initial radius of the discharge is 180 µm, except in Case (C), where the radius is 100 µm. The generator voltage time trace is simplified to a trapezoid shape, as shown in Fig. 14. For Case (D), the maximal amplitude of the trapezoid is increased to 10 kV. All simulations start in pure CH₄, with an initial number density of 10¹³ cm⁻³ for electrons and CH₄⁺ ions. The internal resistance of the generator is set to 1 Ω, so that $\alpha_g \approx 1$, hence $V^{1\rightarrow} = V_g$. On the numerical side, the BDF solver is setup with a maximum order of 5, a first timestep of 10⁻¹⁵ s, a maximum timestep of 10⁻¹⁰ s, an absolute tolerance of 10⁻¹⁶ and a relative tolerance of 10⁻¹⁰. Each case runs in less than ten minutes on a standard laptop computer.

<table><caption>Table 4. Parameters used for each case. Case A considers an isochoric volume, Case B involves an expansion process, Case C reduces the initial radius, and Case D utilizes a higher generator voltage. The generator voltage profiles noted (*) and (**) are presented in Fig. 14.</caption><thead><tr><th rowspan="2">Case</th><th rowspan="2">k</th><th colspan="3">Cable</th><th colspan="2">Generator</th><th colspan="2">Geometry</th><th colspan="4">Initial condition</th></tr><tr><th>Z<sub>c</sub> [Ω]</th><th>c [m.s<sup>-1</sup>]</th><th>L [m]</th><th>V<sub>g</sub> (t)</th><th>R<sub>g</sub> [Ω]</th><th>r<sub>d</sub> [µm]</th><th>L<sub>gap</sub> [mm]</th><th>T<sub>g</sub> [K]</th><th>T<sub>e</sub> [K]</th><th>P [atm]</th><th>n<sub>e</sub> [cm<sup>-3</sup>]</th></tr></thead><tbody><tr><td>(A)</td><td>+∞</td><td>75</td><td>1.9 10<sup>8</sup></td><td>6.2</td><td>(*)</td><td>1</td><td>180</td><td>3.8</td><td>1000</td><td>1000</td><td>1</td><td>10<sup>13</sup></td></tr><tr><td>(B)</td><td><b>5/3</b></td><td>75</td><td>1.9 10<sup>8</sup></td><td>6.2</td><td>(*)</td><td>1</td><td>180</td><td>3.8</td><td>1000</td><td>1000</td><td>1</td><td>10<sup>13</sup></td></tr><tr><td>(C)</td><td>5/3</td><td>75</td><td>1.9 10<sup>8</sup></td><td>6.2</td><td>(*)</td><td>1</td><td><b>100</b></td><td>3.8</td><td>1000</td><td>1000</td><td>1</td><td>10<sup>13</sup></td></tr><tr><td>(D)</td><td>5/3</td><td>75</td><td>1.9 10<sup>8</sup></td><td>6.2</td><td>(**)</td><td>1</td><td>180</td><td>3.8</td><td>1000</td><td>1000</td><td>1</td><td>10<sup>13</sup></td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Fig. 14 Measured incident voltage (taken equal to generator voltage) of Ref. [6] (black line), and trapezoidal signals (*) and (**) used as simulation input (green and magenta lines).**</figcaption></figure>

The simulation gives the evolution in time of the density of each species, the electron and gas temperatures, the pressure, the resistance, the radius, the current through the discharge, and the voltage across the plasma. First, we present results for Case (A). Next, comparisons of electron density are presented with experimental data obtained by Delahaie *et al.* [6] for all cases considered.

At the beginning of the discharge, the gas is mainly composed of neutral methane, and the electrical resistance is high ($R_p \gg Z_c$), as can be seen in Fig. 15. Therefore, according to Eq. (15), the plasma voltage is equal to twice the incident voltage ($V_p(t) \approx 2V_g(t - L/c)$), as shown in Fig. 16. Simultaneously, the resulting electric field accelerates the electrons, whose temperature rise is proportional to the plasma voltage. At the same time, ionization occurs, and the electron number density increases. The conductivity starts to increase as the discharge forms, and the resistance falls between 1 ns and 10 ns in Fig. 15. As the resistance becomes equal to or less than the cable impedance, the plasma voltage is modified accordingly to $V_p(t) \approx 2\frac{R_p(t)}{Z_c}V_g(t - L/c)$, as seen in Fig. 16 after 2 ns. The electron temperature $T_e$ is still proportional to the plasma voltage during this period. Also, as the number of electrons increases, the exchange of energy by collision between electrons and heavy particles increases, resulting in a nanosecond rise of the gas temperature up to 5,000 K. However, thermalization does not occur before the end of the pulse. At $t = 12$ ns, the applied voltage decreases, and the electron temperature also diminishes. For $t > 12$ ns, the electron temperature (Fig. 16) decreases because there is no longer any injected power, and electrons start to recombine (shown later in Fig. 17). The gas temperature remains constant. Then, around 65 ns (corresponding to $2L/c$), a reflected pulse reaches the electrode, reigniting the discharge. The plasma voltage starts inverted due to the reflection off of the generator of the previous pulse. The process occurs again, with a lesser impact on the gas temperature due to less energy being deposited.

<figure><img src="image_16.png" /><figcaption>**Fig. 15 Time evolution of the plasma resistance for Case (A).**</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16 Electron temperature $T_e$ (black), gas temperature $T_g$ (red), generator voltage $V_g$ shifted to account for the propagation delay between the generator and the plasma (blue), and plasma voltage (green), for Case (A). The graph on the left shows a zoom on the first 20 ns on a logarithmic scale. The graph on the right shows the whole simulation in a linear scale.</figcaption></figure>

We compare in Fig. 17 the variations of the electron number density with experimental data. The pressure evolution is shown in Fig. 18. For Case (A), the radius is constant, while the pressure increases. The electron number density of the first pulse is overestimated by a factor of 5 compared to the measurements. The second and third reflections increase the density of electrons, but this time underestimating it by an order of magnitude. The measurements suggest that the discharge transitions from a non-equilibrium spark (first discharge) to a thermal spark (second and third discharges), which is not captured in Case (A).

Let us now turn our attention to the impact of pressure relaxation – Case (A) vs (B) in Table 4 – and its influence on the electron number density shown in Fig. 17. Case (B) considers that the pressure relaxes towards atmospheric conditions, see Fig. 18. We could expect that the decrease of total density induced by the expansion will increase the reduced electric field significantly. However, as shown in Fig. 17, Case (B) is similar to Case (A) in terms of $n_e$, suggesting that the expansion has a moderate impact on the overall kinetics of the non-equilibrium spark.

Other initial conditions were considered, such as changing the initial radius to 100 µm instead of 180 µm for Case (C). The pressure increases to higher values, as energy is deposited in a smaller volume and decompresses more quickly due to the smaller radius. This time, the electron density reaches the high values measured (≈ $10^{18}$ cm⁻³). Nevertheless, the trends do not follow the experimental points. In Case (D), a higher generator voltage is applied compared to Case (B). The physical outputs are consistent: $n_e$, *p* and *r* exhibit higher values and occur more quickly. The decays of $n_e$ after the 2nd and 3rd pulses correspond qualitatively to the trends observed experimentally.

The radius evolution is compared to optical measurements in Fig. 19. Case (A) matches reasonably well the measurements for $t > 75$ ns, corresponding to the second pulse, but misses the larger initial radius between 0 and 10 ns. The isentropic expansion of Case (B) does not correspond to the observed trend, which suggests that the pressure relaxation is too quick.

These simulations demonstrate that reasonable agreement with the experimental data can be achieved with Case (A). Adjusting the initial radius or the applied voltage in the simulation significantly improves the agreement in electron number density and needs to be reviewed thoroughly. Other hypotheses could also be refined. Primarily, the inclusion of excited (electronic) species is known to accelerate the transition from non-equilibrium to thermal spark [10], [11]. The assumption of a Maxwellian distribution should also be relaxed.
<figure><img src="image_18.png" /><figcaption>Fig. 17 Comparison of electron density evolution, for the various cases considered, with experimental data from Delahaie *et al.* [6] represented as black dots.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18 Comparison of pressure evolution for the various cases considered.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19 Comparison of radius evolution, for the various cases considered, with experimental data from Delahaie *et al.* [6] represented as dots.</figcaption></figure>
IV. Conclusion

In this paper, coupled kinetic and electrical simulations were conducted to analyse methane plasmalysis towards
hydrogen and acetylene during the first hundreds of nanoseconds of a spark discharge. A kinetic mechanism for the
conversion of methane to acetylene was constructed, including all H<sub>y</sub>, CH<sub>y</sub>, and C<sub>2</sub>H<sub>y</sub> species in their ground states,
as well as their corresponding ions and electrons. Excited vibrational and electronic states were neglected. Special
care was taken to consider the thermodynamic properties of all these species in order to calculate the reverse reaction
rates and the enthalpy of each species. The balance equations of a two-temperature homogeneous system are solved
under the assumption of a Maxwellian distribution at T<sub>e</sub> for electrons and at T<sub>g</sub> for heavy species. Pressure relaxation
is also considered. Particular attention was given to the injected power, which requires considering the generator
voltage as a wave traveling towards the plasma and reflecting off of it. This results in multiple round-trips between
the generator and electrodes over the time studied.

For now, simulations can reproduce the electron density observed experimentally for the first pulse in the non-equilibrium regime with an accuracy better than an order of magnitude. Secondary and tertiary pulses do produce electrons in the model, but are off by more than an order of magnitude from the experiment, which reaches the high electron densities (10<sup>18</sup> cm<sup>-3</sup>) typical of a thermal spark. Alternative cases were run, demonstrating that higher electron number densities can be achieved with a higher applied voltage or a smaller initial radius.

In future work, excited electronic states of atomic and molecular species will be included, since they are known to
accelerate the transition from non-equilibrium to thermal spark and the electron production. The code will also be
compared with other solvers to assess its validity. For now, within an order of magnitude, the simulation provides an
estimate of the plasma resistance, which can be used to design and regulate a power source. Moreover, the hydrogen
and acetylene yield in the plasma can serve as a source term for a 3D simulation of the afterglow of the discharge.
# Appendix

## A. Computing the electrical conductivity $\sigma_p$

To compute the plasma resistance $R_p$, the electrical conductivity of the plasma $\sigma_p$ is required. According to Ref. [14] (chapter II, section 13), the conductivity has the following expression:

$$ \sigma_p = \frac{n_e e^2}{m_e \overline{v_{eH}}} \qquad (38) $$

where $n_e$ is the electron density, $e$ the electric charge, $m_e$ the electron mass, and $\overline{v_{eH}} = \sum_k \overline{v_{ek}}$ the average momentum transfer collision frequency of an electron with all heavy particles ($\overline{v_{ek}}$ are computed in Appendix B). $\overline{v_{eH}}$ can be decomposed as a sum of collision frequencies with neutrals $\overline{v_{eN}}$ and with ions $\overline{v_{el}}$. In the special case where $\overline{v_{el}} \gg \overline{v_{eN}}$, the gas may be seen as fully ionized (II–8, [14]), and one may use the Spitzer–Härm conductivities:

$$ \sigma_p[\Omega^{-1} \cdot m^{-1}] = 1.98 \frac{n_e e^2}{m_e \overline{v_{el}}} = 1.53 \ 10^{-2} \frac{(T_e[\text{K}])^{3/2}}{\ln \Lambda} \qquad (39) $$

In the code, the threshold at which the plasma is considered fully ionized is when $\overline{v_{el}} > \overline{v_{eN}}$. To check the code implementation, a computation at equilibrium is performed, and compared with the results of Niu et al. [50]:

<figure><img src="image_21.png" /><figcaption>Fig. 20 Electrical conductivity of an equilibrium composition of methane against temperature (dark line), as computed using Eq. (38). A factor 1.98 is introduced (Eq. (39), compared to Eq. (38)) when $\overline{v_{el}} > \overline{v_{eN}}$, hence the sharp increase around 8,000 K. Numerical results computed in local thermodynamic equilibrium by Niu et al. [50] are also plotted (red dot) for comparison.</figcaption></figure>

## B. Computing momentum-transfer collision frequencies $\overline{v_{ek}}$

To compute the rate $\dot{Q}_{loss, vol}^{elas}$ at which electrons lose energy per unit volume as a result of elastic collisions with heavy particles, it is necessary to know $\overline{v_{ek}}$, the energy-weighted average momentum transfer collision frequency between an electron and heavy particles of species $k$, as indicated in Eq. (29). This frequency is determined following the method proposed in section 6 of chapter II of Ref. [14], which is based on an approximation of the Chapman-Enskog theory and whose results are recalled thereafter. The average momentum-transfer collision frequency (between an electron and a heavy species $k$) is given by:

$$ \overline{v_{ek}} = n_k v_{th,e} \overline{Q_{ek}} \qquad (40) $$

where $n_k$ is the number density of species $k$, $v_{th,e} = \sqrt{\frac{8k_B T_e}{\pi m_e}}$ the thermal velocity of the electrons, and $\overline{Q_{ek}}$ the average momentum-transfer cross section defined as:

$$ \overline{Q_{ek}} = \frac{2}{3} \int_{0}^{+\infty} x^2 e^{-x} \overline{Q_{ek}}(k_B T_e x) dx \qquad (41) $$

where $\overline{Q_{ek}}$ is the energy-dependent momentum-transfer cross section. In the case where $\overline{Q_{ek}}$ is constant, Eq. (41) reduces to $\overline{Q_{ek}} = \frac{4}{3} \overline{Q_{ek}}$. Data for the momentum-transfer cross section comes from multiple sources. For stable neutral species, measured values of $\overline{Q_{ek}}$ are reported in database like LXCAT<fnref n="9" /> [34], [35], [36]. For other neutral species, the

<footnote n="9">⁹ On LXCat, the “elastic” cross-section was chosen, not the “effective”.</footnote>
approximation of a rigid sphere with a constant radius $r_k$ yields a constant average cross-section $\overline{Q}_{ek} = \frac{4}{3}\pi r_k^2$. For ions, it can be accurately calculated using the Coulomb potential and the resulting Rutherford cross-section. The computation involves the Coulomb logarithm $\ln \Lambda$, which considers the impact of charged particles to shield themselves from the potential. In the restrictive case of $\ln \Lambda \gg 1$, the average cross-section for a collision between an electron and an ion is given by [14]:

$$ \overline{Q}_{ei} = 6\pi\overline{b}_0^2 \ln \Lambda \qquad (42) $$

where $\Lambda$ is calculated as the ratio of the Debye length $\lambda_D = \sqrt{\frac{\epsilon_0 k_B T_e}{n_e e^2}}$ and the impact parameter $\overline{b}_0 = \frac{ze^2}{12\pi\epsilon_0 k_B T_e}$. Table 5 summarizes the choice of momentum-transfer cross-sections.

<table><caption>Table 5. List of references used for the energy-dependent momentum-transfer cross-section $\tilde{Q}_{ek}$, and the equation used to compute the average momentum-transfer cross section $\overline{Q}_{ek}$, for each species considered.</caption><thead><tr><th>Species</th><th>Q̃<sub>ek</sub></th><th>Q̃<sub>ek</sub></th></tr></thead><tbody><tr><td>CH<sub>4</sub>, C<sub>2</sub>H<sub>6</sub>, C<sub>2</sub>H<sub>4</sub>, C<sub>2</sub>H<sub>2</sub></td><td>Hayashi [51], [52], [53]</td><td>Eq. (41)</td></tr><tr><td>H</td><td>Morgan [43]</td><td>Eq. (41)</td></tr><tr><td>H<sub>2</sub></td><td>Biagi [40]</td><td>Eq. (41)</td></tr><tr><td>C</td><td>BSR [44], [45], [46]</td><td>Eq. (41)</td></tr><tr><td>Neutrals without elastic momentum-transfer cross section on LXCat [34], [35], [36]</td><td>Q̃<sub>ek</sub> = πr<sub>0</sub><sup>2</sup>, with r<sub>0</sub> = 1 Å</td><td>Q̃<sub>ek</sub> = <sup>4</sup>&frasl;<sub>3</sub> Q̃<sub>ek</sub></td></tr><tr><td>Ions (CH<sub>4</sub><sup>+</sup>, CH<sub>3</sub><sup>+</sup>, CH<sub>2</sub><sup>+</sup>, CH<sup>+</sup>, C<sub>2</sub>H<sub>6</sub><sup>+</sup>, C<sub>2</sub>H<sub>5</sub><sup>+</sup>, C<sub>2</sub>H<sub>4</sub><sup>+</sup>, C<sub>2</sub>H<sub>3</sub><sup>+</sup>, C<sub>2</sub>H<sub>2</sub><sup>+</sup>, C<sub>2</sub>H<sup>+</sup>, C<sub>2</sub><sup>+</sup>, H<sub>2</sub><sup>+</sup>, H<sup>+</sup>, C<sup>+</sup>, C<sup>++</sup>, C<sup>+++</sup>)</td><td>Coulomb-potential</td><td>Eq. (42)</td></tr></tbody></table>

## Acknowledgments

This work has received funding from the French National Research Agency (ANR) under the project ANR CH4Kin (ANR grant agreement ANR-24-CE51-7698), as well as funding from Spark Cleantech and ANRT (CIFRE n°2024/0023).

## References

[1] T. Gül, ‘Global Hydrogen Review 2025’, International Energy Agency, Paris, 2025. Accessed: Nov. 12, 2025. [Online]. Available: https://www.iea.org/reports/global-hydrogen-review-2025

[2] L. Fulcheri, V.-J. Rohani, E. Wyse, N. Hardman, and E. Dames, ‘An energy-efficient plasma methane pyrolysis process for high yields of carbon black and hydrogen’, *International Journal of Hydrogen Energy*, vol. 48, no. 8, pp. 2920–2928, Jan. 2023, doi: 10.1016/j.ijhydene.2022.10.144.

[3] P. D. G. Maqueo, S. Coulombe, and J. M. Bergthorson, ‘Energy efficiency of a nanosecond repetitively pulsed discharge for methane reforming’, *J. Phys. D: Appl. Phys.*, vol. 52, no. 27, p. 274002, May 2019, doi: 10.1088/1361-6463/ab199b.

[4] D. Z. Pai, D. A. Lacoste, and C. O. Laux, ‘Transitions between corona, glow, and spark regimes of nanosecond repetitively pulsed discharges in air at atmospheric pressure’, *Journal of Applied Physics*, vol. 107, no. 9, p. 093303, 2010, doi: 10.1063/1.3309758.

[5] N. Minesi, S. Stepanyan, P. Mariotto, G. D. Stancu, and C. O. Laux, ‘Fully ionized nanosecond discharges in air: the thermal spark’, *Plasma Sources Sci. Technol.*, vol. 29, no. 8, p. 085003, Aug. 2020, doi: 10.1088/1361-6595/ab94d3.

[6] S. Delahaie, P.-A. Goutier, N. Q. Minesi, S. D. McGuire, E. Pannier, and C. O. Laux, ‘Experimental characterization of nanosecond discharges in methane’, presented at the 2026 AIAA SciTech Forum (Submitted), Jan. 2026.

[7] E. Roger, P. Mariotto, and C. O. Laux, ‘Origin of the recirculation flow pattern induced by nanosecond discharges and criterion for its development’, *J. Phys. D: Appl. Phys.*, Feb. 2025, doi: 10.1088/1361-6463/adb1f0.

[8] D. Z. Pai, D. A. Lacoste, and C. O. Laux, ‘Nanosecond repetitively pulsed discharges in air at atmospheric pressure—the spark regime’, *Plasma Sources Sci. Technol.*, vol. 19, no. 6, p. 065015, Nov. 2010, doi: 10.1088/0963-0252/19/6/065015.

[9] J. Maillard, ‘Energy Efficiency of CO2 Conversion with Nanosecond Plasma Discharges’, PhD Thesis, Université Paris-Saclay, CentraleSupélec, 2024.

[10] N. A. Popov, ‘Kinetics of plasma-assisted combustion: effect of non-equilibrium excitation on the ignition and oxidation of combustible mixtures’, *Plasma Sources Sci. Technol.*, vol. 25, no. 4, p. 043002, Aug. 2016, doi: 10.1088/0963-0252/25/4/043002.

[11] N. Minesi, P. Mariotto, E. Pannier, G. D. Stancu, and C. O. Laux, ‘The role of excited electronic states in ambient air ionization by a nanosecond discharge’, *Plasma Sources Sci. Technol.*, vol. 30, no. 3, p. 035008, Mar. 2021, doi: 10.1088/1361-6595/abe0a3.
[12] E. Meeks, J. F. Grcar, R. J. Kee, and H. K. Moffat, 'AURORA: A FORTRAN program for modeling well stirred plasma and thermal reactors with gas and surface reactions', SAND--96-8218, 206570, Feb. 1996. doi: 10.2172/206570.

[13] U. S. Inan, A. S. Inan, and R. K. Said, *Engineering electromagnetics and waves*, 2. ed., Global ed. in Always learning. Boston, Mass.: Pearson, 2016.

[14] M. Mitchner and C. H. Kruger, *Partially ionized gases*. in Wiley series in plasma physics. New York NY: Wiley, 1973.

[15] D. G. Goodwin, H. K. Moffat, I. Schoegl, R. L. Speth, and B. W. Weber, *Cantera: An Object-oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes*. (Dec. 14, 2024). Zenodo. doi: 10.5281/ZENODO.14455267.

[16] P. Virtanen et al., 'SciPy 1.0: fundamental algorithms for scientific computing in Python', *Nat Methods*, vol. 17, no. 3, pp. 261-272, Mar. 2020, doi: 10.1038/s41592-019-0686-2.

[17] M. Scapinello, E. Delikonstantis, and G. D. Stefanidis, 'Direct methane-to-ethylene conversion in a nanosecond pulsed discharge', *Fuel*, vol. 222, pp. 705-710, June 2018, doi: 10.1016/j.fuel.2018.03.017.

[18] Y. Gao, S. Zhang, H. Sun, R. Wang, X. Tu, and T. Shao, 'Highly efficient conversion of methane using microsecond and nanosecond pulsed spark discharges', *Applied Energy*, vol. 226, pp. 534-545, Sept. 2018, doi: 10.1016/j.apenergy.2018.06.006.

[19] R. Lotfalipour, A. M. Ghorbanzadeh, and A. Mahdian, 'Methane conversion by repetitive nanosecond pulsed plasma', *J. Phys. D: Appl. Phys.*, vol. 47, no. 36, p. 365201, Aug. 2014, doi: 10.1088/0022-3727/47/36/365201.

[20] E. Morais, E. Delikonstantis, M. Scapinello, G. Smith, G. D. Stefanidis, and A. Bogaerts, 'Methane coupling in nanosecond pulsed plasmas: Correlation between temperature and pressure and effects on product selectivity', *Chemical Engineering Journal*, vol. 462, p. 142227, Apr. 2023, doi: 10.1016/j.cej.2023.142227.

[21] E. Delikonstantis, M. Scapinello, O. Van Geenhoven, and G. D. Stefanidis, 'Nanosecond pulsed discharge-driven non-oxidative methane coupling in a plate-to-plate electrode configuration plasma reactor', *Chemical Engineering Journal*, vol. 380, p. 122477, Jan. 2020, doi: 10.1016/j.cej.2019.122477.

[22] J. R. Fincke, R. P. Anderson, T. A. Hyde, and B. A. Detering, 'Plasma Pyrolysis of Methane to Hydrogen and Carbon Black', *Ind. Eng. Chem. Res.*, vol. 41, no. 6, pp. 1425-1435, Mar. 2002, doi: 10.1021/ie010722e.

[23] B. J. McBride, M. J. Zehe, and S. Gordon, 'NASA Glenn coefficients for calculating thermodynamic properties of individual species', National Aeronautics and Space Administration, Technical report NASA/TP-2002-211556, 2002. [Online]. Available: https://ntrs.nasa.gov/citations/20020085330

[24] M. Gautier, 'Etude de la formation de nanoparticules de carbone au cours de la décomposition thermique d'hydrocarbures : application à la coproduction de noir de carbone et d'hydrogène par craquage thermique du méthane par voie plasma', PhD Thesis, Université Paris Sciences et Lettres, MINES ParisTech, 2016.

[25] S. Heijkers, M. Aghaei, and A. Bogaerts, 'Plasma-Based CH₄ Conversion into Higher Hydrocarbons and H₂: Modeling to Reveal the Reaction Mechanisms of Different Plasma Sources', *J. Phys. Chem. C*, vol. 124, no. 13, pp. 7016-7030, Apr. 2020, doi: 10.1021/acs.jpcc.0c00082.

[26] E. Goos, A. Burcat, and B. Ruscic, 'Extended Third Millennium Thermodynamic Database of New NASA Polynomials with Active Thermochemical Tables update.' http://burcat.technion.ac.il/dir, 2023. doi: 10.2172/925269.

[27] J. Lym, G. R. Wittreich, and D. G. Vlachos, 'A Python Multiscale Thermochemistry Toolbox (pMuTT) for thermochemical and kinetic parameter estimation', *Computer Physics Communications*, vol. 247, p. 106864, Feb. 2020, doi: 10.1016/j.cpc.2019.106864.

[28] A. Kramida and Y. Ralchenko, 'NIST Atomic Spectra Database, NIST Standard Reference Database 78'. National Institute of Standards and Technology, 1999. doi: 10.18434/T4W30F.

[29] Y. Wu et al., 'Calculation of 2-temperature plasma thermo-physical properties considering condensed phases: application to CO₂-CH₄ plasma: part 1. Composition and thermodynamic properties', *J. Phys. D: Appl. Phys.*, vol. 49, no. 40, p. 405203, Oct. 2016, doi: 10.1088/0022-3727/49/40/405203.

[30] M. Scapinello, E. Delikonstantis, and G. D. Stefanidis, 'A study on the reaction mechanism of non-oxidative methane coupling in a nanosecond pulsed discharge reactor using isotope analysis', *Chemical Engineering Journal*, vol. 360, pp. 64-74, Mar. 2019, doi: 10.1016/j.cej.2018.11.161.

[31] X. Yuan et al., 'A kinetic study of nonthermal plasma pyrolysis of methane: Insights into hydrogen and carbon material production', *Chemical Engineering Journal*, vol. 499, p. 156396, Nov. 2024, doi: 10.1016/j.cej.2024.156396.

[32] J. A. Manion, 'NIST Chemical Kinetics Database'. National Institute of Standards and Technology. [Online]. Available: https://kinetics.nist.gov/kinetics

[33] T. J. Millar, C. Walsh, M. Van De Sande, and A. J. Markwick, 'The UMIST Database for Astrochemistry 2022', *A&A*, vol. 682, p. A109, Feb. 2024, doi: 10.1051/0004-6361/202346908.

[34] S. Pancheshnyi et al., 'The LXCat project: Electron scattering cross sections and swarm parameters for low temperature plasma modeling', *Chemical Physics*, vol. 398, pp. 148-153, Apr. 2012, doi: 10.1016/j.chemphys.2011.04.020.

[35] L. C. Pitchford et al., 'LXCat: an Open-Access, Web-Based Platform for Data Needed for Modeling Low Temperature Plasmas', *Plasma Processes & Polymers*, vol. 14, no. 1-2, p. 1600098, Jan. 2017, doi: 10.1002/ppap.201600098.

[36] E. Carbone et al., 'Data Needs for Modeling Low-Temperature Non-Equilibrium Plasmas: The LXCat Project, History, Perspectives and a Tutorial', *Atoms*, vol. 9, no. 1, p. 16, Feb. 2021, doi: 10.3390/atoms9010016.

[37] J. Tennyson et al., 'QDB: a new database of plasma chemistries and reactions', *Plasma Sources Sci. Technol.*, vol. 26, no. 5, p. 055014, Apr. 2017, doi: 10.1088/1361-6595/aa6669.
[38] R. K. Janev and D. Reiter, 'Collision processes of CHy and CHy+ hydrocarbons with plasma electrons and protons', *Physics of Plasmas*, vol. 9, no. 9, pp. 4071-4081, Sept. 2002, doi: 10.1063/1.1500735.

[39] R. K. Janev and D. Reiter, 'Collision processes of C2,3Hy and C2,3Hy+ hydrocarbons with electrons and protons', *Physics of Plasmas*, vol. 11, no. 2, pp. 780-829, Feb. 2004, doi: 10.1063/1.1630794.

[40] 'Biagi database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net

[41] L. L. Alves, 'The IST-LISBON database on LXCat', *J. Phys.: Conf. Ser.*, vol. 565, p. 012007, Dec. 2014, doi: 10.1088/1742-6596/565/1/012007.

[42] 'IST Lisbon database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net

[43] 'Morgan database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net

[44] O. Zatsarinny, 'BSR: B-spline atomic R-matrix codes', *Computer Physics Communications*, vol. 174, no. 4, pp. 273-356, Feb. 2006, doi: 10.1016/j.cpc.2005.10.006.

[45] O. Zatsarinny and K. Bartschat, 'The B-spline R-matrix method for atomic processes: application to atomic structure, electron collisions and photoionization', *J. Phys. B: At. Mol. Opt. Phys.*, vol. 46, no. 11, p. 112001, June 2013, doi: 10.1088/0953-4075/46/11/112001.

[46] 'BSR database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net

[47] H. Tawara and T. Kato, 'Total and partial ionization cross sections of atoms and ions by electron impact', *Atomic Data and Nuclear Data Tables*, vol. 36, no. 2, pp. 167-353, Mar. 1987, doi: 10.1016/0092-640X(87)90014-3.

[48] 'QUANTEMOL database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net

[49] V. Wakelam et al., 'A KINETIC DATABASE FOR ASTROCHEMISTRY (KIDA)', *ApJS*, vol. 199, no. 1, p. 21, Mar. 2012, doi: 10.1088/0067-0049/199/1/21.

[50] C. Niu et al., 'Calculation of 2-temperature plasma thermo-physical properties considering condensed phases: application to CO₂-CH₄ plasma: part 2. Transport coefficients', *J. Phys. D: Appl. Phys.*, vol. 49, no. 40, p. 405204, Oct. 2016, doi: 10.1088/0022-3727/49/40/405204.

[51] M. Capitelli and J. N. Bardsley, Eds., *Nonequilibrium Processes in Partially Ionized Gases*, vol. 220. in NATO ASI Series, vol. 220. Boston, MA: Springer US, 1990. doi: 10.1007/978-1-4615-3780-9.

[52] L. C. Pitchford, B. V. McKoy, A. Chutjian, and S. Trajnar, Eds., *Swarm Studies and Inelastic Electron-Molecule Collisions: Proceedings of the Meeting of the Fourth International Swarm Seminar and the Inelastic Electron-Molecule Collisions Symposium, July 19-23, 1985, Tahoe City, California, USA*. New York, NY: Springer New York, 1987. doi: 10.1007/978-1-4612-4662-6.

[53] 'Hayashi database'. Accessed: Dec. 16, 2024. [Online]. Available: www.lxcat.net