# An engineering model for temperature and pressure evolution during and after nanosecond discharges

Edouard Roger, Pierre Mariotto, Christophe Laux

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

## To cite this version:

Edouard Roger, Pierre Mariotto, Christophe Laux. An engineering model for temperature and pressure evolution during and after nanosecond discharges. Journal of Physics D: Applied Physics, 2026, 10.1088/1361-6463/ae39e6. hal-05474737

HAL Id: hal-05474737

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

Submitted on 23 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An Engineering Model for Temperature and Pressure Evolution during and after Nanosecond Discharges

Edouard Roger¹ᵃ, Pierre Mariotto¹, Christophe O. Laux¹

¹ Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 3 rue Joliot-Curie, 91190 Gif-sur-Yvette, France

**Abstract.** A nanosecond voltage pulse applied between two pin electrodes creates an elongated plasma kernel at high temperature and high pressure. The kernel relaxes through a series of flow dynamics and thermal effects, which affect the temperature and pressure in the inter-electrode region. This work presents an engineering model that predicts the temperature and pressure changes throughout (i) the pulse and (ii) the interpulse period. The pulse model returns the temperature and pressure after the pulse, depending on the discharge regime and the governing parameters (pulse energy, nature of the gas, gap length, and discharge radius). The interpulse model is a sequence of four simple submodels that account for the three dominant flow mechanisms identified in CFD simulations and for vibrational relaxation. The first phase involves the rapid expansion of a blast wave. At the end of this phase, the temperature of the core has been about halved, the pressure has relaxed to the ambient pressure, and the gas has acquired vorticity in a narrow region close to the electrodes. However, this vorticity has not yet had time to disturb the core of the kernel. The second phase corresponds to the vibrational relaxation of the gas towards thermal equilibrium in the case of non-thermal sparks in molecular gases. The third phase is the cooling of the kernel by radial heat diffusion. The fourth phase corresponds to the development of a large recirculation zone from the localized vorticity seeds deposited during the first phase. Two vortices entrain fresh gases into the inter-electrode region. The interpulse model is successfully validated by comparison with a large database of post-discharge CFD simulations with varying inter-electrode gaps and kernel temperatures. In addition, the model shows good agreement with the cumulative heating shown by the temporal temperature profiles measured in literature experiments. It is computationally efficient (a few seconds vs. $10^3–10^6$ CPU hours) and can be incorporated into zero-dimensional kinetic solvers to study the chemistry of the post-discharge. The model can also be used to get criteria on the discharge configuration (temperature, gap, frequency) that ensures full cooling between successive pulses.

E-mail: edouard.roger@centralesupelec.fr, pierre.mariotto@centralesupelec.fr, christophe.laux@centralesupelec.fr

## I. Introduction

Nanosecond Repetitively Pulsed (NRP) discharges are a promising solution to increase the efficiency of chemical processes in many industrial applications. For instance, they are beneficial in dissociating CO₂ into CO [1–3], a fundamental step to synthesize carbon-neutral fuels by the Fischer-Tropsch process [4] or H₂ through the water-gas-shift reaction. NRP discharges are also developed to assist combustion in order to improve

<footnote>ᵃ Author to whom any correspondence should be addressed.</footnote>
ignition [5], control thermoacoustic instabilities [6–8], or extend Lean Blow Out (LBO) limits [9–11] at a low energy cost. Such applications are promising to reduce the environmental footprint of combustion systems.

NRP discharges are typically generated by short-duration (~10 ns) pulses of a few mJ applied between two pin electrodes separated by a few mm. These pulses are repeated at a frequency of tens of kHz. This pattern is illustrated in Figure 1 and can be divided into two phases: the spark phase and the interpulse phase. Both phases are described below.

<figure><img src="image_3.png" /><figcaption>Figure 1. Applied voltage during NRP discharges. Taken from [12].</figcaption></figure>

The pulse phase corresponds to the application of a voltage pulse across the electrodes during $\tau_{pulse}$ (see Figure 1). The voltage pulse induces an electric field in the inter-electrode region. Thus, electrons accelerate and transfer their kinetic energy to the gas by excitation and ionization. This leads to ultra-fast heating and the production of a plasma with radicals and charged species in an elongated volume between the electrodes. Several pulse regimes can occur, namely the glow regime [13], the non-thermal spark [14], and the thermal spark [15] in order of increasing deposited energy. Right after the pulse, the kernel can reach high pressures $P_k$ (typically 3–100 atm) and high temperatures $T_k$ (typically 1000–30,000 K) compared to the ambient gases, which are typically at $P_{out} \sim 1$ atm and $T_{out} \sim 300–2000$ K. The kernel formation is approximately isochoric because it occurs on time scales shorter than the rarefaction time scale (also called hydrodynamics time scale).

The interpulse phase corresponds to the time interval between the end of the pulse and the beginning of the following pulse ($0 < t < \tau_{interpulse}$). In a previous article [12], we presented a detailed study of the flow dynamics phenomena during the interpulse. We
identified three phases. First, a leading shock leaves the plasma kernel and generates vorticity near the electrodes via the "Instantaneous Shock-Induced Vorticity" (ISIV) mechanism. ISIV leads to the formation of vortices near the tips of the electrode. Subsequently, a rarefaction wave propagates towards the center of the kernel. The pressure relaxes to ambient conditions, and the kernel temperature is approximately halved. Second, the kernel cools mainly by thermal diffusion, and the vortices slow down due to viscous effects. Third, the vortices move towards the middle of the inter-electrode gap due to self-advection (unless they vanish because of viscous dissipation). This flow pattern leads to the radial ejection of hot gases along the symmetry plane [12,16–19]. In addition to these three effects of flow dynamics, energy relaxation towards thermal equilibrium may occur. This is the case for molecular gases that would have been excited by a non-thermal spark up to a non-equilibrium state.

The efficiency of heat and radical production by NRP discharges depends on key parameters such as the peak voltage $V_{max}$, the pulse duration $\tau_{pulse}$, the repetition frequency $f = 1/\tau_{interpulse}$, and the inter-electrode gap length $L$. Optimizing these parameters requires simulations of the pulse phase, the interpulse phase, and their cumulative effects over successive pulses. For the pulse phase, efficient zero-dimensional (0D) simulation tools, usually coupling the kinetic equations with a Boltzmann solver for electrons, are already available [20–22]. In contrast, numerical tools to simulate the interpulse phase and the multipulse cumulative effects are less available. Most studies on the interpulse phase are based on multidimensional CFD simulations that solve the reactive Navier-Stokes equations [23–25]. Although this method accounts for the multiphysics of NRP discharges, they are time-consuming, which limits the number of simulations that can be performed and makes it difficult to optimize the discharge parameters ($V_{max}$, $\tau_{pulse}$, $f$, $L$). In addition, the number of species and reactions is limited by the computational cost (typically ~10–100 species and ~10–1000 reactions) [26,27], which prevents the use of the most advanced state-to-state plasma kinetic mechanisms (~100–10,000 species in ground or excited states, and 100,000s of reactions) [28–32].

For these reasons, a computationally efficient engineering model of nanosecond discharges is needed to study the coupled effects of flow dynamics and plasma kinetics. The present work extends the results presented in [33], which focused on the discharge volume and temperature growth and decay during the pulse and blast wave phases.

This article presents a comprehensive model that can predict the temporal evolution of temperature and pressure during and after the pulse for different gases. These laws can be used in various ways later on. First, they can be input into a 0D chemistry solver to predict the evolution of species densities. Second, they can be used to design a discharge, for example, to control the degree of cooling between two pulses. Third, they can be used through model iterations to predict the multipulse cumulative heating between the electrodes in the case of repetitive pulses.
The complete model combines a model for the pulse phase and another model for the interpulse, as illustrated in Figure 2. First, in Section II, we introduce the pulse model, which determines the temperature $T_k$ and pressure $P_k$ immediately after the pulse for both the non-thermal and thermal spark regimes. Second, we present the interpulse model, which takes $T_k$ and $P_k$ as inputs and predicts the variations in temperature $T(t)$ and pressure $P(t)$ of the hot gases (Section III for CFD simulations, Sections IV, V, and VI for model development, validation, and extension, respectively). In Section VII, we utilize the interpulse model to derive simple relationships that establish the fastest cooling mechanism and determine whether it will have sufficient time to return the plasma to ambient conditions before the next pulse begins. In Section VIII, we test the pulse and interpulse models against two experimental datasets. On the one hand, we select an experiment from the literature performed in ambient air and compare the predictions of our model to the temperature profile of the kernel during both the pulse and the interpulse. On the other hand, we combine the pulse and interpulse models iteratively to study the cumulative heating and test the coupled models against two different Plasma-Assisted Combustion (PAC) experiments.

Figure 2 illustrates the structure of the article to improve readability. Furthermore, for each model, the different assumptions made are labeled (a1), (a2), and so on, to enhance clarity. For readers interested in implementing the various models described in this article, we also provide a summary of the assumptions and equations of the engineering model in Appendix B.
<figure><img src="image_4.png" /><figcaption>Figure 2. Structure of the article.</figcaption></figure>
## II. Pulse model

This Section presents the engineering model for the pulse phase (see Figure 2). It is tailored for the two spark regimes commonly encountered in nanosecond discharges: the non-thermal [14] and the thermal spark regimes [15]. For each case, we establish a set of formulas to calculate the trans-rotational temperature $T_k$ and pressure $P_k$ of the kernel after the pulse. We also calculate the vibrational temperature $T_k^ν$ after the pulse in the case of non-thermal sparks of molecular gases. This model solely relies on the pulse energy, thermodynamic properties of the interelectrode gas prior to the pulse ($T_{ini}$, $P_{ini}$), and the kernel dimensions.

### II.A. Non-thermal sparks

The pulse model for non-thermal sparks is inspired by the model of Castela *et al.* [25], initially built for air discharges in order to compute the production of radicals and fast heating in plasma-assisted combustion. In that model, the pulse energy $\Delta E_{pulse}$ deposited into the plasma is decomposed into different fractions. A fraction $\theta_{uh}$ of the energy goes to trans-rotational ultrafast heating (Eqn. (1)). Another fraction $\theta_{sh}$ of the pulse energy<fnref n="2" /> goes to vibrational excitation in the case of molecular gases (the subscripts “uh” and “sh” stand for ultrafast and slow heating). The rest of the energy goes to ionization and molecular dissociation. After the pulse, the plasma is neither in thermal equilibrium nor in chemical equilibrium. Equilibrium is only reached during the interpulse, on a timescale much longer than the pulse duration.

In this pulse model, we first focus on the translational temperature (assumed equal to the rotational temperature), so we only consider the energy $\theta_{uh} \times \Delta E_{pulse}$. We make the assumption (a1) that the ultrafast heating is assumed to be isochoric because it occurs on a timescale shorter than the acoustic one [28,34,35]. The acoustic (or rarefaction) timescale can be estimated by dividing the discharge radius by the speed of sound in the hot gases. The acoustic timescale represents the time at which the head of the rarefaction wave, traveling at the speed of sound, reaches the center of the discharge and initiates the decrease of the temperature and pressure within the core of the plasma (see Section 7.7 of [36] for more details). In the case of non-thermal sparks, the kernel temperature is below 6000 K, and its radius is at least 100 µm. Thus, we can minimize the acoustic timescale by 66 ns for discharges in N₂, 88 ns for discharges in CO₂, and 70 ns for discharges in O₂. These values are all higher than the typical duration of ultrafast heating processes (10 – 20 ns). Thus, the kernel density is constant and equal to the density of the gases between the electrodes before the pulse $\rho_{ini}$. We similarly define the initial temperature, pressure, density, and molar mass $T_{ini}$, $P_{ini}$, $\rho_{ini}$, $M_{ini}$. The translational

<footnote n="2">² $\theta_{uh}$, $\theta_{sh}$ are equivalent to $\alpha - \eta$ and $1 - \alpha$ with the notations of [25], but $\alpha$ will be used in this article as thermal diffusivity.</footnote>
temperature $T_k$ of the gas at the end of the pulse can be derived from the analytical formula given by Eqn. (2):

$$ \theta_{\text{uh}} \times \Delta E_{\text{pulse}} = \pi R_d^2 L \times \rho_{\text{ini}} \times c_v^{\text{tr-rot}} (T_k - T_{\text{ini}}) \quad (1) $$

$$ T_k = T_{\text{ini}} + \frac{\theta_{\text{uh}} \times \Delta E_{\text{pulse}}}{\pi R_d^2 L \times \rho_{\text{ini}} \times c_v^{\text{tr-rot}}} \quad (2) $$

where $\rho_{\text{ini}} = P_{\text{ini}}M_{\text{ini}}/(RT_{\text{ini}})$, $R_d$ is the mean kernel radius, $R$ the ideal gas constant, $c_v^{\text{tr-rot}}$ the specific trans-rotational heat capacity given by $c_v^{\text{tr-rot}} = (n_{\text{freedom}}/2)R/M$ for a gas with $n_{\text{freedom}}$ trans-rotational degrees of freedom (3 for atoms, 5 for diatomic or linear molecules, 6 in other cases). The change of molar mass by dissociation is negligible during non-thermal sparks. Thus, we can assume (a2) that the molar mass $M$ is equal to the ambient molar mass $M_{\text{ini}}$ in the formula of $c_v^{\text{tr-rot}}$.

The pressure of the plasma kernel is calculated by combining the isochoric heating assumption, the ideal gas law, and the principle of mass conservation:

$$ P_k = P_{\text{ini}} \frac{M_{\text{ini}} T_k}{M_k T_{\text{ini}}} \quad (3) $$

where $M_k$ is the molar mass of the plasma kernel immediately after the pulse. Once again, in the usual cases of non-thermal sparks, the dissociation and ionization degrees are low, and we have approximately $M_k = M_{\text{ini}}$. Eqn. (3) simplifies to:

$$ P_k = P_{\text{ini}} \frac{T_k}{T_{\text{ini}}} \quad (4) $$

We now focus on the potential vibrational heating during the pulse, and we calculate, with Eqn. (5), the vibrational temperature $T_k^v$ of the gas at the end of the pulse using the same assumptions as above:

$$ \frac{\theta_{\text{sh}} \times \Delta E_{\text{pulse}}}{\pi R_d^2 L \times \rho_{\text{ini}}} = \int_{T_{\text{ini}}}^{T_k^v} c_v^{\text{vib}}(T) dT \quad (5) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. (a) Specific heat capacity (normalized by its maximal value) and (b) vibrational heat energy as a function of the gas temperature (normalized by the characteristic vibrational temperature of the gas).</figcaption></figure>

We make an additional assumption on $c_v^{\text{vib}}(T)$ to calculate $T_k^v$ explicitly (which is the essence of our engineering model). (a3) We consider that the vibrational specific heat capacity is linear for temperatures below the characteristic vibrational temperature $\Theta_{\text{vib}}$ and constant at $c_v^{\text{vib,max}}$ above (see Table 18.2 in [37] for typical values of $\Theta_{\text{vib}}$):

$$c_v^{\text{vib}}(T) = \begin{cases} c_v^{\text{vib,max}} T / \Theta_{\text{vib}} & \text{if } T \le \Theta_{\text{vib}} \\ c_v^{\text{vib,max}} & \text{if } T > \Theta_{\text{vib}} \end{cases} \qquad (6)$$

This model is particularly well-suited for diatomic molecules. Figure 3a compares the $c_v^{\text{vib}}$ model with the heat capacity of a typical diatomic gas (see Section 18 in [37]). Figure 3b compares the corresponding vibrational energies. The differences are largest when $T/\Theta_{\text{vib}} < 0.2$, but in that case $c_v^{\text{vib}}(T)$ is negligibly small. The error induced by the approximation of the heat capacity in the estimation of $T_k^v$ is shown in Figure 4. For a given vibrational energy, we calculate the vibrational temperatures returned by the two heat capacities and plot the relative error. The error is always below 30%, and below 10% when the vibrational heating exceeds $0.7 \times \Theta_{\text{vib}}$.

<figure><img src="image_6.png" /><figcaption>Figure 4. Relative error on temperature determination using the linear-constant model (Eqn. (6)) of specific heat capacity compared to the typical diatomic gas model.</figcaption></figure>
Using the simplified heat capacity model, we can then compute $T_k^v$ starting from Eqn. (5):

$$ T_k^v = \begin{cases} \frac{\Delta e_{\text{sh}}}{c_v^{\text{vib,max}}} + T_{\text{ini}} & \text{if } T_{\text{ini}} \ge \Theta_{\text{vib}} \\ \frac{\Delta e_{\text{sh}}}{c_v^{\text{vib,max}}} + \frac{\Theta_{\text{vib}}}{2} + \frac{T_{\text{ini}}^2}{2\Theta_{\text{vib}}} & \text{if } T_{\text{ini}} < \Theta_{\text{vib}} \text{ and } \frac{\Delta e_{\text{sh}}}{c_v^{\text{vib,max}}} > \frac{\Theta_{\text{vib}}^2 - T_{\text{ini}}^2}{2\Theta_{\text{vib}}} \\ \sqrt{\frac{2\Theta_{\text{vib}}\Delta e_{\text{sh}}}{c_v^{\text{vib,max}}} + T_{\text{ini}}^2} & \text{otherwise} \end{cases} \quad (7) $$

with $\Delta e_{\text{sh}} \equiv \theta_{\text{sh}} \Delta E_{\text{pulse}} / (\pi R_d^2 L \rho_{\text{ini}})$.

In summary, the engineering model for non-thermal sparks takes $\theta_{\text{uh}}$, $\Delta E_{\text{pulse}}$, $R_d$, $L$, $T_{\text{ini}}$, $P_{\text{ini}}$, $\rho_{\text{ini}}$ as inputs and returns $P_k$ and $T_k$ from Eqns. (2) and (4). In the case of molecular gases, the vibrational temperature of the kernel, $T_k^v$, is obtained from Eqn. (7) using $\theta_{\text{sh}}$, $\Theta_{\text{vib}}$, $c_v^{\text{vib,max}}$ as inputs.

## II.B. Thermal sparks

Thermal sparks have been studied in detail by Minesi et al. [15] and Maillard et al. [38]. They showed that, at the end of a pulse with sufficient energy, the plasma kernel can reach thermal equilibrium conditions, and sometimes even chemical equilibrium conditions (i.e., Local Thermodynamic Equilibrium, LTE), at temperatures typically between 20,000 and 50,000 K. In such cases, all deposited energy is converted into internal energy (translational, rotational, vibrational, electronic, and chemical energy) during the pulse.

An assumption (a1) regarding the chemical state is required to calculate $T_k$ and $P_k$. We consider two limit cases to determine an upper and lower value for these quantities. On the one hand, we assume that the gas chemistry is infinitely fast, i.e., the gas is at LTE. On the other hand, we assume that the chemistry is frozen. For these two scenarios, the conservation of energy gives:

$$ \Delta E_{\text{pulse}} = \begin{cases} \pi R_d^2 L \times \rho_{\text{ini}} \times (u_{\text{eq}}(T_k, P_k) - u_{\text{eq}}(T_{\text{ini}}, P_{\text{ini}})) & (\text{LTE}) \\ \pi R_d^2 L \times \rho_{\text{ini}} \times (u_f(T_k, P_k) - u_{\text{ini}}) & (\text{frozen}) \end{cases} \quad (8) $$

where $u_{\text{eq}}(T, P)$ is the specific internal energy of the plasma kernel in LTE, $u_f(T, P)$ is the specific internal energy of the plasma kernel whose composition is frozen to that at the pre-pulse initial conditions ($T_{\text{ini}}$, $P_{\text{ini}}$), and $u_{\text{ini}}$ is the initial internal specific energy of the kernel. The specific LTE internal energy reads $u_{\text{eq}}(T, P) = h_{\text{eq}}(T, P) - P/\rho_{\text{eq}}(T, P)$ where $h_{\text{eq}}$ and $\rho_{\text{eq}}$ are the specific enthalpy and density of the gas in LTE.
$h_{eq}$ and $\rho_{eq}$ are taken from Chapter 10 of Ref. [39]. $u_f(T, P)$ and $u_{ini}$ are computed using the code CANTERA [40] and the NASA9 coefficients given in [41].

During the formation of the thermal spark, (a2) gas heating is isochoric because it occurs on a timescale shorter than the acoustic time [15,34]. Indeed, the kernel temperature is typically one order of magnitude higher than for non-thermal sparks, and the radius of the discharge is one order of magnitude lower. The molar mass is divided by a factor of 4 to 6 due to dissociation and ionization (for diatomic and triatomic gases, respectively), and the adiabatic index is 5/3. Thus, the minimum acoustic timescales are about one nanosecond for discharges in $N_2$, $CO_2$, or $O_2$. In contrast, the heating timescale for thermal sparks is about 100 ps [42]. The criterion for isochoricity is therefore ensured, and the density remains constant:

$$ \rho_{eq}(T_k, P_k) = \frac{P_k M_k}{RT_k} = \rho_{ini} \qquad (9) $$

Consequently, $T_k$ and $P_k$ can be determined by solving Eqns. (8)–(9) with a two-dimensional root-finding solver. Contrary to the non-thermal spark case, the plasma’s molar mass may change due to dissociation and ionization. We emphasize that Eqn. (9) has no simple analytical solution $P_k$ that could be inserted into Eqn. (8) because $M_k$ varies with the composition, which may depend on $T_k$ and $P_k$. Therefore, there is no closed-form expression for $T_k$ and $P_k$.

In summary, the engineering model for thermal sparks takes $\Delta E_{pulse}$, $R_d$, $L$, $T_{ini}$, $P_{ini}$ as inputs, solves Eqns. (8)–(9) with a root-finding algorithm, and returns $P_k$ and $T_k$.

We solved Eqns. (8)–(9) for thermal sparks in air, $CO_2$, and $H_2$ (note that the values of $h_{eq}$ and $\rho_{eq}$ for these gases are assumed to be the same as those given by [39] for the atmospheres of Earth, Mars, and Jupiter, which are composed of these respective gases by more than 99, 95, and 90 %). Figures 5–7 show the resulting $T_k$ and $P_k$ as a function of $\Delta E_{pulse}/(\pi R_d^2 L)$ for different initial conditions ($T_{ini}$, $P_{ini}$) and under the assumption that chemistry is either frozen (dashed lines) or infinitely fast (solid lines). $T_k$ is higher with the frozen-chemistry assumption than with the LTE assumption because no input energy is consumed to break molecular bonds or to ionize.

These figures can be used as look-up charts to run the model for thermal sparks in air, $CO_2$, and $H_2$, respectively. They avoid the use of a root-finding algorithm to solve Eqns. (8)–(9).

For the sake of completeness, we mention that a relation between (i) the input lineic energy $\Delta E_{pulse}/L$, (ii) the ambient gas thermodynamic properties ($T_{out}$, $P_{out}$), and (iii) the radius of the discharge $R_d$ is proposed in the literature for strong blast waves (see Refs. [43,44]).
<figure><img src="image_7.png" /><figcaption>Figure 5. Chart to relate the volumetric pulse energy $\Delta E_{pulse}/(\pi R_d^2 L)$ to the temperature $T_k$ and pressure $P_k$ of the gas immediately after the pulse for a thermal spark in air under different ambient conditions. Dashed lines: frozen chemistry. Solid lines: infinitely fast chemistry (LTE).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. Chart to relate the volumetric pulse energy $\Delta E_{pulse}/(\pi R_d^2 L)$ to the temperature $T_k$ and pressure $P_k$ of the gas immediately after the pulse for a thermal spark in CO$_2$ under different ambient conditions. Dashed lines: frozen chemistry. Solid lines: infinitely fast chemistry (LTE).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7. Chart to relate the volumetric pulse energy $\Delta E_{pulse}/(\pi R_d^2 L)$ to the temperature $T_k$ and pressure $P_k$ of the gas immediately after the pulse for a thermal spark in H$_2$ under different ambient conditions. Dashed lines: frozen chemistry. Solid lines: infinitely fast chemistry (LTE).</figcaption></figure>
### III. Interpulse model: CFD simulations

In this Section, we perform several CFD simulations of the post-discharge following a nanosecond pulse. We will use these simulations to validate the interpulse model and its three flow dynamics submodels. The working gas is air, considered in LTE. Unlike the pulse models developed above, this Section (along with Sections IV and V) focuses on the interpulses of non-thermal sparks for the sake of simplicity. We will propose an extension of our interpulse model to thermal sparks in Section VI.

In Subsection III.A, we detail the numerical models we use. In Subsection III.B, we analyze the temperature and pressure evolution inside the kernel and identify several flow phases (expansion, diffusion, and recirculation). These observations will later guide the design of the post-discharge engineering model presented in Section IV.

#### III.A. Numerical model

The fluid dynamics effects caused by the discharge are simulated using the Ansys® Fluent® program, release 22.1.0 [45]. The solvers, boundary conditions, initial conditions, and the mesh are the same as those used in Ref. [12]. In summary, we solve the 2D non-reactive Navier-Stokes equations on the axisymmetric computational domain (pin-to-pin configuration) shown in Figure 8. We do not solve the Maxwell equations because there is no external electromagnetic field during the post-discharge, and because the self-induced electromagnetic field is negligible (low ionization degree below $10^{-3}$). The average cell size is 2 µm and is refined by a factor of 10 in the boundary layers near the electrodes. An adiabatic no-slip boundary condition is applied at the electrode walls. The ambient gas is initially at rest at $T_{\text{out}} = 300$ K and $P_{\text{out}} = 1$ atm. The hot kernel (elliptic red zone in Figure 8) has a radius $R_d = 100$ µm and a gap length $L$. In this kernel, the temperature is initialized at $T_k$ and the pressure at $P_k = T_k \times P_{\text{out}}/T_{\text{out}}$, in line with the non-thermal spark models presented in Section II. In addition, we assume that the gas remains at LTE, because our previous study [12] showed that the choice of the chemical model (frozen, LTE, and finite-rate chemistry) does not significantly alter the flow dynamic predictions of non-thermal sparks in air with kernel temperatures up to $T_k = 6000$ K. We select the thermodynamic and transport properties of equilibrium air computed by Gupta *et al.* [46].

In addition to this code, we use two Python in-house 1D finite-volume solvers locally in the article to solve the Euler and heat equations, respectively (in a cylindrical domain). The first one employs a 7th-order Roe-type Riemann solver [47,48]. The second one uses a 2nd-order central scheme. Without explicit mention of one of these two solvers, the reader should assume that the simulations are performed using the Fluent® program described before within a 2D domain.
<figure><img src="image_10.png" /><figcaption>Figure 8. Computational domain. R_d and L represent the radius and the length of the plasma kernel, respectively. Taken from [12].</figcaption></figure>

### III.B. Identification of three flow phases

This Subsection focuses on two simulations conducted using the numerical model, which will serve as test cases in Subsection V.A. These configurations correspond to the two main post-discharge regimes: the diffusion regime (Case 1) and the recirculation regime (Case 2). The gap and kernel temperature of these two cases are as follows:

* Case 1: L = 1 mm, T_k = 500 K

* Case 2: L = 1 mm, T_k = 5500 K (we use this case as an example throughout Section IV)

For each case and at any instant t, we compute the mean temperature $\bar{T}$ and the mean pressure $\bar{P}$ along the inter-electrode axis with Eqns. (10) and (11):

$$ \bar{T}(t) = \frac{1}{L/2} \int_{z=0}^{z=L/2} T(r=0, z, t) dz \qquad (10) $$

$$ \bar{P}(t) = \frac{1}{L/2} \int_{z=0}^{z=L/2} P(r=0, z, t) dz \qquad (11) $$
The evolution of $\bar{T}$ and $\bar{P}$ are shown in Figure 9 and Figure 10. Three phases can be identified:

1. Blast wave expansion (0.03–0.1 µs): the temperature decreases to approximately half its initial value, and the kernel pressure reaches ambient pressure. This phase is represented by the blue-shaded area in Figures 9–10. This phenomenon is sometimes followed by a subatmospheric depression as well as overcooling due to the blast wind (marked by the dagger symbol). It can also be followed by a jump in pressure and temperature due to a secondary shock (marked by an asterisk). Additional slight perturbations of the pressure and temperature profiles are due to the rebounds of rarefaction waves (marked by a double dagger symbol in both cases).

2. Recirculation (3–30 µs): the kernel collapses into a torus, and hot gases are expelled radially along the symmetry plane. The plasma cools down to the ambient temperature. This phase is represented by the green-shaded area. It occurs in Case 2 (Figure 9b) but not in Case 1 (Figure 9a).

3. Diffusion (10–1000 µs): the kernel begins to cool down by exchanging energy with the environment through diffusion. This phenomenon may occur in conjunction with recirculation, as in Case 2. It is indicated by the gray-shaded area in Figures 9–10.

More details about these three phases can be found in Subsection III.A of Ref. [12].
<figure><img src="image_11.png" /><figcaption>Figure 9. Time history of the average temperature computed along the pin-to-pin axis. (a) Case 1: $L = 1$ mm, $T_k = 500$ K. (b) Case 2: $L = 1$ mm, $T_k = 5500$ K. The dagger (†) marks the blast wind overcooling. The asterisk (*) marks the temperature peak due to the secondary shock. The double dagger (‡) marks the additional perturbation due to the rebounds of rarefaction waves (see text).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. Time history of the average pressure computed along the pin-to-pin axis. (a) Case 1: $L = 1$ mm, $T_k = 500$ K. (b) Case 2: $L = 1$ mm, $T_k = 5500$ K. The dagger (†) marks the blast wind subatmospheric depression. The asterisk (*) marks the pressure peak due to the secondary shock. The double dagger (‡) marks the additional perturbation due to the rebounds of rarefaction waves (see text).</figcaption></figure>
## IV. Interpulse model: modeling with four submodels in series

In this Section, we develop an engineering model that can predict the temperature $\bar{T}$ (Eqn. (10)) and pressure $\bar{P}$ (Eqn. (11)) of the plasma kernel during the interpulse. We aim to capture the intensity and timescale of the different phases, as well as their dependency on the discharge parameters ($L, T_k, P_k, R_d, T_{out}, P_{out}$ and nature of the gas). Subsection IV.A outlines the modeling strategy. Subsections IV.B, IV.C, IV.D, and IV.E present a simplified submodel for each phase. Finally, Subsection IV.F summarizes our whole model.

### IV.A. Principles of the interpulse model

The interpulse model is based on the functional decomposition method. The evolution of the discharge and the flow field is divided into separate phases associated with a subscript $i \in \{1, 2, 3, 4\}$: expansion of the heated channel ($i = 1$), relaxation of the nonequilibrium energy storage ($i = 2$), thermal diffusion ($i = 3$), and gas recirculation ($i = 4$). Phases $i = 1, 3$, and $4$ are related to flow dynamics phenomena. Phase $i = 2$ is related to plasma kinetics and is specific to non-thermal sparks in molecular gases.

We quantitatively analyze and model each phase with a simplified set of equations. For each phase $i$, we derive analytical expressions for the temperature $T_i(t)$ and the pressure $P_i(t)$ of the kernel. These expressions depend on discharge parameters (initial plasma temperature, ambient temperature, heat capacity, thermal diffusivity, kernel radius, adiabatic index). The expressions are intricate, as the expression for step $i$ involves those for step $i-1$. This results in a global interpulse model that takes all phenomena into account.

In parallel to the derivation of these laws, we illustrate them step by step with a sample case: an air discharge with $T_k = 5500$ K, $T_{out} = 300$ K, $P_k = 18.3$ atm, $P_{out} = 1$ atm, $\gamma = 9/7$, $M = 28$ g/mol, and $R_d = 100$ µm. This is the configuration of Case 2 (see Subsection III.B). As the CFD simulations are based on a one-temperature thermodynamic model that does not separate trans-rotational and vibrational energies, it is not possible to test the engineering model of phase $i = 2$ with the CFD results. This will be done later (in Section VIII) through comparison with experiments.

### IV.B. Expansion submodel

In this Subsection, we focus on the first flow dynamics phenomenon that influences the temperature and pressure of the discharge, i.e., the blast wave expansion. Blast wave models were developed in Refs. [43,49–52] using various strategies, but with different input and output parameters than ours. Thus, it would have been challenging to relate
them to our pulse models and the following flow dynamics submodels. For these reasons, we introduce another cylindrical shock wave submodel here.

We make several assumptions to derive closed-form expressions for the evolution of temperature and pressure:

(a1) The early expansion of the shock and the rarefaction wave are 1D in the radial direction, which is particularly relevant for nanosecond filamentary discharges since their radius is small compared to the gap length (the ratio gap/radius is typically about 10).

(a2) The fluid can be considered inviscid compared to compressible effects during the expansion. We calculated the minimum viscous timescale $\tau_{\text{visc}} = R_d^2/\nu$ ($\nu$ being the kinematic viscosity) for argon, air, oxygen, and hydrogen. For all gases, $\tau_{\text{visc}} \ge 10$ µs, which is above the typical acoustic timescale for NTS (100 ns–1 µs).

(a3) The gas is considered calorically perfect, which is approximate, but necessary to derive an analytical solution. We validate this assumption a posteriori.

(a4) The adiabatic index $\gamma$ and molar mass $M$ are constant throughout the blast wave expansion and equal to their values $\gamma_k$ and $M_k$ at the post-pulse thermodynamic conditions. This assumption was assessed in Appendix B of [12], where its impact on the initial Mach number of the shock was shown to be limited. Besides, the impact of dissociation on the molar mass is negligible in NTS.

(a5) The gas is ideal. For example, Gupta et al. [46] show that the compressibility factor $Z$ remains below 1.5 for air at temperatures between 300 and 6000 K and pressures between 0.01 and 10 bar. Hilsenrath et al. [53] provides data for many other gases (argon, carbon dioxide, hydrogen) up to 1000–3000 K, and $Z$ remains between 0.8 and 1.2.

(a6) The gas is initially at rest.

(a7) Heat diffusion is negligible. As we did for (a2), we estimated the minimum heat diffusion timescale $\tau_{\text{heat}} = R_d^2/\alpha$ ($\alpha$ being the thermal diffusivity) for argon, air, and carbon dioxide. We obtained $\tau_{\text{heat}} \ge 5$ µs, which is also above the typical acoustic timescale for NTS (100 ns–1 µs).

(a8) The gas has the following discontinuous temperature and pressure profiles:

$$ T = \begin{cases} T_k & \text{if } 0 \le r \le R_d \\ T_{\text{out}} & \text{otherwise} \end{cases} \qquad (12) $$

$$ P = \begin{cases} P_k & \text{if } 0 \le r \le R_d \\ P_{\text{out}} & \text{otherwise} \end{cases} \qquad (13) $$

The system is then described with the 1D cylindrical Euler equations. The configuration matches that of a blast wave, which is extensively treated in textbooks on compressible fluid dynamics [36]. The specificities of the non-planar configuration are detailed in [54], and the propagation of shock waves in nonequilibrium gases is detailed in [55] or Chapter 16 of [56].
The Euler equations are solved using the 1D in-house code mentioned at the end of Subsection III.A. An example of a solution is shown in Figure 11 with the configuration of Case 2 ($T_k = 5500$ K, $R_d = 100$ µm). It shows the radial profiles of temperature (top left), pressure (top right), density (bottom left), and radial velocity (bottom right) at different instants. Several characteristic features of a cylindrical blast wave are visible. These include the Leading Shock Front (LSF), the Contact Surface between the plasma and the ambient gas (CS), and an inward-propagating Rarefaction Wave (RW). After a hundred nanoseconds, we also see a secondary Shock Front (2^nd^ SF) propagating inward and then reflecting on the center.

<figure><img src="image_13.png" /><figcaption>Figure 11. Radial profiles of the plasma/ambient properties at six instants: temperature (top left), pressure (top right), density (bottom left), radial velocity (bottom right). The main flow dynamics phenomena are indicated with markers, and their propagation directions are shown with arrows: (LSF) leading shock front, (CS) contact surface between the plasma and the ambient gas, (RW) rarefaction wave, (2^nd^ SF) secondary shock front. The initial configuration corresponds to Case 2 ($T_k = 5500$ K, $R_d = 100$ µm).</figcaption></figure>

For the same example, Figure 12 shows the temperature $T(r = 0, t)$ and pressure $P(r = 0, t)$ of the core. The temporal profiles can be decomposed into three main steps. First, the core stays at $T_k$ and $P_k$ until the rarefaction wave reaches the inter-electrode axis ($t \le 0.06$ µs). Second, between $t = 0.06$ µs and $t \approx 0.2$ µs, the rarefaction induces an abrupt fall of temperature and pressure toward a plateau at temperature $T_{rarefaction}$ and pressure $P_{rarefaction}$ close to $P_{out}$. Third, the rebounds of rarefaction waves and/or secondary shocks slightly perturb the temporal profiles. We studied tens of other configurations (varying $\gamma$, $M$, $T_{out}$, $T_k/T_{out}$, $R_d$) and the expansion can always be
divided into these three steps. Thus, they will be the basis of our engineering expansion submodel.

<figure><img src="image_14.png" /><figcaption>Figure 12. Temporal profiles of temperature and pressure at the center of the kernel. The vertical dotted lines indicate the stepwise decomposition used in the expansion submodel. The initial configuration corresponds to Case 2 ($T_k = 5500$ K, $R_d = 100$ µm).</figcaption></figure>

The stepwise submodel of pressure and temperature in the core of the plasma has three unknowns: the temperature plateau ($T_{rarefaction}$), the instant at which the pressure and temperature start to fall ($\tau_{rarefaction}$), and the duration of the abrupt fall ($\delta_{rarefaction}$). To complete the submodel, we need to estimate these terms as a function of the discharge parameters.

The characteristic time $\tau_{rarefaction}$ can be derived from compressible gas theory:

$$ \tau_{rarefaction} \equiv \frac{R_d}{\sqrt{\gamma_k RT_k / \mathcal{M}_k}} \qquad (14) $$

The calculation of the temperature plateau $T_{rarefaction}$ requires an additional assumption. (a9) It is hypothesized that the evolution of the plasma kernel, enclosed within the contact surface and represented in orange in Figure 13, is reversible and adiabatic. Thus, using the isentropic relation $T_{rarefaction}/T_k = (P_{rarefaction}/P_k)^{1-1/\gamma}$ (see Section 1.4.6 of [36]) and given that $P_{rarefaction} \approx P_{out}$ at the end of the rarefaction phase (as demonstrated with the previous 1D simulation), the ideal gas law leads to:

$$ \frac{T_{rarefaction}}{T_{out}} = \left( \frac{T_k}{T_{out}} \right)^{1/\gamma} \qquad (15) $$
<figure><img src="image_15.png" /><figcaption>Figure 13. Volume expansion of the plasma kernel during the rarefaction phase. (a) corresponds to the plasma ideal distribution after the pulse. (b) corresponds to the plasma ideal distribution after the expansion phase.</figcaption></figure>

We identified three non-isentropic processes that may differ from (a9). The secondary shock waves, the heat transfer through the electrodes, and the heat exchange across the contact surface (or within the kernel). The last two are negligible via (a1) and (a7), respectively. Regarding the secondary shock wave, since the temperature $T_{\text{rarefaction}}$ calculated with Eqn. (15) is within 5% of the temperature obtained from the previous 1D CFD simulations (Figure 11) that include secondary shock waves, we conclude that the isentropic approximation (a9) is excellent to estimate $T_{\text{rarefaction}}$. Figure 14a shows the temperature ratio $T_{\text{rarefaction}}/T_{\text{out}}$ as a function of $T_k/T_{\text{out}}$ for different values of $\gamma$.

$\delta_{\text{rarefaction}}$ cannot be derived in closed form. Thus, we performed many 1D CFD simulations, inferred their $\delta_{\text{rarefaction}}$, and fit the dataset with a simple function of temperature. The various simulations employed different values of $T_{\text{out}}$ (within 300–1000 K), $T_k/T_{\text{out}}$ (within 1.05–20), $R_d$ (within 50–300 µm), and gas (N₂, CO₂, CH₄, N, N⁺, with typical molar mass *M* within 2–44 g/mol). For each gas, we also vary $\gamma$ to consider the cases of frozen or infinitely fast vibrational kinetics. Table 1 summarizes the selected range of parameters.
<figure><img src="image_16.png" /><figcaption>Figure 14. (a) Post-rarefaction temperature $T_{\text{rarefaction}}$ as a function of the initial temperature $T_k$, both normalized by the ambient temperature $T_{\text{out}}$. Each color corresponds to a different value of the adiabatic index. (b) Ratio between the temperature/pressure decrease time and the instant of the beginning of this decrease as a function of the plasma-ambient temperature ratio. The colored points are data extracted from solutions of the 1D Euler system of equations; the grey lines are the power-law fits for different values of $\gamma$: 5/3, 7/5, 8/6, 9/7, 15/13, 26/24 (from top to bottom).</figcaption></figure>

For all configurations, we plot in Figure 14b the ratio $\delta_{\text{rarefaction}}/\tau_{\text{rarefaction}}$ as a function of $T_k/T_{\text{out}}$. We note that $M$, $R_d$ and $T_{\text{out}}$ have a negligible impact, unlike $\gamma$ and $T_k/T_{\text{out}}$. We fit $\delta_{\text{rarefaction}}$ with the rational function (16), whose coefficients are $\gamma$-dependent:

$$ \delta_{\text{rarefaction}} = \tau_{\text{rarefaction}} \left( \beta_1 + \frac{(1 - \beta_1)}{1 + \beta_2 \left( \frac{T_k}{T_{\text{out}}} - 1 \right)} \right) \quad \text{with} \quad \begin{cases} \beta_1 = 1.418 + 1.194\gamma \\ \beta_2 = 0.716 - 0.242\gamma \end{cases} \qquad (16) $$
<table><caption>Table 1. Properties of the different gases tested in the 1D CFD simulations of Subsection IV.B and the Global Shock Dynamics (GSD) simulations, and their properties of interest. “w/ vib” and “w/o vib” indicate whether the vibrational degrees of freedom of the molecules stay at equilibrium or are frozen during the expansion.</caption><thead><tr><th>Name of the configuration</th><th>Mean molar mass <i>M</i> [g/mol]</th><th>Adiabatic index <i>γ</i></th></tr></thead><tbody><tr><td>CO<sub>2</sub> w/ vib</td><td>44</td><td>15/13 ≈ 1.154</td></tr><tr><td>CO<sub>2</sub> w/o vib</td><td>44</td><td>7/5 = 1.4</td></tr><tr><td>Ar</td><td>40</td><td>5/3 ≈ 1.667</td></tr><tr><td>N<sub>2</sub> w/ vib</td><td>28</td><td>9/7 ≈ 1.286</td></tr><tr><td>N<sub>2</sub> w/o vib</td><td>28</td><td>7/5 = 1.4</td></tr><tr><td>CH<sub>4</sub> w/ vib</td><td>16</td><td>26/24 ≈ 1.083</td></tr><tr><td>CH<sub>4</sub> w/o vib</td><td>16</td><td>8/6 ≈ 1.333</td></tr><tr><td>N</td><td>14</td><td>5/3 ≈ 1.667</td></tr><tr><td>N<sup>+</sup>/e<sup>-</sup></td><td>7</td><td>5/3 ≈ 1.667</td></tr><tr><td>H<sub>2</sub> w/ vib</td><td>2</td><td>9/7 ≈ 1.286</td></tr><tr><td>H<sub>2</sub> w/o vib</td><td>2</td><td>7/5 = 1.4</td></tr></tbody></table>

The expansion submodel is now complete. Eqns. (17) and (18) summarize the evolution of temperature $T_1(t)$ and pressure $P_1(t)$ as a function of the discharge parameters:

$$ T_1(t) = \begin{cases} T_k & \text{if } t < \tau_{\text{rarefaction}} \\ T_k + (T_{\text{rarefaction}} - T_k) \frac{(t - \tau_{\text{rarefaction}})}{\delta_{\text{rarefaction}}} & \text{if } 1 < \frac{t}{\tau_{\text{rarefaction}}} < 1 + \frac{\delta_{\text{rarefaction}}}{\tau_{\text{rarefaction}}} \\ T_{\text{rarefaction}} & \text{otherwise} \end{cases} \quad (17) $$

$$ P_1(t) = \begin{cases} P_k & \text{if } t < \tau_{\text{rarefaction}} \\ P_k + (P_{\text{out}} - P_k) \frac{(t - \tau_{\text{rarefaction}})}{\delta_{\text{rarefaction}}} & \text{if } 1 < \frac{t}{\tau_{\text{rarefaction}}} < 1 + \frac{\delta_{\text{rarefaction}}}{\tau_{\text{rarefaction}}} \\ P_{\text{out}} & \text{otherwise} \end{cases} \quad (18) $$
IV.C. Submodel for the vibrational relaxation of molecules

In this Subsection we discuss the second phenomenon of our interpulse model, the only one that is not directly related to flow dynamics. After non-thermal sparks in molecular gases, there might be a nonequilibrium distribution of energy between the molecular degrees of freedom. Typically, the vibrational temperature may be higher than the trans- rotational temperature. This results in the transfer of energy between the different energy modes until they reach Thermal Equilibrium (TE) (see Figure 15). During Vibration- Translation (VT) relaxation, the trans-rotational temperature increases. The present submodel is designed to capture this effect.

<figure><img src="image_17.png" /><figcaption>Figure 15. Typical temporal profiles of the translational-rotational and vibrational temperatures in the case of non-thermal sparks in molecular gases. The solid curve stands for the trans-rotational temperature. The dashed line stands for the vibrational temperature.</figcaption></figure>

We make two assumptions in this submodel.

(a1) The typical timescale of VT energy relaxation, $\tau_{VT}$, is significantly higher than the rarefaction timescale, $\tau_{rarefaction}$. This assumption, necessary to decouple the VT relaxation phase from the expansion phase, is justified for many discharge configurations. For instance, in the case of a non-thermal spark in atmospheric pressure air, previous simulations showed that $\tau_{rarefaction} \sim 100$ ns. In [57], the author showed that $\tau_{VT}$ is of the order of a few microseconds, due to the quenching of N$_2$ by atomic oxygen. More generally, Figure 16 gathers values of $\tau_{VT}$ for N$_2$, O$_2$ and CO$_2$
and different colliders. This figure compares models and experimental correlations
from [56–60]. In most cases, $\tau_{VT} \gg 100$ ns. The only exceptions are CO$_2$ discharges
with a translational temperature above 3000 K and O$_2$ discharges with H$_2$ as the main
quencher.

(a2) The characteristic time $\tau_{VT}$ for VT relaxation process is constant. The choice of $\tau_{VT}$ is determined by the user of the model and depends on the composition of the gas, the ambient pressure, and the temperature. We provide Figure 16 to help the user with this choice. We recommend using a temperature close to $T_{\text{rarefaction}}$ and a pressure $P_{\text{out}}$. This assumption is the main limitation of the submodel and may lead to underestimating the duration of the VT relaxation. However, this is necessary to keep the model analytical and easy to use.
<figure><img src="image_18.png" /><figcaption>Figure 16. Characteristic VT relaxation timescales as a function of the temperature, for different couples of molecule (s) – collider (s’). $\tau_{VT}$ is multiplied by the collider partial pressure $P_{s}$, and normalized by $P_{atm} = 101325$ Pa. From top to bottom, s = N₂, O₂, CO₂.</figcaption></figure>
We define $T_{TE}$ as the temperature of the gas after it has reached thermal equilibrium. $T_{TE}$ can be related to the amount of energy deposited into slow heating, $\theta_{sh}\Delta E_{pulse}$, via Eqn. (19):

$$ \frac{\theta_{sh} \times \Delta E_{pulse}}{\pi R_d^2 L \times \rho_{ini}} = c_v^{tr-rot}(T_{TE} - T_{rarefaction}) + \int_{T_{ini}}^{T_{TE}} c_v^{vib}(\theta) d\theta \quad (19) $$

Similar to the non-thermal pulse model (Subsection II.B), we assume that the vibrational specific heat capacity follows Eqn. (6). This allows us to derive the following closed-form expression for $T_{TE}$:

$$ T_{TE} = \begin{cases} A & \text{if } T_{ini} \ge \Theta_{vib} \\ B & \text{if } T_{ini} < \Theta_{vib} \text{ and } K > 0 \\ C & \text{otherwise} \end{cases} \quad (20) $$

with:

$$ A = (\Delta e_{sh} + c_v^{tr-rot} T_{rarefaction} + c_v^{vib,max} T_{ini}) / (c_v^{tr-rot} + c_v^{vib,max}) $$

$$ B = \left( \Delta e_{sh} + c_v^{tr-rot} T_{rarefaction} + \frac{c_v^{vib,max}}{2} \left( \Theta_{vib} + \frac{T_{ini}^2}{\Theta_{vib}} \right) \right) / \left( c_v^{tr-rot} + c_v^{vib,max} \right) $$

$$ C = \sqrt{T_{ini}^2 + \left( \frac{c_v^{tr-rot} \Theta_{vib}}{c_v^{vib,max}} \right)^2 + \frac{2 \Theta_{vib}}{c_v^{vib,max}} (\Delta e_{sh} + c_v^{tr-rot} T_{rarefaction}) - \frac{c_v^{tr-rot} \Theta_{vib}}{c_v^{vib,max}}} $$

$$ K = T_{ini}^2 + \frac{2\Theta_{vib}}{c_v^{vib,max}} (\Delta e_{sh} + c_v^{tr-rot} T_{rarefaction}) - \Theta_{vib}^2 \left( 1 + \frac{2c_v^{tr-rot}}{c_v^{vib,max}} \right) $$

In the first approximation, the evolution of the gas temperature $T(t)$ due to expansion and VT relaxation is modeled by performing a weighted average of $T_{TE}$ and $T_{rarefaction}$:

$$ T(t) = T_{rarefaction} \times (1 - F_{TE}(t)) + T_{TE} \times F_{TE}(t) \quad (21) $$

where $F_{TE}(t)$ and its complementary are the weighting factors. $F_{TE}(t)$ is a simple activation function whose steepness involves $\tau_{rarefaction}$, $\delta_{rarefaction}$ and $\tau_{VT}$:

$$ F_{TE}(t) = \begin{cases} 0 & \text{if } t < \tau_{rarefaction} + \delta_{rarefaction} \\ 1 - \exp\left(-\frac{t - (\tau_{rarefaction} + \delta_{rarefaction})}{\tau_{VT}}\right) & \text{otherwise} \end{cases} \quad (22) $$
Taking advantage of the time splitting between the different phenomena, we substitute $T_{\text{rarefaction}}$ in Eqn. (21) with $T_1(t)$ from Eqn. (17). Thus, the temperature evolution $T_2(t)$ accounting for expansion and VT relaxation reads:

$$T_2(t) = T_1(t) \times (1 - F_{\text{TE}}(t)) + T_{\text{TE}} \times F_{\text{TE}}(t) \quad (23)$$

In addition, since the pressure relaxes to the ambient pressure during the expansion, we have:

$$P_2(t) = P_1(t) \quad (24)$$

In the case of thermal sparks, the gas temperature is already $T_{\text{TE}}$ after the pulse. Thus, VT relaxation does not occur and $F_{\text{TE}}(t)$ must be set to 0 in Eqn. (21).

## IV.D. Diffusion submodel

At the end of the blast wave expansion and VT relaxation, the plasma kernel is relatively uniform at the temperature $T_{\text{TE}}$ and pressure $P_{\text{out}}$. Furthermore, its radius has increased because of the displacement of the contact surface, as sketched in Figure 13. Figure 11 (top left panel) illustrates this radius change quantitatively for Case 2 (configuration described in Subsection IV.A). The plasma channel radius was 100 µm at $t = 0$ µs and is more than 300 µm at $t = 1$ µs. Figure 17 shows the evolution of the discharge radius $R_{\text{plasma}}(t)$ normalized by its initial value $R_d$. This calculation was done using the same in-house 1D gas dynamics code as in Subsection IV.B (solver of the Euler system of equations). We see that the radius converges to the value $\xi_{\text{expand}}R_d$ with $\xi_{\text{expand}} \approx 3.2$ in this example.

<figure><img src="image_19.png" /><figcaption>Figure 17. Evolution of the discharge radius (contact surface) during the blast wave expansion phase. The initial configuration corresponds to Case 2 ($T_k = 5500$ K, $R_d = 100$ µm).</figcaption></figure>
The value of $\xi_{expand}$ is critical to model radial diffusion properly. Thus, we aim to predict accurately $\xi_{expand}$ and its dependence on the discharge parameters. Considering that:

(a1) the plasma remained enclosed within the contact surface (represented by the shaded orange volume $V = \pi R_d^2 L$ in Figure 13) during the expansion,

(a2) the expansion was a reversible adiabatic process,

We can deduce that the isentropic relation $PV^{\gamma} = \text{constant}$ becomes $P_k \times (\pi R_d^2 L)^{\gamma} = P_{rarefaction} \times (\pi \xi_{expand}^2 R_d^2 L)^{\gamma}$ in which we can then substitute $P_{rarefaction}$ with $P_{out}$ (Subsection IV.B). We also make the assumptions:

(a3) that the total number of particles $N$ inside the plasma is constant during the expansion (i.e., that there is no recombination),

(a4) that the gas follows the ideal gas law.

Thus, we derive $\xi_{expand}$ in closed-form:

$$ \xi_{expand} = \left( \frac{T_k}{T_{out}} \right)^{1/(2\gamma)} \qquad (25) $$

Equation (25) is illustrated in Figure 18 for different values of the adiabatic index. We compared the values of $\xi_{expand}$ obtained from Eqn. (25) with those extracted from the 1D CFD simulations of Subsection IV.B and Figure 11. There were slight differences attributed to difficulties in accurately locating the contact surface, which is thickened by numerical diffusion, as well as the presence of irreversible processes like secondary shocks within the plasma kernel. However, in all instances, the discrepancies did not exceed a few percent, which supports the validity of Eqn. (25).

<figure><img src="image_20.png" /><figcaption>Figure 18. Expansion ratio of the discharges. The colored points are solutions of the 1D Euler system; the grey lines are the power-law fits for different values of $\gamma$: 5/3, 7/5, 8/6, 9/7, 15/13, 26/24 (from bottom to top).</figcaption></figure>
Having obtained simple expressions for $\xi_{\text{expand}}$, we now discuss the radial diffusion submodel itself, its physical assumptions, and the closed-form formulas describing the evolution of P and T in the kernel.

In addition to assumptions (a1)-(a4) made before to compute $\xi_{\text{expand}}$, the submodel relies on the following assumptions:

(a5) The system is 1D and axisymmetric with respect to the inter-electrode axis.

(a6) Radial diffusion begins when the entire plasma kernel has already relaxed to the ambient pressure $P_{\text{out}}$ and behaves as an incompressible gas.

(a7) The fluid is calorically perfect and has a constant thermal diffusivity $\alpha$.

(a8) At the start of the diffusion phase, the gas is at rest.

(a9) The gas is also at TE because VT relaxation is already complete.

(a10) The initial temperature profile follows the step function given by Eqn. (26).

(a11) We consider that the temperature remains constant (equal to $T_{\text{out}}$) for $r \ge \xi_{\text{expand}}R_{\text{d}}$.

$$T(r, t = 0) = \begin{cases} T_{\text{TE}} & \text{if } r < \xi_{\text{expand}}R_{\text{d}} \\ T_{\text{out}} & \text{if } r \ge \xi_{\text{expand}}R_{\text{d}} \end{cases} \qquad (26)$$

(a6), (a8), and (a10) are justified by the velocity and temperature profiles at $t = 1 - 2 \mu s$ shown in Figure 11. In this simulation, the Mach number is below 0.05 after 1 $\mu s$. Note that (a5) implies that we neglect heat losses due to the flux towards the electrodes, so cooling is underestimated, especially for short gap lengths. (a9) is made to develop the reasoning, but is removed at the end of the discussion. In addition, (a7) and (a11) will be removed in a refined submodel. Finally, the sensitivity of the model to (a10) will also be estimated.

Assuming (a4)-(a10), the energy equation of the plasma kernel reduces to the 1D axisymmetric heat equation:

$$\frac{\partial T}{\partial t} = \frac{1}{r} \frac{\partial}{\partial r} \left( \alpha r \frac{\partial T}{\partial r} \right) \quad \text{with} \quad \begin{cases} T(\xi_{\text{expand}}R_{\text{d}}, t) = T_{\text{out}} \quad \forall t \ge 0 \\ \partial T / \partial r \Big|_{r=0} = 0 \quad \forall t \ge 0 \\ T(r, t=0) \quad \text{is the function} \quad (26) \end{cases} \qquad (27)$$

Equation (27) has an analytical solution [61] which reduces to Eqn. (28) at $r = 0$:

$$T(r = 0, t) = T_{\text{out}} + \frac{2}{\xi_{\text{expand}} R_{\text{d}}} (T_{\text{TE}} - T_{\text{out}}) \sum_{n=1}^{+\infty} \frac{1}{\omega_n J_1(\omega_n \xi_{\text{expand}} R_{\text{d}})} e^{-\alpha \omega_n^2 t} \qquad (28)$$

$J_{i \in \mathbb{N}}$ are the Bessel functions of the first kind and $\omega_{n \in \mathbb{N}^*}$ are the roots of Eqn. (29):

$$J_0(\omega_n \xi_{\text{expand}} R_{\text{d}}) = 0 \qquad (29)$$
Table 2 gives the first six roots $\Omega_n$ of $J_0$ (corresponding to the first six values of the product $\omega_n \xi_{\text{expand}} R_d$) [62]. These first six roots are sufficient to calculate the temperature on the scale of microseconds. However, calculating the temperature at the nanosecond scale requires several hundred roots to reduce the error caused by truncating the sum in Eqn. (28).

<table><caption>Table 2. The first six roots of the first kind Bessel function $J_0$.</caption><thead><tr><th>n</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Ω<sub>n</sub></td><td>2.4048</td><td>5.5201</td><td>8.6537</td><td>11.7915</td><td>14.9309</td><td>18.0711</td></tr></tbody></table>

As an example, we consider Case 2, for which the kernel cools from $T_k = 5500$ K to $T_{\text{rarefaction}} = 2850$ K during expansion. In parallel, its radius increased by a factor $\xi_{\text{expand}} \approx 3.2$. Due to the LTE hypothesis in this specific example, $T_{\text{rarefaction}} = T_{\text{TE}}$. Figure 19 illustrates the radial temperature profiles at different instants (a), and the temporal evolution of the kernel temperature (b) computed with Eqn. (28).

<figure><img src="image_21.png" /><figcaption>Figure 19. (a) Radial temperature profiles during the diffusion cooling phase for the configuration of Case 2 after the expansion studied in Section IV.B ($T_k = 5500$ K, $T_{\text{TE}} = 2850$ K). The profiles are computed with formulas from Ref [61]. (b) Temperature at the center of the plasma kernel calculated with Eqn. (28).</figcaption></figure>

In summary, the diffusion submodel takes $T_{\text{TE}}$ as input, and returns the temperature $T_3(t)$ with Eqn. (28) and the pressure with $P_3(t) = P_{\text{out}}$.

As the onset of diffusion cooling may occur before the end of VT relaxation (Subsection IV.C) phases, we propose removing the assumption (a9) and loosely coupling the three phenomena by substituting $T_{\text{TE}}$ in Eqn. (28) with $T_2(t)$, and setting $P_3(t)$ to $P_2(t)$.
Consequently, the coupled expansion, vibrational relaxation, and diffusion submodels reduce to:

$$T_3(t) = T_{\text{out}} + \sum_{n=1}^{+\infty} \frac{2(T_2(t) - T_{\text{out}})}{\Omega_n J_1(\Omega_n)} e^{-\left(\frac{\Omega_n}{\xi_{\text{expand}} R_d}\right)^2 \alpha t} \quad (30)$$

$$P_3(t) \equiv P_2(t) \qquad (31)$$

One of the main limitations of this submodel arises from the key assumption (a7): the constant thermal diffusivity. While it is acceptable within the temperature range 300–2000 K, its validity becomes questionable at higher temperatures. To address these limitations while preserving the analytical nature of the solutions, we divide the time range into several intervals $[t_m; t_{m+1}]$. During each time interval, the thermal diffusivity of the gas is assumed constant, equal to its value at the temperature $T(t_m)$. Thus, we estimate that the temperature of the plasma is:

$$T_3(t) = T_{\text{out}} + \sum_{n=1}^{+\infty} \frac{2(T_2(t) - T_{\text{out}})}{\Omega_n J_1(\Omega_n)} e^{-\left(\frac{\Omega_n}{\xi_{\text{expand}} R_d}\right)^2 \Xi(t)} \quad (32)$$

with

$$\Xi(t) = \begin{cases} \alpha_0 t & \text{if } 0 \le t \le t_1 \\ \alpha_0(t_1 - t_0) + \alpha_1(t - t_1) & \text{if } t_1 \le t \le t_2 \\ \vdots & \vdots \\ \sum_{m=0}^{j-1} \alpha_m (t_{m+1} - t_m) + \alpha_j(t - t_j) & \text{otherwise} \end{cases} \quad (33)$$

where $\alpha_m = \alpha(T_3(t_m))$ is the thermal diffusivity on the time interval $[t_m; t_{m+1}]$ and $j$ is the index of the time interval such that $t_j < t \le t_{j+1}$. We refer to this model as the intermediate diffusion submodel.

Additionally, we can remove the assumption of a Dirichlet boundary condition at a fixed radius (a11). Indeed, it tends to overestimate the cooling effect. Thus, we now assume that it applies at a boundary that moves outwards, with a penetration depth proportional to $\sqrt{\alpha t}$. This approach provides a straightforward method for accounting for the penetration of heat into the surrounding gas. Consequently, Eqn. (32) transforms into Eqn. (34).

$$T_3(t) = T_{\text{out}} + \sum_{n=1}^{+\infty} \frac{2(T_2(t) - T_{\text{out}})}{\Omega_n J_1(\Omega_n)} \Upsilon_n(t) \quad (34)$$
with

$$ \Upsilon_n(t) = \begin{cases} e\left(-\left(\frac{\Omega_n}{\xi_{\text{expand}} R_d}\right)^2 \alpha_0 t\right) & \text{if } 0 \le t \le t_1 \\ e\left(-\left(\frac{\Omega_n}{R_{\text{eq},0}}\right)^2 \alpha_1(t_1-t_0)\right) \times e\left(-\left(\frac{\Omega_n}{R_{\text{eq},1}}\right)^2 \alpha_1(t-t_1)\right) & \text{if } t_1 \le t \le t_2 \\ \vdots & \vdots \\ \prod_{m=0}^{j-1} e\left(-\left(\frac{\Omega_n}{R_{\text{eq},m}}\right)^2 \alpha_m(t_{m+1}-t_m)\right) \times e\left(-\left(\frac{\Omega_n}{R_{\text{eq},j}}\right)^2 \alpha_j(t-t_j)\right) & \text{otherwise} \end{cases} \quad (35) $$

Where $R_{\text{eq},m} = \xi_{\text{expand}} R_d + \sqrt{\Xi(t_m)}$ is the extended radius on the time interval $[t_m; t_{m+1}]$. We refer to this model as the refined diffusion submodel.

Considering the heat equation (28) with temperature-dependent diffusivity $\alpha(T)$, we test the refined model against reference solutions obtained with the 1D finite-volume solver described at the end of Subsection III.A. To illustrate the influence of the initial temperature profile, we conducted three simulations with different initial temperature profiles: a step function, a Gaussian function, and a Bessel function of the first kind $J_0$. The latter two profiles were scaled to match the energy contained in the first profile. The three corresponding temperature evolutions at the center within the configuration of Case 2, displayed in green in Figure 20a, are similar and serve as our baseline.

Figure 20a illustrates the predictions of different versions of the analytical diffusion submodel:

i. The simple submodel based on Eqn. (30) assumes a constant thermal diffusivity fixed at $\alpha(T_{\text{out}})$ (gray squares) or $\alpha(T_{\text{TE}})$ (gray triangles) and no heat penetration into the surrounding gas.

ii. An intermediate submodel based on (32), which includes the temperature variation of thermal diffusivity, but still excludes heat penetration. The predictions are represented by the orange solid line.

iii. The refined submodel based on Eqn. (34). The predictions are displayed in purple.

The time interval 0–200 µs is divided into 512 intervals. At the end of the interpulse, the simple submodel differs from the baseline model by approximately 100%. The deviation is reduced to ~50% with the intermediate submodel. The refined submodel is in good agreement with the baseline solution during the interpulse, with differences that do not exceed 20% compared to the simulation using a Bessel function as an initial profile.

Finally, we verified that the predictions of the refined submodel are not sensitive to the time step. Figure 20b shows the difference between the predictions of $T(r = 0, t = 200 \text{ µs})$ obtained with the refined submodel and the simulation (initialized with a Bessel-like initial profile). Regardless of the method for time subdivision (linear or logarithmic),
the refined submodel agrees within 20% with the simulation when the 0–200 µs time
interval is split into more than 50 segments.

<figure><img src="image_22.png" /><figcaption>Figure 20. (a) Evolution of the temperature at the center of the kernel for the configuration of Case 2 after the expansion studied in Section IV.B ($T_k = 5500$ K, $T_{TE} \approx 2850$ K). The baseline solution obtained with the finite difference method is shown in green with different line styles depending on the initial profile $T(r, t = 0)$. The simple submodel (constant $\alpha$, without heat penetration), intermediate submodel (temperature-dependent $\alpha$, without heat penetration), and refined submodel (temperature-dependent $\alpha$, with heat penetration) are depicted in gray, orange, and purple, respectively. (b) Relative difference in the temperature $T(r = 0, t = 200$ µs) predicted with the refined model versus the number of time divisions.</figcaption></figure>

Another limitation of the diffusion submodel is the coupling with VT relaxation. Indeed, the range of VT timescales is very broad (100 ns–1 ms) and may superpose with the diffusion timescale. Thus, a splitting error may arise if the phenomena occur simultaneously. The formalism chosen in this model addresses this issue. Diffusion depends at any instant on the value of $T_2(t)$. Thus, any increase in $T_2(t)$ due to VT relaxation will increase $T_3(t)$. The rarefaction process occurring significantly earlier than these two processes, this intricate coupling method is particularly useful for diffusion–VT relaxation coupling.

IV.E. Recirculation submodel

The fourth interpulse mechanism is the shock-induced recirculation discussed in Ref. [12]. When the shock wave detaches from the plasma kernel, vorticity is created near the tips of the electrodes. If the induced vorticity is strong enough, the flow starts to rotate, and the hot gases may be expelled radially. Otherwise, vorticity dissipates through friction, and the plasma cools by radial diffusion (see Figure 21). We model this phenomenon using a two-step approach. First, we determine the amount of vorticity induced by the
shock wave. Second, we use this estimate to evaluate the inter-electrode volume cooled by the vortex ring.

<figure><img src="image_23.png" /><figcaption>Figure 21. Main flow dynamics processes during the post-discharge.</figcaption></figure>

IV.E.1. Calculation of Instantaneous Shock-Induced Vorticity (ISIV)

In Ref. [12], we showed that the vorticity associated with the recirculating flow is mainly produced by the “Instantaneous Shock Induced Vorticity” (ISIV) mechanism. Details and assumptions can be found in Ref. [12], but we summarize the main elements for convenience here. A non-uniform shock traveling through a quiescent gas generates post-shock vorticity $\Delta\omega$ on the timescale of a few mean free times [63]. The expression of $\Delta\omega$ reads:

$$ \Delta\omega = \left[ -\vec{n} \times (c_{\text{up}} \vec{\nabla}_t M) \left( \frac{\rho_{\text{down}}}{\rho_{\text{up}}} - \frac{\rho_{\text{up}}}{\rho_{\text{down}}} - 2 \right) \right] \cdot \vec{e}_{\theta} \quad (36) $$
$$ \frac{\rho_{\text{down}}}{\rho_{\text{up}}} = \frac{(\gamma + 1)M^2}{(\gamma - 1)M^2 + 2} \qquad (37) $$

where $M$ is the local shock Mach number, subscripts “up” and “down” stand for the upstream and downstream flows (in the reference frame of the shock), and vectors $\vec{n}$ and $\vec{t}$ are the normal and tangential vectors to the shock front. The non-uniformity of the post-discharge shock wave, which is responsible for generating vorticity, arises from geometric factors. The leading shock of the blast wave takes on an elongated ovoid shape as it detaches from the ovoid plasma kernel. Near the electrodes, the shock front expands spherically, whereas it expands cylindrically in the plane of symmetry $z = 0$ (see Figure 8). Consequently, the shock front decelerates more rapidly near the electrodes than it does near mid-gap. This results in a tangential gradient $\vec{\nabla}_t M$ of the shock Mach number $M$ and to non-zero values of $\Delta\omega$ (Eqn. (36)).

Eqns. (36) and (37) indicate that the total amount of vorticity induced by the ISIV mechanism depends on the shape of the shock front, its strength (as determined by the local Mach number), the adiabatic index $\gamma$, and the ambient sound speed $c_{\text{up}}$. We use the Geometrical Shock Dynamics (GSD) method [64,65] to calculate the temporal evolution of the shape and local Mach number of the shock. Then, we substitute the results into Eqns. (36)–(37) to obtain $\Delta\omega$ across the shock. By ignoring the evolution of the gas downstream, the GSD method significantly reduces the computational time compared to traditional CFD simulations.

Figure 22a illustrates the GSD method. The flow field is divided into stream tubes, which in turn are divided into segments. The inflow surface of a segment corresponds to the shock front at a specified time $t$. Its area is $A(t)$ and the shock Mach number is $M(t)$. The lateral boundaries of the stream tube are perpendicular to the shock surface at $t$ and their length $dn$ is the distance traveled by the shock during $dt$, i.e. $dn = M(t) \times c_{\text{up}} \times dt$. The outflow surface of the segment corresponds to the shock front at time $t + dt$. Its area is $A + dA$, where $dA$ is determined by geometric considerations. The shock Mach number $M(t + dt)$ is equal to $M(t) + dM$, where $dM$ is given by the Chester-Chisnell-Whitham relation [66]:

$$ \frac{M}{M^2 - 1} \lambda_{\text{GSD}}(M, \gamma) \frac{dM}{dn} = -\frac{1}{A} \frac{dA}{dn} \qquad (38) $$

$\lambda_{\text{GSD}}$ is a function of the shock Mach number and the adiabatic index $\gamma$. The GSD method operates as a time-marching technique, advancing the shock front from one time increment to the next.

The method is valid under several key assumptions.

(a1) The flow can be described by the Euler system of equations during the early expansion of the shock. This was verified in Subsection IV.B.
(a2) The upstream flow velocity is perpendicular to the shock front. In our scenario, this condition is satisfied since the shock speed is significantly greater than the velocity of the ambient gas, which we can consider to be at rest (at least during the first microsecond, when the vorticity is created).

(a3) The gas is calorically perfect with constant molar mass and adiabatic index. This assumption may break under the most stringent conditions where $T_k/T_{out} \approx 20$ because the gas gets internally excited or dissociates at high temperatures $T_k$. However, Ref. [67] showed that the thermochemical effects remain mild when the Mach number of the leading shock does not exceed 10, hence our strategy remains valid.

(a4) The multiple wave reflections due to the non-planar geometry have a negligible impact on the leading shock, a condition that was verified in Section 8.1 of Ref. [66].

We add assumptions to previous ones, specifically tailored to our use of GSD.

(a5) The geometry is 2D and cylindrical.

(a6) The kernel has an elliptical shape. We show below that this assumption has a negligible impact on the model.
<figure><img src="image_24.png" /><figcaption>Figure 22. (a) Schematic representation of the GSD method to predict the blast wave leading shock propagation. The fluid domain is divided into 150 stream tubes with varying cross-sectional areas. The violet contour depicts the initial shock front ($t_1 = 0^+$ s) whose local Mach number $M_s$ is uniform along the entire shock profile. The blue contour shows the evolved shock front ($t_2 > 0$ s) whose local Mach number $M$ is non-uniform. (b) The right panel presents a magnified view of the left panel at $t_1 = 0^+$ s and outlines the analogy with the shock tube problem.</figcaption></figure>

The GSD method is initialized as follows. First, the initial shock front is assumed to have the same ovoid shape as those defined in Ref. [12] and used in the CFD calculations of Section III. Second, the local Mach number $M_s(t = 0^+)$ of the shock is assumed to be uniform along the entire shock profile. During the first instants of blast wave relaxation, an analogy can be drawn between the stream tubes and a fictitious shock tube with 1D planar geometry. This is illustrated in Figure 22b. The initial Mach number $M_s(t = 0^+)$ is then retrieved by solving the usual shock tube relation (obtained from equations 7.13 and 7.94 of Ref. [36]):

$$ \frac{P_{\text{out}}}{P_{\text{k}}} \left( \frac{2\gamma_{\text{out}}M_{\text{s}}^2}{\gamma_{\text{out}} + 1} - \frac{\gamma_{\text{out}} - 1}{\gamma_{\text{out}} + 1} \right) = \left( 1 - \frac{\gamma_{\text{k}} - 1}{\gamma_{\text{out}} + 1} (M_{\text{s}} - M_{\text{s}}^{-1}) \frac{c_{\text{out}}}{c_{\text{k}}} \right)^{2\gamma_{\text{k}}/(\gamma_{\text{k}} - 1)} \quad (39) $$

where $c_k = \sqrt{\gamma_k(R/M_k)T_k}$ and $c_{\text{out}} = \sqrt{\gamma_{\text{out}}(R/M_{\text{out}})T_{\text{out}}}$ correspond to the speed of sound in the plasma kernel and the ambient gas. $P_k$ is equal to $T_k \times P_{\text{out}}/T_{\text{out}} \times$
$M_{\text{out}}/M_{\text{k}}$ because the pulse heating was assumed to be isochoric in Section II. As for the other submodels, we consider that the adiabatic index and molar mass are the same in the plasma and in the ambient gases.<fnref n="3" /> Equation (39) simplifies into:

$$ \frac{T_{\text{out}}}{T_{\text{k}}} \left( \frac{2\gamma M_{\text{s}}^2}{\gamma + 1} - \frac{\gamma - 1}{\gamma + 1} \right) = \left( 1 - \frac{\gamma - 1}{\gamma + 1} (M_{\text{s}} - M_{\text{s}}^{-1}) \sqrt{\frac{T_{\text{out}}}{T_{\text{k}}}} \right)^{2\gamma/(\gamma-1)} \quad (40) $$

Then, the GSD equations are discretized in space using a fourth-order centered finite difference scheme and integrated in time using an eighth-order Runge-Kutta method [68]. Based on convergence studies, the initial spatial segmentation of the shock is reduced to $dx_{\min} = 0.5$ µm (near the tip of the ellipsoid), and the time step is set to $dt = 0.2$ ns. The simulations are continued up to ~1 µs until the shock is too weak to contribute to the increase in vorticity by the ISIV mechanism. For example, we initialize a calculation with the post-pulse conditions of Case 2, introduced in Subsection IV.B. Figure 23 shows the shock front propagation at different instants. Then, we inject the shock profile and its local Mach number at each timestep into Eqns. (36) and (37) to compute the jump of vorticity created by the shock at each point of the domain. Figure 24 shows a two-dimensional map of the instantaneous shock-induced vorticity $\Delta\omega(r, z)$ over a quarter of the domain.

<figure><img src="image_25.png" /><figcaption>Figure 23. Propagation of the shock front computed with the Geometrical Shock Dynamic (GSD) method. The calculation was initialized with the conditions of Case 2. The horizontal axis represents the inter-electrode axis, and the vertical axis is the radial direction.</figcaption></figure>

<footnote n="3">³ The influence of $\gamma_{\text{k}}$, $\gamma_{\text{out}}$, $M_{\text{k}}$, $M_{\text{out}}$ on $M_{\text{s}}$ has been compared to the one of $T_{\text{k}}/T_{\text{out}}$ in Appendix B of Ref. [12]. We showed that the temperature ratio is significantly more influent.</footnote>
Finally, we integrate $\Delta\omega(r, z)$ over the entire domain to determine the total integrated vorticity $\Gamma$ deposited by the leading shock (ISIV mechanism), see Eqn. (41). $\Gamma$ represents the strength of the vortices appearing in the post-discharge recirculation regime.

$$ \Gamma \equiv \int_{r,z} \Delta\omega(r,z) dzdr \quad [\text{m}^2/\text{s}] \quad (41) $$

<figure><img src="image_26.png" /><figcaption>Figure 24. Instantaneous Shock Induced Vorticity (ISIV) deposited by the leading shock of the blast wave. The vorticity is calculated by applying Eqn. (36) to the shock front predicted by the Geometrical Shock Dynamics (GSD) method. The simulation was initialized with the conditions of Case 2.</figcaption></figure>

$\Gamma$ depends on the discharge parameters. To determine correlations, we evaluate $|\Gamma|$ for thousands of discharge configurations, varying the gap length $L$ from 500 µm to 6 mm, the ambient temperatures $T_{\text{out}}$ from 300 K to 1500 K, the kernel radius $R_d$ from 50 µm to 300 µm, the heating ratio $T_k/T_{\text{out}}$ from 1.05 to 20, and the gas properties $M$ and $\gamma$ among those of Table 1. All simulations are initialized with an ovoid shape whose major and minor axes are $L$ and $2R_d$, respectively. The resulting $\Gamma$ are reported in Figure 25. We observe that $|\Gamma|$ correlates to first order with $T_k/T_{\text{out}}$, and to second order with $\gamma$. The sensitivity to $L$, $R_d$, $M$ and $T_{\text{out}}$ is very low. The results can be fit with Eqn. (42):

$$ |\Gamma| = \exp \left[ \beta_3 \left( \frac{T_k}{T_{\text{out}}} - 1 \right)^{\beta_4} + \beta_5 \left( \frac{T_k}{T_{\text{out}}} - 1 \right)^{2\beta_4} + \beta_6 \right] \\ \text{with} \quad \begin{cases} \beta_3 = 88.553 + 85.929\sqrt{\gamma - 1} \\ \beta_4 = -0.020 - 0.014\sqrt{\gamma - 1} \\ \beta_5 = -100 - 22.388\sqrt{\gamma - 1} \\ \beta_6 = 4.143 - 65.078\sqrt{\gamma - 1} \end{cases} \quad (42) $$

In Appendix A, we show that the sensitivity of $\Gamma$ to the initial kernel shape is mild (~30%) compared with the sensitivity to $T_k/T_{\text{out}}$. Thus, Eqn. (42) holds to a good approximation for any post-discharge, regardless of the shape of the initial plasma kernel (assumption (a6)). The present evidence for the strong dependence of $|\Gamma|$ on $\Pi^* \equiv T_k/T_{\text{out}}$ is consistent with the preliminary demonstration provided in Section 4.3 of Ref. [12].
<figure><img src="image_27.png" /><figcaption>Figure 25. Total integrated vorticity deposited in the ambient gas by the ISIV mechanism for different discharge configurations (L, R<sub>d</sub>, T<sub>out</sub>, T<sub>k</sub>/T<sub>out</sub>). The data are fit with a second-order polynomial. The grey lines are the power-law fits for different values of γ: 5/3, 7/5, 8/6, 9/7, 15/13, 26/24 (from bottom to top).</figcaption></figure>

IV.E.2. Effect of the ISIV mechanism on the temperature profile

The total integrated vorticity |Γ|, defined in Eqn. (42), is the key quantity of the recirculation model. We now present simple arguments for estimating the evolution of the mean temperature T<sub>4</sub>(t) and pressure P<sub>4</sub>(t) along the inter-electrode axis, considering the three main hydrodynamics phenomena (expansion, diffusion, and recirculation) as well as any potential vibrational relaxation.

The generation of vorticity by the ISIV mechanism occurs in a narrow region near the electrode tip, as shown in Figure 24 and Ref. [12]. A vortex ring forms and then self-advects towards the inter-electrode plane [16,69]. The vortex ring brings cold gases into the inter-electrode region and contributes to its cooling. We assume that the temperature of the advected cold gas is T<sub>out</sub>. Thus, the temperature of the inter-electrode region is a weighted average (see Eqn. (43)) between the discharge temperature T<sub>3</sub>(t) (gases already cooled by expansion and diffusion) and the ambient temperature T<sub>out</sub>:

$$T_4(t) = T_{\text{out}} \times F_{\text{advection}}(t) + T_3(t) \times (1 - F_{\text{advection}}(t)) \quad (43)$$

The weighting factor F<sub>advection</sub>(t) is proportional to the volume of gas introduced into the inter-electrode region by the vortex. Considering the situation sketched in Figure 26 and the characteristic time of self-advection τ<sub>vortex</sub>, we have:

$$\tau_{\text{vortex}} \approx L^2 / |\Gamma| \quad (44)$$
<figure><img src="image_28.png" /><figcaption>Figure 26. Schematic of the vortex self-advection of cold gas into the inter-electrode gap.</figcaption></figure>

This advection process starts when the shock has deposited most of the shock-induced vorticity, i.e., when it has moved by a radial distance on the order of $R_d$ (see Figure 24). We suggest the weighting factor $F_{\text{advection}}(t)$, a simple differentiable function that smoothly connects the values 0 and 1:

$$F_{\text{advection}}(t) = \begin{cases} 0 & \text{if } t < \tau_{\text{shock}} \\ \frac{1}{2} \left(1 - \cos\left(\pi \frac{t - \tau_{\text{shock}}}{\tau_{\text{vortex}}}\right)\right) & \text{if } \tau_{\text{shock}} \le t \le \tau_{\text{shock}} + \tau_{\text{vortex}} \\ 1 & \text{if } t > \tau_{\text{shock}} + \tau_{\text{vortex}} \end{cases} \quad (45)$$

$\tau_{\text{shock}}$ is the duration needed by the shock to travel a distance $R_d$:

$$\tau_{\text{shock}} = \frac{R_d}{M_s c_{\text{out}}} \quad (46)$$

We calculate a large set of $M_s(t = 0^+)$ by solving Eqn. (39) for different values of $\gamma$ and $M$. Figure 27 displays the results and the Eqn. (47) obtained by fitting a power law.

$$M_s = \beta_7 + (1 - \beta_7) \left(\frac{T_k}{T_{\text{out}}}\right)^{\beta_8} \quad \text{with} \quad \begin{cases} \beta_7 = 0.452 + 0.084\gamma \\ \beta_8 = 0.597 - 0.020\gamma \end{cases} \quad (47)$$
<figure><img src="image_29.png" /><figcaption>Figure 27. Initial Mach number of the shock wave leaving the kernel. The colored points are solutions of Eqn. (39); the grey lines are the power-law fits for different values of $\gamma$: 5/3, 7/5, 8/6, 9/7, 15/13, 26/24 (from bottom to top).</figcaption></figure>

Finally, we assume that recirculation cooling is isobaric at $P_{out}$ because $\tau_{rarefaction} \ll \tau_{vortex}$. This translates into:

$$P_4(t) \equiv P_3(t) \qquad (48)$$

This strategy may overestimate the impact of the ISIV mechanism in the recirculating regimes because of damping mechanisms, such as the viscous effect or secondary sources of vorticity.

## IV.F. Summary

In Subsections IV.B, IV.C, IV.D, and IV.E, we presented an engineering model to predict the mean thermodynamic state along the inter-electrode axis by considering the three main flow dynamics phenomena (blast wave expansion, radial diffusion, and shock-induced recirculation) and the VT energy exchange (of molecular gases). This model couples four intricate submodels, each of which contributes to the analytical expressions $T_4(t)$ (Eqn. (43)) and $P_4(t)$ (Eqn. (48)).

Furthermore, in a recirculating flow regime, the thermodynamic state of the advected gases remains unaffected by the recirculation movement. Therefore, the expressions $T_3(t)$ (Eqn. (34)) and $P_3(t)$ (Eqn. (31)), which account for expansion, VT energy transfer, and radial diffusion but exclude the contribution of recirculation, are well-suited to estimate the thermodynamic state inside the expelled torus of hot matter. In such flows, the
chemistry of interest may be located inside the torus rather than in the inter-electrode
region.

Figure 28 summarizes our engineering model for the post-discharge fluid dynamics presented in Section III. For the sake of completeness, we have also included the pulse model introduced in Section II.
<figure><img src="image_30.png" /><figcaption>Figure 28. Engineering model for nanosecond pulse and post-discharge.</figcaption></figure>
V. Interpulse model: validation of the flow dynamics submodels

This section tests and validates the engineering interpulse model using a three-step
approach. First, we compare its predictions of temperature $T_4(t)$ and pressure $P_4(t)$ with
the results of CFD simulations in the two typical air post-discharge configurations
detailed in Subsection III.B. Second, we analyze the temperature of the inter-electrode
region at the end of the interpulse ($t = 1/f$) for a wide range of air discharge
configurations, including various gap lengths $L$, and plasma kernel temperatures $T_k$. We
compare the predictions of the engineering model to the results obtained from CFD
simulations. Third, we challenge the reduced-order model by examining extreme cases of
air discharges ($T_k \approx 20,000$ K), where some of the initial assumptions may not hold. The
detailed comparison with CFD simulations will help assess the limits of the engineering
model.

For the validation cases below, we initialize the gas kernel at LTE for a given pressure
and temperature ($P_k, T_k$). The interpulse simulation involves the different post-discharge
flow dynamics phenomena: blast wave expansion, radial diffusion, and shock-induced
recirculation. This section does not test the energy relaxation submodel because the CFD
code does not include a multi-temperature thermodynamic model.

V.A. Comparison of temperature and pressure profiles versus
time

In this Subsection, we demonstrate that the reduced-order flow dynamics model correctly
captures the temporal evolution of pressure and temperature inside the filament over a
wide range of timescales. We consider the two cases presented in Subsection III.B because
they are typical of the recirculation and diffusive regimes encountered in post-discharges.
We compare $T_4(t)$ and $P_4(t)$, defined in Eqns. (43) and (48), with the corresponding
predictions of two-dimensional CFD simulations ($\bar{T}(t)$ and $\bar{P}(t)$, defined in Eqns. (10)–
(11)). The results are reported in Figure 29 and Figure 30. In both figures, the upper
panel (1) corresponds to a post-discharge in a diffusion regime (Case 1), whereas the
lower panel (2) corresponds to a recirculation regime (Case 2).

Regarding temperature, Figure 29 shows that the engineering model agrees with the
baseline CFD method within ~20%:

* In Cases 1 and 2, the fall time and the intermediate temperature plateau $T_{\text{rarefaction}}$ after the blast wave expansion (phase ①) are predicted within ~20%. The differences are mainly explained by numerical diffusion. However, these discrepancies are mild in comparison to the temperature fall from $T_k$ to $T_{\text{rarefaction}}$. The onset of temperature and pressure decreases occurs twice as early as in the CFD simulations. This is due to the choice of kernel shape used in the simulations, which is elliptic and not purely cylindrical. Thus, the rarefaction waves starting
near the tips of the electrodes cross a distance shorter than $R_d$ and make the mean pressure and temperature fall earlier. Nevertheless, the model catches the main trends and the characteristic time well.

* In Case 1, the engineering model correctly predicts the absence of recirculation (Figure 29a) and the dominant contribution of diffusion cooling (phase ③). The engineering model agrees with the CFD simulations within ~10%. Both cooling rates of diffusion are similar, but start from a higher temperature with the model. The difference mainly stems from the deviation introduced by the rarefaction submodel.

* In Case 2 (Figure 29b), the temperature drop caused by recirculation (phase ④) is relatively well captured, both in terms of onset and fall times. The observed delay of 4 µs between both recirculation start time predictions remains small compared to the fall time required for the vortex to cool the kernel from $T_{\text{rarefaction}}$ to $T_{\text{out}}$. This delay is likely due to two phenomena that were neglected in the engineering model: (i) small vortices induced by the shock require an induction time to coalesce into a large self-advecting vortex, (ii) the vortex slows down due to counter-rotating vorticity sources (such as baroclinic torque, expansion) or viscous dissipation.

Regarding pressure evolution, the engineering and CFD models are compared in Figure 30. Although the reduced model neglects the overpressure caused by the secondary shock (marked by * in the figure), both models agree regarding the decay time and the steady-state pressure value. The onset of rarefaction has the same issue as what we remarked in the temperature study: the CFD post-processed data fall twice as early as with the engineering model, due to the non-cylindrical initial shape of the kernel in the CFD simulation.

In summary, the engineering model developed in Section IV can effectively predict the temperature and pressure evolution during the interpulse, as calculated by the baseline CFD method. Although this model relies on approximations, the resulting pressure and temperature differences are small compared to the variations observed throughout a complete interpulse. Its level of accuracy is adequate for conducting studies on chemical kinetics or preliminary reactor designs.
<figure><img src="image_31.png" /></figure>

Figure 29. Mean temperature of the plasma kernel along the pin-to-pin axis returned by the CFD method (solid line) and by the engineering model (dashed line). (a) Case 1: $L = 1$ mm, $T_k = 500$ K. (b) Case 2: $L = 1$ mm, $T_k = 5500$ K. The dagger (†) marks the blast wave sub-atmospheric depression. The asterisk (*) marks the temperature peak due to the secondary shock (see text).

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

Figure 30. Mean pressure of the plasma kernel along the pin-to-pin axis returned by the CFD method (solid line) and by the engineering model (dashed line). (a) Case 1: $L = 1$ mm, $T_k = 500$ K. (b) Case 2: $L = 1$ mm, $T_k = 5500$ K. The dagger (†) marks the blast wave sub-atmospheric depression. The asterisk (*) marks the temperature peak due to the secondary shock (see text).

## V.B. Comparison of temperature at the end of the interpulse

We now analyze the temperature of the inter-electrode region at the end of the interpulse ($t = 1/f$) for a wide range of NRP discharge configurations, including various gap lengths $L$ and plasma kernel temperatures $T_k$. We compare the predictions of the engineering
model with those obtained from CFD simulations. This strategy extends the one presented in Section V.A, which was only carried out in two discharge configurations.

First, we select the CFD database of Ref. [12], which contains the temperature field of tens of CFD simulations for a variety of L (1–5 mm) and T_k (500–6000 K). The ambient temperature is 300 K in all simulations. For each element in the CFD dataset, we compute the average final inter-electrode temperature $\bar{T}_{end}$ (see Eqn. (49)) with a frequency of f = 20 kHz:

$$ \bar{T}_{end} \equiv \bar{T}(t = 1/f) \qquad (49) $$

Second, we calculate $T_1(t = 1/f)$, $T_3(t = 1/f)$, and $T_4(t = 1/f)$ using the engineering model (Eqns. (17), (30), and (43)). This splitting allows to analyze the separate contributions of blast wave expansion, radial diffusion, and shock-induced recirculation:

* $T_1$ accounts for the effect of expansion only,

* $T_3$ includes the effects of expansion and diffusion,

* $T_4$ includes the effects of expansion, diffusion, and recirculation.

Finally, we compare $\bar{T}_{end}$ (obtained from CFD simulations) to $T_{i \in \{1,3,4\}}(t = 1/f)$. This comparison is shown in Figure 31. The outputs ($T_1(t = 1/f)$, $T_3(t = 1/f)$, $T_4(t = 1/f)$, and $\bar{T}_{end}$) are plotted on a color map for each combination of inputs (L, T_k) given by the x and y entries.

The upper left panel reveals that the blast wave expansion approximately halves the temperature in the conditions studied, regardless of the gap length. The upper right panel shows that radial diffusion, whose effect is added to that of the expansion, contributes to the temperature decrease below 1500 K, regardless of the gap. The lower left panel displays the combined effects of the blast wave expansion, diffusion, and recirculation. By comparison with the upper right panel, where recirculation is ignored, we see that circulation has a strong cooling effect in the region labeled A. In this region, the shock-induced recirculation completely replaces the hot plasma kernel at $T_k$ by fresh gas at temperature $T_{out}$. Similarly, we identify region D, dominated by expansion and diffusion cooling. Region B includes the conditions where the contribution of diffusion is small compared to the blast wave expansion and the recirculation. In region C, the gas is partially cooled by expansion and diffusion. Note that the extent of these regions depends on the pulse frequency f: increasing the interpulse period T = 1/f gives more time for recirculation and diffusion to act. Further details about regions A–D are given in Section III of Ref. [12].

As a final validation test, we compare the lower left panel ($T_4$) obtained with the engineering model and the lower right panel ($\bar{T}_{end}$) obtained from CFD simulations. A good agreement is obtained, both in terms of regime boundaries and quantitative
temperatures. The temperature discrepancies do not exceed 20%, which is small compared to the difference between the initial ($T_k$) and final ($\bar{T}_{end}$) temperatures. The deviation is mainly caused by the assumption of 1D isentropic expansion (Subsection IV.B) and the approximation of the ovoid kernel by a straight cylinder in which the edge effects of diffusion cooling are neglected (Subsection IV.D).

In summary, the agreement with CFD simulations, as shown in Figures 29–30 and the lower panels of Figure 31, demonstrates the predictive capabilities of the engineering model. The three major flow dynamics mechanisms and their dependence on the discharge parameters are well captured by the model. The model offers two primary advantages: it has a very short computation time (a couple of seconds) and is straightforward to use. In contrast, a CFD simulation, such as the one presented in Section III, demands significant expertise and requires around 1,500 CPU hours on a high-performance computing (HPC) cluster.

<figure><img src="image_33.png" /><figcaption>Figure 31. Average temperature of the inter-electrode gas at the end of the interpulse $t = 50~\mu s$ ($f = 20~\text{kHz}$) for different gap lengths $L$ and kernel temperatures $T_k$. The temperatures in the upper left, upper right, and lower left panels were obtained with Eqns. (17), (30), and (43) of the engineering model. They include the cumulative effects of the blast wave expansion ($T_1$), diffusion ($T_3$), and recirculation ($T_4$). The lower right panel shows $\bar{T}_{end} = \bar{T}(t = 1/f)$ obtained by applying Eqn. (10) to the CFD database of Ref. [12].</figcaption></figure>

## V.C. Limits of the engineering model

In Subsections V.A and V.B, we validated the post-discharge engineering model within the parametric domain $T_k \le 6000~\text{K}$ and $P_k \le 20~\text{bar}$, where $T_k$ and $P_k$ are the temperature and pressure of the plasma kernel at the start of the post-discharge. These parameters are typical of non-thermal sparks. The current Subsection seeks to test the model under the interpulse conditions following thermal sparks, characterized by a
temperature $T_k$ and pressure $P_k$ that can reach 30,000 K and 150 bar, respectively. This test will help quantify the errors that may arise when using the model beyond its applicability limits.

We conduct CFD simulations that take high-temperature effects into account, compare these results with the predictions of the engineering model, and analyze the differences in light of the simplifications made in the engineering model (calorically perfect gas, constant molar masses, constant adiabatic indices).

We consider three CFD models:

i. The first model is similar to the CFD model presented in Section III. Chemistry is assumed to be infinitely fast: the thermodynamic and transport properties are those at LTE [46]. As suggested in Ref. [15], the effect of radiation is neglected. The kernel radius $R_d$ is set to 100 µm, the ambient conditions $T_{\text{out}}$ and $P_{\text{out}}$ are set to 300 K and 1 atm, the gap length is 1 mm. Unlike in Section III, we calculate the initial kernel pressures $P_k$ using Eqn. (3) instead of Eqn. (4) because dissociation or ionization are significant.

ii. The second model (baseline model) has the same assumptions as the first one, except that it includes the non-LTE effects. Chemistry is treated with Park's finite rate model [70]. Ordinary diffusion is accounted for with the multicomponent Maxwell-Stefan diffusion model [71]. The diffusion coefficients of the neutral species are taken from Ref. [72]. The diffusion coefficients of the ions are computed with the ambipolar diffusion approximation, and no charge separation is allowed [73].

iii. The third model uses the same assumptions as the LTE CFD model, except that chemistry is frozen. The heat capacities and the transport properties are fixed at their values at $t = 0$ s [46].

The three CFD models are initialized with air at LTE.

We focus on two typical post-discharges:

* In Case 2,<fnref n="4" /> the plasma kernel is initialized at $T_k = 5500$ K (non-thermal spark),

* In Case 3, the plasma kernel is initialized at $T_k = 20,000$ K (thermal spark).

In both cases, we compute the average inter-electrode pressure and temperature (Eqns. (10)–(11)) using the three CFD models (i–iii) and the engineering model (Eqns. (43) and (48)). Figure 32 compares the CFD models and the engineering model for Cases 2 and 3.

<footnote n="4">⁴ Case numbering continues from where it left off in Subsection V.A, with Case 2.</footnote>
In the low-heating Case 2 ($T_k$ = 5500 K, Figures 32a and 32b), the temperature and pressure predictions of the CFD models overlap. For readability reasons, they are represented with a unique solid black line. The chemistry has a negligible effect on the temperature and pressure evolution after a non-thermal spark. This conclusion is consistent with the findings of Ref. [12]. In addition, the engineering model (dashed red curve) agrees with the CFD models (solid black line), exhibiting the same errors as discussed in Subsection V.A.

In the high-heating Case 3 ($T_k$ = 20,000 K, Figures 32c and 32d), the pressure predictions of the three CFD models overlap. For readability reasons, they are also represented with a unique solid black line. However, the chemistry now significantly influences the temperature evolution. The frozen chemistry model predicts a cooling down to 5000 K ~50 times faster than the equilibrium chemistry model. This is attributed to the large thermal conductivity induced by the high mole fractions of non-recombined electrons. The frozen and LTE models are two limiting cases whose validity can be assessed by comparison with the baseline finite rate chemistry model. The close agreement between the temperature evolutions of the finite rate and LTE models (Figure 32c) proves that the gas is likely in LTE. This also leads to a preference for the LTE model, which is much less computationally intensive than the finite rate model.

The three CFD models predict an abrupt temperature fall at $t \approx 10$ µs, marked by the symbol ⌖ in Figure 32c. Owing to its stiffness, the transition is likely due to recirculation rather than diffusion. The onset of cooling is not influenced by the chemistry, which is consistent with the generation of vorticity by the ISIV mechanism [12]: the vorticity is induced behind the leading shock in a gas whose temperature and pressure are slightly higher than in the surrounding environment (owing to the passage of the shock) but much lower than in the plasma kernel. In this “gentle” post-shock environment, the chemical effects are mild.

The temperature prediction of the engineering model is shown in red in Figure 32c. The model overpredicts the temperature fall during the blast wave expansion (marker ⌖). It is impossible to evaluate its ability to reproduce diffusion cooling, as recirculation cooling occurs beforehand. Furthermore, the engineering model predicts the onset of recirculation (marker ⌖) ~10 µs earlier than the baseline CFD model. It overestimates the intensity of the recirculation vortex and predicts $\tau_{vortex} \approx 1$ µs instead of 10 µs. On the other hand, the pressure evolution calculated with the model is in good agreement with the CFD simulations.

In conclusion, the engineering model is valid for conditions $T_k \le 6000$ K and $P_k \le 20$ bar, typical of non-thermal sparks. In situations where $T_k \ge 6000$ K and $P_k \ge 20$ bar, such as thermal sparks, the engineering model suffers from shortcomings regarding the extent of the temperature fall due to rarefaction () and the onset of the recirculation phase (et). This second shortcoming is addressed in Section VI.
<figure><img src="image_34.png" /><figcaption>Figure 32. Mean temperature (a, c) and pressure (b, d) of the plasma kernel along the pin-to-pin axis obtained with three CFD models (solid lines) and the engineering model (dashed line). The three CFD models differ by their chemistry models: LTE, finite rate, and frozen chemistry. In panels (a), (b), and (d), the predictions of the CFD models overlap and are represented by a single solid dark line. Panels (a) and (b) correspond to non-thermal spark conditions (Case 2). Panels (c) and (d) correspond to thermal spark conditions (Case 3). Markers ● and ⚫ refer to the blast wave expansion and recirculation phases.</figcaption></figure>

## VI. Interpulse model: prospects for refinement

We now suggest a few modifications to improve the engineering submodels in situations not addressed in Section IV. We focus on the influence of non-quiescent surrounding environments (Subsection VI.A) and on the effects of high plasma/ambient temperature ratios $T_k/T_{out}$ on the recirculation model (Subsection IV.B). A refined model is proposed and validated.

### VI.A. External bulk flow

In most practical cases, discharges are not produced in quiescent conditions but in flowing mixtures. Fresh gases are brought into the inter-electrode region by the recirculation cell and/or by the bulk flow (velocity $u_{flow}$), to which we associate a residence time $\tau_{flow}$:

$$ \tau_{flow} = u_{flow}/R_{d} \qquad (50) $$

In practice, the plasma reactor will be controlled by the phenomenon having the shortest residence time ($\tau_{vortex}$ or $\tau_{flow}$). In high-speed flow systems, such as power production
gas turbines [6], CO₂ conversion reactors [74], or scramjet combustors [75], the convection of fresh gas by the bulk flow (at speed $u_{flow}$) into the inter-electrode region may exceed that caused by the shock-induced recirculation vortex. In that case, we use $F_{advection}^{refined}(t)$ instead of $F_{advection}(t)$ inside Eqn. (43) of the recirculation submodel (Subsection IV.E):

$$F_{advection}^{refined}(t) = \min \left( F_{advection}(t), \frac{t}{\tau_{flow}} \right) \qquad (51)$$

where $F_{advection}$ is given by Eqn. (45).

<figure><img src="image_35.png" /><figcaption>Figure 33. Plasma-assisted ignition mechanism. Adapted from [76].</figcaption></figure>

## VI.B. Saturation of the shock-induced vorticity (air discharges)

We now address the shortcoming of the recirculation submodel underscored in Subsection V.C: in a situation typical of thermal spark discharges, where the initial plasma-ambient temperature ratio $T_k/T_{out}$ is approximately 70, the predicted onset of recirculation (at $t = \tau_{vortex}$) by the engineering model was ~10 times earlier than the predictions of CFD models. This suggests that Eqn. (42) overpredicts the integrated vorticity $\Gamma$ by an order of magnitude when $T_k/T_{out} \approx 70$, i.e. well outside the range $T_k/T_{out} = 1-20$ used to fit the correlation (Eqn. (42)). We propose a simple modification to this equation that extends its predictive capabilities to the range $T_k/T_{out} = 20-150$.
We consider four post-discharge conditions:<fnref n="5" />

* Case 3 corresponds to that of Subsection V.C, where $T_k = 20,000$ K and $T_{out} = 300$ K ($T_k/T_{out} = 67$).

* Cases 4, 5, and 6 are identical to Case 3, except that $T_k$ is set to 10,000 K, 35,000 K, and 45,000 K ($T_k/T_{out} = 33, 117$, and $150$).

These cases are simulated with the baseline CFD finite-rate model (listed as (ii) in Subsection V.C). We analyze the temperature evolution, identify the characteristic time $t = \tau_{vortex}^{CFD}$ for recirculation cooling, and convert $\tau_{vortex}$ into an integrated vorticity $|\Gamma|^{CFD}$ using Eqn. (44). For these four cases, $|\Gamma|^{CFD}$ was about $10^{-1}$ m².s⁻¹. This value is similar to the one obtained from Case 2, where $T_k/T_{out} \approx 18$. There is a saturation of the shock-induced vorticity for $T_k/T_{out} \ge 20$. Thus, we suggest replacing the correlation (42) by the correlation (52) to align with the predictions of the CFD simulations.

$$|\Gamma| = \min \left\{ \exp \left[ \beta_3 \left( \frac{T_k}{T_{out}} - 1 \right)^{\beta_4} + \beta_5 \left( \frac{T_k}{T_{out}} - 1 \right)^{2\beta_4} + \beta_6 \right], 10^{-1} \right\} \quad (52)$$

The physical origin of this saturation remains an open question. We suggest that this could be due to other vorticity sources competing with the shock-induced vorticity and slowing down recirculation. One possible source is the baroclinic torque distributed behind the shock, as discussed by Wang et al. in Section 4.2 of their study on laser-induced optical breakdown [77]. However, this would need further investigation.

We took into account the modification and tested this refined model to simulate the high-heating Case 3 discussed in Subsection V.C. The results are shown in Figure 34. As expected, the cooling due to recirculation is consistent between the refined model and the CFD validation cases. However, the limitation of the interpulse model related to the temperature fall during the rarefaction process remains visible in Case 3.

<footnote n="5">⁵ Case numbering continues from where it left off in Subsection V.C, with Case 3.</footnote>
<figure><img src="image_36.png" /><figcaption>Figure 34. Mean temperature (a, c) and pressure (b, d) of the plasma kernel along the pin-to-pin axis obtained with three CFD models (solid lines) and the engineering refined model (dashed line). The three CFD models differ by their chemistry models: LTE, finite rate, and frozen chemistry. In panels (a), (b), and (d), the predictions of the CFD models overlap and are represented by a single solid dark line. Panels (a) and (b) correspond to non-thermal spark conditions: $L = 1$ mm, $T_k = 5500$ K (Case 2). Panels (c) and (d) correspond to thermal spark conditions: $L = 1$ mm, $T_k = 20,000$ K (Case 3).</figcaption></figure>

VII. Simple cooling criteria for practical use by engineers

In this Section, we use the engineering model developed in this article to optimize nanosecond discharges for plasma processing applications.

* In some cases, the process efficiency may be increased by promoting the transport of plasma-created radicals towards the environment. A quick transport maximizes the volume of ambient gases impacted by the plasma-related chemistry. The post-pulse recirculation plays a key role in this transport.

* In other cases, the process efficiency may depend on the ability to keep the inter-electrode zone cold and filled with fresh gas before the next pulse occurs. This approach ensures that the energy from the next pulse is not wasted on already processed gas. Additionally, in the context of Plasma-Assisted Combustion (PAC), efficient cooling may reduce the production of nitrogen oxides by thermal processes, such as the Zeldovich reactions. To optimize carbon dioxide dissociation, interpulse cooling may help reduce the recombination of $CO_2$ during the interpulse.
Thus, we split this Section into two parts. First, we discuss vorticity optimization to maximize recirculation in the interpulse flow regime. Second, we study how to optimize cooling, and we build criteria that may help engineers or experimentalists design nanosecond discharges.

## VII.A. Optimization of vorticity creation and flow regime

In Ref. [12], we demonstrated that the vorticity creation mainly depends on the non-dimensional number $\Pi^*$ defined as:

$$ \Pi^* \equiv \frac{T_k}{T_{\text{out}}} \qquad (53) $$

The $\Pi^*$ number quantifies the heating of the kernel during the pulse, compared to the ambient temperature. We showed, through order-of-magnitude reasoning, that $\Pi^*$ was likely to be the main scaling factor for the vorticity generated by the ISIV mechanism. The shock-induced recirculation submodel presented in Subsection IV.E fully supports this statement, as shown in Eqn. (42) and Figure 25. The integrated shock-induced vorticity $\Gamma$ mainly depends on $\Pi^*$. Going beyond Ref. [12], the relationship between the generated vorticity and $\Pi^*$ is now derived from first principles and in closed form.

The practical consequence is that increasing the amount of shock-induced vorticity requires increasing $\Pi^*$. This can be achieved by increasing the pulse energy and reducing the gap length of the discharge, as suggested by the pulse model of Section II (Eqn. (2) for non-thermal sparks and Figures 5–7 for thermal sparks).

In Ref. [12], we already showed that, when $\Pi^* \lesssim 8-12$, the interpulse flow is in the diffusive regime. Otherwise, the regime is recirculating.

However, this $\Pi^*$-criterion does not predict the state of the gas in the inter-electrode region at the end of the post-discharge. It only predicts whether recirculation or diffusion dominates the cooling of the inter-electrode region in the ideal case where the interpulse duration is infinite. In practice, the interpulse duration in NRP discharges is finite. In some cases, to improve the efficiency of NRP discharges, the next pulse should be applied to fresh gas rather than already processed gas. Therefore, it is important to determine which phenomena (blast wave expansion, recirculation, or radial diffusion) influence the interpulse fluid dynamics, and whether any of these can effectively bring the gas close to ambient conditions before the next pulse occurs. This question is the core of the following Subsection.
## VII.B. Post-discharge cooling efficiency

A simple criterion to quantify the cooling efficiency is the temperature at the end of the interpulse, $T(t = 1/f)$. Cooling is efficient if $T(t = 1/f)$ is close to the ambient temperature $T_{out}$. This is expressed through the following equation:

$$T(t = 1/f) < \zeta T_{out} \qquad (54)$$

The choice of $\zeta$ is free and left to the model user, but $\zeta$ is typically $\ge 1$. Based on this equation, we will use the flow dynamics model to infer cooling criteria in Subsection VII.B.1. These criteria involve the discharge parameters ($T_k$, $T_{out}$, $f$, $L$, $R_d$). Then, we will test them against CFD simulations in Subsection VII.B.2.

### VII.B.1. Three cooling criteria

In this Subsection, we build simple criteria to determine which phenomenon (expansion, diffusion, or recirculation) is responsible for returning the inter-electrode region close to ambient conditions, as given by Eqn. (54). First and foremost, the following formulas require the blast wave to have relaxed before the end of the interpulse. This condition reads:

$$\tau_{\text{rarefaction}} < 1/f \qquad (55)$$

$$T_k > \frac{(fR_d)^2 M_k}{\gamma_k R} \approx \frac{(fR_d)^2 M_{\text{out}}}{\gamma_{\text{out}} R} \qquad (56)$$

and is generally met for typical NRP discharges, since they operate at typical repetition frequencies in the range 1–300 kHz.

The conditions under which the *blast wave expansion* nearly restores the inter-electrode region to the surrounding conditions are derived by inserting Eqn. (15) into Eqn. (54). The criterion reads:

$$\Pi^* = \frac{T_k}{T_{\text{out}}} < \zeta \gamma \qquad (57)$$

The conditions under which the combination of *expansion* and *radial diffusion* nearly restores the inter-electrode region to the surrounding conditions are obtained by inserting Eqn. (30) into Eqn. (54). We cannot directly get an analytical cooling criterion from the refined diffusion submodel. We made four simplifications to achieve this. We postulate a value of $\xi_{\text{expand}}$. We consider that the gases are in thermal equilibrium at the beginning
of the interpulse, as in the CFD simulations. We neglect heat penetration, and we assume a constant thermal diffusivity. Then, the criterion reads:

$$ \Pi^* = \frac{T_k}{T_{\text{out}}} < \left[ \frac{(\zeta - 1)}{2 \sum_{n=1}^{+\infty} \frac{\exp(-\alpha \omega_n^2 / f)}{\Omega_n J_1(\Omega_n)}} + 1 \right]^\gamma \quad (58) $$

In practice, $f$ is approximately tens of kHz. Thus, the infinite series in Eqn. (58) is well approximated by its first term only and simplifies to:

$$ \Pi^* = \frac{T_k}{T_{\text{out}}} < \left[ \frac{(\zeta - 1)\Omega_1 J_1(\Omega_1)}{2 \exp(-\alpha \omega_1^2 / f)} + 1 \right]^\gamma \quad (59) $$

where $\Omega_1 = 2.405$, $J_1(\Omega_1) = 0.519$.

Finally, the conditions under which the combination of expansion and recirculation nearly returns the inter-electrode region to its surrounding conditions are given by the following equation:

$$ 1/f > \tau_{\text{vortex}} \quad (60) $$

In other words, the interpulse duration should be longer than the recirculation timescale. Inserting Eqns. (42) and (44) into (60) leads to:

$$ \Pi^* = \frac{T_k}{T_{\text{out}}} > 1 + \left( \frac{-\beta_3 - \sqrt{\beta_3^2 - 4\beta_5(\beta_6 - \ln(fL^2))}}{2\beta_5} \right)^{\frac{1}{\beta_4}} \quad (61) $$

In summary, Eqns. (57), (59), and (61) are three simple inequalities that indicate which flow dynamics effect (blast wave expansion, diffusion, or recirculation) is responsible for returning the inter-electrode region to the surrounding conditions. They involve the non-dimensional parameter $\Pi^* = T_k/T_{\text{out}}$ and indicate which parameters must be modified to favor cooling by expansion, diffusion, or recirculation.

## VII.B.2. Comparison with CFD simulations

The three “cooling criteria” (Eqns. (57), (59), and (61)) are now compared with results from CFD simulations. We start by setting the pulse frequency to $f = 5$ kHz. For each discharge configuration ($T_k, L$), where $T_k$ and $L$ range over 300–6000 K and 1–5 mm, we compute the inter-electrode mean temperature $\bar{T}_{\text{end}}$ at the end of the post-discharge (Eqns. (10) and (49)). We repeat this procedure for six different pulse frequencies (5, 10, 15, 20, 30, 50 kHz). We obtain the colormaps of temperature $\bar{T}_{\text{end}}$ shown in Figure 35. The boundaries corresponding to the three “cooling criteria” (Eqns. (57), (59), (61)) with $\zeta = 2$ (i.e., $\bar{T}(1/f) \le 600$ K) are superposed on the graph. Here, we assumed $\xi_{\text{expand}} \approx 3$ in order to compute the diffusion cooling criteria.
First, we note that the cold region A, which characterizes complete recirculation cooling, is consistent with the prediction of the “recirculation cooling criterion” across all pulse frequencies because region A<fnref n="6" /> lies above the white solid line representing Eqn. (61). The slope of this line, with respect to L, is strictly positive. Indeed, the longer L, the higher $T_k/T_{out}$ is required to favor recirculation by the ISIV mechanism during the same interpulse period. When f increases (the pulse period decreases), the solid line shifts toward smaller L values because the vortex has less time to expel the hot gas out of the inter-electrode region. The slope with respect to L is also steeper, as predicted by Eqn. (61).

Second, we observe that the “diffusion cooling criterion” Eqn. (59) (white dashed line) correctly separates the region D, which characterizes diffusion cooling, from the remaining parametric space. This is particularly clear in the subpanels $f = 15, 20, 30, 50$ kHz of Figure 35.<fnref n="7" /> Diffusion being radial, the corresponding dashed line is horizontal (no dependence on the gap length).

Third, we remark that the expansion criterion (dotted line, Eqn (57)) is difficult to satisfy; a very low initial temperature is necessary to cool the kernel only by expansion. In addition, it does not depend on the frequency because the rarefaction wave crossed the kernel way earlier than the interpulse end.

In conclusion, Figure 35 demonstrates the ability of the “cooling criteria” (Eqns. (57), (59), (61)) to predict whether the inter-electrode region has returned to ambient conditions before the next pulse. It also predicts the flow dynamics phenomenon, whether blast wave expansion, radial diffusion, or shock-induced recirculation, is responsible for the cooling.

<footnote n="6">⁶ Labels A–D are defined in Section V.B and in Ref. [12].</footnote>

<footnote n="7">⁷ This is particularly obvious in the panels of Figure 35 which involve configurations where the plasma kernel cools by diffusion and regions where it does not, as this results in a strong color contrast. It is less obvious for the panels $f = 5$ and 10 kHz: the period is long enough for almost the entire parametric subspace to verify the “diffusion cooling criterion”, resulting in a low contrast colormap.</footnote>
<figure><img src="image_37.png" /><figcaption>Figure 35. Mean temperature along the inter-electrode axis at the end of the interpulse ($t = 1/f$). Dotted, dashed and solid lines represent the “cooling criteria” (57), (59), and (61) for $\zeta = 2$.</figcaption></figure>

## VIII. Comparison of the model predictions with experiments and CFD results from the literature

In this Section, we test the pulse and interpulse models against different results from the literature. We select:

i. the temperature measurements of Lo *et al.* [78] throughout the nanosecond pulse and the subsequent interpulse,

ii. the temperature measured at the end of successive interpulses by Blanchard *et al.* [79] and later by Gablier *et al.* [80] in a laboratory-scale burner,

iii. the temperature profiles extracted from the CFD simulation of Castela *et al.* [25], which reproduce the NRP discharges studied by Rusterholtz *et al.* [35].

This Section complements Section V, in which we only validated the interpulse model without VT relaxation.

### VIII.A. Comparison to the experiment of Lo *et al.* (2014)

In the experiment described in [78], the temporal profile of temperature was measured by Optical Emission Spectroscopy and Spontaneous Raman Scattering during a plasma pulse and an interpulse. The plasma discharge was generated in air between a pin and a plane. The presence of a plane electrode could perturb the recirculation and result in heat loss to the wall, two phenomena that are not included in our engineering model,
which is designed for pin-to-pin discharges. We verified that, under the condition of [78], pulse heating is not high enough to generate recirculation ($\Pi^* \approx 3$ does not reach the threshold of 8–12 developed in [12]). We also selected temperatures measured at the two locations closest to the pin (3.25 and 5 mm from the plane). Thus, the impact of the thermal diffusion on the plane is expected to be negligible compared to the radial one for the selected data, and our model is applicable in the first approximation.

The discharge parameters used for our model are given in Table 3. We adjusted the ultrafast and slow heating fractions proposed in [25] by 5% to improve the agreement with the experiment. The VT relaxation timescale is taken from Figure 16 for air at 3000 K.

<table><caption>Table 3. Main input parameters of the discharge in the experiment of Lo et al.</caption><thead><tr><th>Input parameters</th><th>Case from [78]</th></tr></thead><tbody><tr><td>θ<sub>uh</sub>, θ<sub>sh</sub></td><td>0.25, 0.45</td></tr><tr><td>ΔE<sub>pulse</sub></td><td>20 mJ</td></tr><tr><td>T<sub>out</sub>, P<sub>out</sub></td><td>300 K, 1 bar</td></tr><tr><td>M, γ</td><td>28 g/mol, 7/5</td></tr><tr><td>L, R<sub>d</sub></td><td>6.5 mm, 700 µm</td></tr><tr><td>α</td><td>Temperature-dependent relation (see [46])</td></tr><tr><td>Θ<sub>vib</sub></td><td>~3400 K (see [37])</td></tr><tr><td>τ<sub>VT</sub></td><td>~10 µs (see Figure 16)</td></tr></tbody></table>

The model's prediction is plotted in red in Figure 36. The red swath illustrates the effect of changing the discharge radius by ±50 µm, highlighting the impact of the discharge radius on temperature. This is mainly due to the presence of $R_d^2$ in the volumetric heating and diffusion equations. The blue and cyan symbols correspond to the measurements of [78].

As shown in Figure 36, the model captures the global trend of the temperature profile. We considered that the ultrafast heating is linear and lasts 30 ns to fit the experimental profile. The time resolution of measurements between 1 and 10 µs does not permit the capture of any potential decrease in temperature due to the expansion. The second increase in temperature, due to VT relaxation, fits in terms of amplitude (+800 K) and timescale (10 µs). Finally, the model underestimates the decrease in temperature due to diffusion. At $t = 1$ ms, the model predicts a temperature of 1000 K, whereas the
measured temperature is between 700 and 850 K. The difference between the model and the experimental data is higher near the wall ($z = 3.25$ mm), which could be due to heat losses at the wall.

<figure><img src="image_38.png" /><figcaption>Figure 36. Comparison between the model prediction and data from Lo et al. [78]. The red line is the model prediction for the temporal gas temperature profile (pulse model and $T_4$ law, i.e., average temperature of the gap). The blue dots and the cyan triangles correspond to the measurements of Lo et al. [78] at two different locations $z$ along the discharge axis. The red swath illustrates the dependence of the model on the discharge radius, when changing its value by $\pm 50$ µm.</figcaption></figure>

### VIII.B. Comparison to the simulation of Castela et al. (2016)

In this Subsection, we model the NRP discharges experiment of Rusterholtz et al. [35], which was simulated by Castela et al [25]. Castela et al used a reactive and compressible DNS code, coupled with a plasma model. We collected the discharge parameters used in [25] in Table 4, and injected these parameters into our analytical model. The pulse and interpulse models are applied iteratively. To this end, we couple the pulse and interpulse models at each step $n$ and $n + 1$. The initial temperature of the channel for pulse $n + 1$ is the final temperature of the interpulse $n$, i.e. $T_{\text{ini}}^{n+1} = T_4^n(t = 1/f)$.
<table><caption>Table 4. Main input parameters of the discharge in the simulation of Castela et al.</caption><thead><tr><th>Input parameters</th><th>Case from [25]</th></tr></thead><tbody><tr><td>θ<sub>uh</sub>, θ<sub>sh</sub></td><td>0.2, 0.45</td></tr><tr><td>ΔE<sub>pulse</sub></td><td>0.7 mJ</td></tr><tr><td>f</td><td>10 kHz</td></tr><tr><td>T<sub>out</sub>, P<sub>out</sub></td><td>300 K, 1 bar</td></tr><tr><td>M, γ</td><td>28 g/mol, 7/5</td></tr><tr><td>L, R<sub>d</sub></td><td>4 mm, 225 μm</td></tr><tr><td>α</td><td>Temperature-dependent relation (see [46])</td></tr><tr><td>Θ<sub>vib</sub></td><td>~3400 K (see [37])</td></tr><tr><td>τ<sub>VT</sub></td><td>~3 μs (see [57])</td></tr></tbody></table>

Figure 37 shows the prediction of our analytical model (in red) and the results of the literature CFD simulation (in blue). Here again, the red swath indicates the model's sensitivity to the discharge radius, whose value is increased or decreased by 25 µm.

At first glance, our analytical model captures correctly the temperature trends while being significantly faster to execute than the CFD simulation ($\tau_{\text{model}}^{\text{CPU}} < 1$ mn and $\tau_{\text{CFD}}^{\text{CPU}} \sim 10^6$ h). The cyclic pattern is reproduced with errors on maxima and minima below 10%. However, we can note several differences. (i) The gas heating during the first pulses is stronger in the model than in the simulation. (ii) The VT transfer is quicker in the model than in the simulation. It is due to a significant difference between the model and the simulation, $\tau_{\text{VT}}$ is fixed in the model, whereas it results from temperature and densities in each cell during the CFD simulation. We also notice that the $\Pi^*$ number is low (<10), so cooling is mainly driven by diffusion.
<figure><img src="image_39.png" /><figcaption>Figure 37. Comparison of the model prediction with the DNS data of Castela *et al.* [25]. The red line is the model predictions for the temporal temperature profile of the gas (pulse model and $T_4$ law, i.e., the average temperature of the gap). The blue line indicates the result of the simulation of [25]. The red swath illustrates the dependence of the model on the discharge radius, when changing its value by ±25 µm.</figcaption></figure>

### VIII.C. Comparison to the two experiments of Blanchard *et al* (2022) and Gablier *et al* (2025)

In this Subsection, we use our model to reproduce the conditions of two experiments made in a laboratory-scale burner called Mini-PAC (see [81] for more details about this facility). Two series of experiments were conducted in methane-air flames with two different inter-electrode gap distances (1 and 5 mm) but the same deposited energy (1.8 mJ/pulse). The temperature was measured by emission spectroscopy at the beginning of each pulse. In the experiments of Blanchard *et al.* [79], conducted with a gap distance of 5 mm, the temperature midway between the electrodes increased progressively and stabilized after a few tens of pulses. The discharge radius was measured by OES and was approximately 350 µm (see Fig. 4.4 of [82]). In the experiments of Gablier *et al.* [80], the inter-electrode gap was 1 mm, the measured radius was about 250 µm, and the degree of gas heating at midgap was found to stay constant at a value close to the adiabatic temperature of the flame without plasma.

The different discharge parameters required to run the model in the two configurations are detailed in Table 5. Regarding the volumetric energy in the second case, the non-thermal spark pulse model predicted values of $T_k$ above 10,000 K. The value of $T_k$ is indicative of a thermal spark, and this is further confirmed by the fact that the emission spectra measured by Gablier *et al.* features $N^+$ lines and no $N_2(C-B)$ emission. For these reasons, we used the thermal spark model (assuming infinitely fast chemistry).
<table><caption>Table 5. Input parameters of the model for the two experimental cases simulated in this Subsection.</caption><thead><tr><th>Input parameters</th><th>Blanchard et al. [79]</th><th>Gablier et al. [80]</th></tr></thead><tbody><tr><td>θ<sub>uh</sub>, θ<sub>sh</sub></td><td>0.2, 0.45<br>(see Tab. 6.2 and 6.3 of [79])</td><td>Not defined for a<br>thermal spark model</td></tr><tr><td>ΔE<sub>pulse</sub></td><td colspan="2">1.8 mJ</td></tr><tr><td>f</td><td colspan="2">20 kHz</td></tr><tr><td>T<sub>out</sub>, P<sub>out</sub></td><td colspan="2">1992 K (adiabatic flame temperature with an equivalence ratio<br>of 0.8), 1 bar</td></tr><tr><td>M, γ</td><td>28 g/mol, 7/5</td><td>Computed at LTE from T<sub>k</sub><br>and P<sub>k</sub></td></tr><tr><td>L, R<sub>d</sub></td><td>5 mm, 350 μm</td><td>1 mm, 250 μm</td></tr><tr><td>α</td><td colspan="2">Temperature-dependent relation (see [46])</td></tr><tr><td>Θ<sub>vib</sub></td><td>~3400 K (see [37])</td><td>Not needed for a thermal<br>spark</td></tr><tr><td>τ<sub>VT</sub></td><td>~3 μs (see [57])</td><td></td></tr></tbody></table>

Figure 38 and Figure 39 compare the modeled mean inter-electrode temperature $T_4$ with the measured temperature for the two cases studied by Blanchard *et al.* [79] and Gablier *et al.* [80]. In the limit cycle, the temperature increases rapidly during the pulse as a result of ultrafast heating, up to 5000 or 13,000 K (depending on the case). This difference of $T_k$ is due to the difference in radius between the two cases. Then, the temperature falls to $T_{rarefaction} \approx 4000$ K or 9,500 K. In the first case (non-thermal spark), the temperature increases again due to the slow heating up to 5500 K. Then, the temperature decreases due to diffusion in the first case and recirculation in the second case.<fnref n="8" />

Similarly to Subsection VIII.B, the increase in temperature during the first pulses is not correctly predicted for the non-thermal spark case. This may be due to the smaller energy injected during the first pulse (the initial temperature between electrodes being smaller,

<footnote n="8">⁸ For the second case, the diffusion and expansion models used in the interpulse model are outside their range of validity. However, recirculation is dominant during the cooling phase, and its submodel is valid, even for thermal sparks.</footnote>
the reduced electric field is also lower). The value of 1.8 mJ chosen here is a mean value measured during the asymptotic temperature cycle.

<figure><img src="image_40.png" /><figcaption>Figure 38. Comparison between the model prediction and data from Blanchard *et al.* [79]. The red line is the model predictions for the temporal temperature profile of the gas (pulse model and $T_4$ law, i.e., the average temperature of the gap). The blue symbols indicate the measurements of [79]. The red swath illustrates the dependence of the model on the discharge radius, when changing its value by ±25 µm.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 39. Comparison between the model prediction and data from Gablier *et al.* [80]. The red line is the model predictions for the temporal temperature profile of the gas (pulse model $T_4$ law, i.e., the average temperature of the gap). The blue symbols indicate the measurements of [80]. The red swath illustrates the dependence of the model on the discharge radius, when changing its value by ±25 µm.</figcaption></figure>

Using this iterative approach, the pulse and interpulse models reach their limit because the inherent errors of such a simple model accumulate over time. This error increases when the pulse frequency increases. In addition, the discharge parameters may change between pulses, but we keep them constant along the iterations. For instance, residual heat and radical species may affect the breakdown process, altering the pulse energy and the fast/slow heating fractions. The structure of our model enables users to couple it iteratively with a breakdown/electrical model, thereby improving its fidelity. However, modeling the breakdown and the electrical circuit is beyond the scope of this work.
The pulse and interpulse models provide a very good and rapid estimate of the discharge temperature evolution over time for nanosecond discharges.

## IX. Conclusions

We developed an engineering model to predict the mean temperature and pressure in the inter-electrode region during and following a nanosecond pulse, creating a cylindrical filamentary hot channel. The model comprises a pulse model and a post-discharge flow dynamics model.

The pulse model has two versions, one for non-thermal sparks and the other for thermal sparks, which provide estimates of the pressure and temperature of the plasma kernel right after the pulse ($P_k, T_k$).

The post-discharge model is divided into four components in series: blast wave expansion, VT energy relaxation, radial heat diffusion, and shock-induced recirculation. For each component, simple analytical relations are provided for the temporal evolutions of temperature and pressure along the inter-electrode axis. First, ($P_1, T_1$) isolates the role of the expansion that occurs during the first hundreds of nanoseconds after the discharge. Second, ($P_2, T_2$) adds a correction linked to any thermal imbalance between the trans-rotational temperature (on which we focus in this paper) and a vibrational temperature. Third, ($P_3, T_3$) takes into account the expansion, thermal relaxation towards equilibrium, and heat diffusion. They can be used to model the properties of the kernel if there is no significant recirculation due to the shock wave. It is also an approximation of the properties of the expelled gases if recirculation is significant. Fourth, ($P_4, T_4$) takes into account the three main flow dynamics phenomena: expansion, heat diffusion, and recirculation, as well as thermal relaxation towards thermal equilibrium. It is suited to compute the condition inside the gap at the end of an interpulse, just before the following pulse, for repetitive discharges.

The predictions of the reduced-order model match those of the CFD simulations for a wide range of non-thermal discharge conditions ($T_k < 6000$ K, $P_k < 20$ bar). A refinement for thermal sparks is also proposed.

We derived three “cooling criteria” that predict whether the return of the inter-electrode region to ambient conditions at the end of the interpulse is partial or nearly complete. These cooling criteria were tested against CFD simulations and offer additional insights beyond the previously established $\Pi^*$-criteria used to predict the recirculation or diffusive regimes.

Finally, predictions of the pulse and interpulse models for temporal temperature profiles and cumulative heating were validated by comparisons with laboratory experiments and simulations from the literature.
The model offers two major advantages: it is easy to use and has a very short computation time (a few seconds). In contrast, CFD simulations such as those presented in Section III require significant expertise and $10^3-10^6$ CPU hours on a high-performance computing (HPC) cluster.

The model is not intended to replace such simulations but rather to capture the coupled mechanisms affecting temperature and pressure during nanosecond discharges and post-discharges, and to explore the dependencies on many parameters rapidly. Additionally, each brick of the model can be used separately for different purposes, such as quantifying shock-induced vorticity or estimating temperature increases during thermal sparks.

## Acknowledgments

The authors are grateful to Prof. Aymeric Vié and Dr. Christian Tenaud for valuable discussions. We also would like to express our thanks to Dr. Christian Tenaud for providing the 1D Euler solver from which the expansion relations were derived in the expansion submodel (Subsection IV.B). This work has received funding from the Ministerial Allocation for Polytechnique students (AMX) and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC GreenBlue grant agreement No.101021538). A CC-BY 4.0 public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions: https://creativecommons.org/licenses/by/4.0/.

## Appendices

### A Sensitivity of the recirculation model to the initial shape of the plasma kernel

In this Appendix, we investigate the sensitivity of the amount of vorticity calculated by GSD to the initial shape of the plasma kernel. We show that the engineering model is not significantly sensitive to the initial shape of the kernel. This is consistent with the CFD simulations performed in Appendix A of [12].

We selected a reference GSD simulation case ($L = 3$ mm, $T_k = 6000$ K, $T_{out} = 300$ K, $R_d = 100$ µm). We tested four kernel shapes with the same dimensions ($R_d, L/2$) but different contours, as drawn in Figure 40. The reference case has an elliptic shape; the other cases have hyperbolic, cubic, and parabolic shapes. Additionally, two other simulations were performed, in which the discharge radius was either divided or multiplied by a factor of 2. GSD simulations were performed for all these cases, allowing
the shock to propagate until it generated no more vorticity, and then the vorticity integral was computed.

<figure><img src="image_42.png" /><figcaption>Figure 40. Shapes of the plasma kernel used to initialize the GSD simulations. Shapes with various profiles (elliptic, cubic, hyperbolic, parabolic) are shown in green. Elliptic shapes with various radii (50, 100, and 200 µm) are shown in purple. The main graph zooms in on the near-electrode region. The full profiles are shown in the inset. Taken from [12].</figcaption></figure>

We compare the recirculation $|Γ|$ of the six cases in Figure 41. We see that the values of $|Γ|$ have the same order of magnitude. It is interesting to compare this result with Figure 25, which shows that a small increase of +1 in the temperature ratio $Π^* \equiv T_k/T_{out}$ can affect this integral by several orders of magnitude, whereas a change in kernel shape has practically no impact. The temperature ratio is by far the dominant factor of vorticity creation, and therefore, a fit of the recirculation based on $Π^*$ is convenient for engineering models.
<figure><img src="image_43.png" /><figcaption>Figure 41. Variation of the vorticity jump integral computed from the GSD simulations in a reference configuration for different discharge parameters (kernel shape, kernel radius). The discharge configuration corresponds to $L = 3$ mm, $T_k = 6000$ K.</figcaption></figure>

## B Overview of the assumptions and formulas of the pulse and interpulse models

The two following Tables summarize the assumptions, equations, and tables for the different models and submodels.

<table><caption>Table 6. Assumptions and equations of the pulse model.</caption><thead><tr><th>Submodel</th><th>Assumptions</th><th>Main equations</th></tr></thead><tbody><tr><td>Non-thermal sparks</td><td>(a1) Isochoric pulse heating.<br>(a2) The gas molar mass remains unchanged during the pulse.<br>(a3) The vibrational heat capacity has a linear-constant profile.</td><td>Eqns. (2) and (4).<br>Eqn. (7) also returns the vibrational temperature.</td></tr><tr><td>Thermal sparks</td><td>(a1) Isochoric pulse heating.<br>(a2) The pulse energy is provided to the gas, which is either at LTE or chemically frozen.</td><td>Eqns. (8) and (9).<br>The system is solved for air, CO<sub>2</sub>, and H<sub>2</sub> in Figure 5, Figure 6, and Figure 7.</td></tr></tbody></table>
<table><caption>Table 7. Assumptions and equations of the interpulse model.</caption><thead><tr><th>Submodel</th><th>Assumptions</th><th>Main equations</th></tr></thead><tbody><tr><td>Expansion<br>submodel</td><td>(a1) 1D cylindrical configuration.<br>(a2) Inviscid flow.<br>(a3) Calorically perfect gas.<br>(a4) Constant adiabatic index and<br>molar mass inside and outside the<br>kernel.<br>(a5) Ideal gas.<br>(a6) Gas initially at rest.<br>(a7) Negligible heat diffusion.<br>(a8) Stepwise initial profiles for T<br>and P.<br>(a9) Reversible and adiabatic<br>expansion in the enclosed domain.</td><td>Eqns. (17) and (18),<br>with the use of<br>Eqns. (14), (15), (16).</td></tr><tr><td>Energy relaxation<br>submodel for<br>molecular gases</td><td>(a1) Isobaric VT energy exchange.<br>(a2) Constant τ<sub>VT</sub>.</td><td>Eqns. (23) and (24)<br>with the use of<br>Eqns. (7), (20), and<br>(22).</td></tr></tbody></table>
<table><tr><td>Diffusion submodel</td><td>(a1) Plasma remained within an enclosed surface during the expansion.<br/>(a2) Isentropic expansion.<br/>(a3) Constant number of particles.<br/>(a4) Ideal gas.<br/>(a5) 1D cylindrical configuration.<br/>(a6) Uniform pressure (equal to P<sub>out</sub>), incompressible flow.<br/>(a7) Calorically perfect gas and constant thermal diffusivity.<br/>(a8) Gas initially at rest.<br/>(a9) Splitting of the timescales associated with VT energy transfer and expansion. (loosely coupling only).<br/>(a10) Rectangular initial profiles for T.<br/>(a11) Dirichlet boundary condition, T(ξ<sub>exp</sub>R<sub>d</sub>, t) = T<sub>out</sub>.<br/><br/>The refined model does not presuppose the assumptions (a7) and (a11).<br/><br/>Eqns. (30) and (31), with the use of Eqn. (25) and Table 2.<br/><br/>Eqns. (31) and (34) with the use of Eqns. (25), (33), (35), and Table 2.</td></tr><tr><td>Recirculation submodel</td><td>(a1) All necessary hypotheses to build the Euler system of equations.<br/>(a2) Negligible external flow velocity (upstream flow perpendicular to the shock wave).<br/>(a3) Calorically perfect gas with constant adiabatic index and molar mass.<br/>(a4) Low impact of the rarefaction waves on the leading shock wave.<br/>(a5) 2D cylindrical configuration.<br/>(a6) Elliptic kernel with stepwise uniform temperature/profile in the kernel and outside.</td></tr></table>
<table><tr><td>The refined model does not presuppose (a2).</td><td>Same and adding Eqns. (50) and (51)</td></tr></table>

## References

[1] Ceppelli, M., Salden, T. P. W., Martini, L. M., Dilecce, G., and Tosi, P., “Time-Resolved Optical Emission Spectroscopy in CO2 Nanosecond Pulsed Discharges,” *Plasma Sources Science and Technology*, Vol. 30, No. 11, 2021, p. 115010. https://doi.org/10.1088/1361-6595/ac2411

[2] Post, A., Budde, M., Vervloedt, S. C. L., Engeln, R., and Huiskamp, T., “Nanosecond Repetitively Pulsed Plasmas with MHz Bursts for CO2 Dissociation,” *Journal of Physics D: Applied Physics*, Vol. 57, No. 30, 2024, p. 305203. https://doi.org/10.1088/1361-6463/ad40bc

[3] Yong, T., Zhong, H., Pannier, E., Laux, C., and Cappelli, M. A., “High-Pressure CO2 Dissociation with Nanosecond Pulsed Discharges,” *Plasma Sources Science and Technology*, Vol. 32, No. 11, 2023, p. 115012. https://doi.org/10.1088/1361-6595/ad066e

[4] de Klerk, A., “Fischer-Tropsch Process,” *Kirk-Othmer Encyclopedia of Chemical Technology*, John Wiley & Sons, Ltd, 2013, pp. 1-20. https://doi.org/10.1002/0471238961.fiscdekl.a01

[5] Starikovskiy, A., and Aleksandrov, N., “Plasma-Assisted Ignition and Combustion,” *Progress in Energy and Combustion Science*, Vol. 39, No. 1, 2013, pp. 61-110. https://doi.org/10.1016/j.pecs.2012.05.003

[6] Malé, Q., Shcherbanev, S., Impagnatiello, M., and Noiray, N., “Stabilization of a Thermoacoustically Unstable Sequential Combustor Using Non-Equilibrium Plasma: Large Eddy Simulation and Experiments,” *Proceedings of the Combustion Institute*, Vol. 40, No. 1, 2024, p. 105277. https://doi.org/10.1016/j.proci.2024.105277

[7] Pavan, C. A., Shanbhogue, S. J., Weibel, D. E., Gomez Del Campo, F., Ghoniem, A. F., and Guerra-Garcia, C., “Dynamic Response of Nanosecond Repetitively Pulsed Discharges to Combustion Dynamics: Regime Transitions Driven by Flame Oscillations,” *Plasma Sources Science and Technology*, Vol. 33, No. 2, 2024, p. 025016. https://doi.org/10.1088/1361-6595/ad227d

[8] Moeck, J., Lacoste, D., Laux, C., and Paschereit, C., “Control of Combustion Dynamics in a Swirl-Stabilized Combustor with Nanosecond Repetitively Pulsed Discharges,” In *51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition*, Vol. AIAA 2013-565, Grapevine (Dallas/Ft. Worth Region), Texas, 2013. https://doi.org/10.2514/6.2013-565
[9] Pilla, G., Galley, D., Lacoste, D. A., Lacas, F., Veynante, D., and Laux, C. O., "Stabilization of a Turbulent Premixed Flame Using a Nanosecond Repetitively Pulsed Plasma," *IEEE Transactions on Plasma Science*, Vol. 34, No. 6, 2006, pp. 2471-2477.
https://doi.org/10.1109/TPS.2006.886081

[10] Sabatino, F. D., and Lacoste, D. A., "Enhancement of the Lean Stability and Blow-off Limits of Methane-Air Swirl Flames at Elevated Pressures by Nanosecond Repetitively Pulsed Discharges," *Journal of Physics D: Applied Physics*, Vol. 53, No. 35, 2020, p. 355201. https://doi.org/10.1088/1361-6463/ab8f54

[11] Blanchard, V. P., Scouflaire, P., Laux, C. O., and Ducruix, S., "Combustion Performance of Plasma-Stabilized Lean Flames in a Gas Turbine Model Combustor," *Applications in Energy and Combustion Science*, Vol. 15, 2023, p. 100158.
https://doi.org/10.1016/j.jaecs.2023.100158

[12] Roger, E., Mariotto, P., and Laux, C. O., "Origin of the Recirculation Flow Pattern Induced by Nanosecond Discharges and Criterion for Its Development," *Journal of Physics D: Applied Physics*, 2025. https://doi.org/10.1088/1361-6463/adb1f0

[13] Pai, D. Z., Stancu, G. D., Lacoste, D. A., and Laux, C. O., "Nanosecond Repetitively Pulsed Discharges in Air at Atmospheric Pressure—the Glow Regime," *Plasma Sources Science and Technology*, Vol. 18, No. 4, 2009, p. 045030.
https://doi.org/10.1088/0963-0252/18/4/045030

[14] Pai, D. Z., Lacoste, D. A., and Laux, C. O., "Nanosecond Repetitively Pulsed Discharges in Air at Atmospheric Pressure—the Spark Regime," *Plasma Sources Science and Technology*, Vol. 19, No. 6, 2010, p. 065015. https://doi.org/10.1088/0963-0252/19/6/065015

[15] Minesi, N., Stepanyan, S., Mariotto, P., Stancu, G. D., and Laux, C. O., "Fully Ionized Nanosecond Discharges in Air: The Thermal Spark," *Plasma Sources Science and Technology*, Vol. 29, No. 8, 2020, p. 085003. https://doi.org/10.1088/1361-6595/ab94d3

[16] Singh, B., Rajendran, L. K., Zhang, J., Vlachos, P. P., and Bane, S. P. M., "Vortex Rings Drive Entrainment and Cooling in Flow Induced by a Spark Discharge," *Physical Review Fluids*, Vol. 5, No. 11, 2020, p. 114501.
https://doi.org/10.1103/PhysRevFluids.5.114501

[17] Kono, M., Niu, K., Tsukamoto, T., and Ujiie, Y., "Mechanism of Flame Kernel Formation Produced by Short Duration Sparks," *Symposium (International) on Combustion*, Vol. 22, No. 1, 1989, pp. 1643-1649. https://doi.org/10.1016/S0082-0784(89)80176-6
[18] Thiele, M., Warnatz, J., and Maas, U., “Geometrical Study of Spark Ignition in Two Dimensions,” *Combustion Theory and Modelling*, Vol. 4, No. 4, 2000, pp. 413–434.
https://doi.org/10.1088/1364-7830/4/4/303

[19] Adams, S., Dumitrache, C., and Murray, C., “Gas Density Evolution and Recirculation after a Nanosecond Pulsed Discharge: Experimental and Computational Results,” *Experiments in Fluids*, Vol. 64, No. 12, 2023, p. 192.
https://doi.org/10.1007/s00348-023-03724-6

[20] Hagelaar, G. J. M., and Pitchford, L. C., “Solving the Boltzmann Equation to Obtain Electron Transport Coefficients and Rate Coefficients for Fluid Models,” *Plasma Sources Science and Technology*, Vol. 14, No. 4, 2005, pp. 722–733.
https://doi.org/10.1088/0963-0252/14/4/011

[21] Pancheshnyi, S., Eismann, B., Hagelaar, G., and Pitchford, L., “ZDPlasKin: A New Tool for Plasmachemical Simulations,” 2008.

[22] Shao, X., Lacoste, D. A., and Im, H. G., “ChemPlasKin: A General-Purpose Program for Unified Gas and Plasma Kinetics Simulations,” *Applications in Energy and Combustion Science*, Vol. 19, 2024, p. 100280.
https://doi.org/10.1016/j.jaecs.2024.100280

[23] Barleon, N., Cheng, L., Cuenot, B., Vermorel, O., and Bourdon, A., “Investigation of the Impact of NRP Discharge Frequency on the Ignition of a Lean Methane-Air Mixture Using Fully Coupled Plasma-Combustion Numerical Simulations,” *Proceedings of the Combustion Institute*, 2022. https://doi.org/10.1016/j.proci.2022.07.046

[24] Bechane, Y., and Fiorina, B., “A Numerical Investigation of Plasma-Assisted Ignition by a Burst of Nanosecond Repetitively Pulsed Discharges,” *Combustion and Flame*, Vol. 259, 2024, p. 113106. https://doi.org/10.1016/j.combustflame.2023.113106

[25] Castela, M., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., and Laux, C. O., “Modelling the Impact of Non-Equilibrium Discharges on Reactive Mixtures for Simulations of Plasma-Assisted Ignition in Turbulent Flows,” *Combustion and Flame*, Vol. 166, 2016, pp. 133–147. https://doi.org/10.1016/j.combustflame.2016.01.009

[26] Smith, G. P., Golden, D. M., Frenklach, M., Moriarty, N. W., Eiteneer, B., Goldenberg, M., Bowman, C. T., Hanson, R. K., Song, S., Gardiner, Jr, W. C., Lissianski, V. V., and Zhiwei, Q., “GRI-Mech Home Page.” Retrieved 29 August 2024.
http://combustion.berkeley.edu/gri-mech/

[27] Lefkowitz, J. K., Guo, P., Rousso, A., and Ju, Y., “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*, Vol. 373, No. 2048, 2015, p. 20140333.
https://doi.org/10.1098/rsta.2014.0333

[28] Popov, N. A., "Pulsed Nanosecond Discharge in Air at High Specific Deposited Energy: Fast Gas Heating and Active Particle Production," *Plasma Sources Science and Technology*, Vol. 25, No. 4, 2016, p. 044003. https://doi.org/10.1088/0963-0252/25/4/044003

[29] Guerra, V., Sá, P. A., and Loureiro, J., "Kinetic Modeling of Low-Pressure Nitrogen Discharges and Post-Discharges," *The European Physical Journal Applied Physics*, Vol. 28, No. 2, 2004, pp. 125-152. https://doi.org/10.1051/epjap:2004188

[30] Adamovich, I. V., Li, T., and Lempert, W. R., "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*, Vol. 373, No. 2048, 2015, p. 20140336. https://doi.org/10.1098/rsta.2014.0336

[31] Mariotto, P., "Kinetics of Nonequilibrium Recombining Nitrogen-Argon Plasmas," Université Paris-Saclay, PhD thesis, 2023.

[32] Kossyi, I. A., Kostinsky, A. Y., Matveyev, A. A., and Silakov, V. P., "Kinetic Scheme of the Non-Equilibrium Discharge in Nitrogen-Oxygen Mixtures," *Plasma Sources Science and Technology*, Vol. 1, No. 3, 1992, p. 207. https://doi.org/10.1088/0963-0252/1/3/011

[33] Meyer, G., and Wimmer, A., "A Thermodynamic Model for the Plasma Kernel Volume and Temperature Resulting from Spark Discharge at High Pressures," *Journal of Thermal Analysis and Calorimetry*, Vol. 133, No. 2, 2018, pp. 1195-1205. https://doi.org/10.1007/s10973-018-7169-z

[34] Xu, D. A., Shneider, M. N., Lacoste, D. A., and Laux, C. O., "Thermal and Hydrodynamic Effects of Nanosecond Discharges in Atmospheric Pressure Air," *Journal of Physics D: Applied Physics*, Vol. 47, No. 23, 2014, p. 235202. https://doi.org/10.1088/0022-3727/47/23/235202

[35] Rusterholtz, D. L., Lacoste, D. A., Stancu, G. D., Pai, D. Z., and Laux, C. O., "Ultrafast Heating and Oxygen Dissociation in Atmospheric Pressure Air by Nanosecond Repetitively Pulsed Discharges," *Journal of Physics D: Applied Physics*, Vol. 46, No. 46, 2013, p. 464010. https://doi.org/10.1088/0022-3727/46/46/464010

[36] Anderson, John. D., "Modern Compressible Flow With Historical Perspective, Third Edition," The McGraw-Hill Companies, Inc., 2003.
[37] McQuarrie, D. A., and Simon, J. D., “Physical Chemistry: A Molecular Approach,” University Science Books, Sausalito, Calif, 1997.

[38] Maillard, J., “Energy Efficiency of CO2 Conversion with Nanosecond Plasma Discharges,” Université Paris-Saclay, PhD thesis, 2024.

[39] Capitelli, M., Colonna, G., and D'Angola, A., “Fundamental Aspects of Plasma Chemical Physics: Thermodynamics,” Springer New York, New York, NY, 2012. https://doi.org/10.1007/978-1-4419-8182-0

[40] Goodwin, D. G., Moffat, H. K., Schoegl, I., Speth, R. L., and Weber, B. W., “Cantera: An Object-Oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes,” Aug 22 2023. https://doi.org/10.5281/zenodo.8137090

[41] McBride, B. J., Zehe, M. J., and Gordon, S., “NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species,” NASA/TP-2002-211556, 2002.

[42] Minesi, N. Q., Mariotto, P. B., Pannier, E., Vincent-Randonnier, A., Stancu, G. D., and Laux, C. O., “Kinetic Mechanism and Sub-Ns Measurements of the Thermal Spark in Air,” *Plasma Sources Science and Technology*, Vol. 32, No. 4, 2023, p. 044005. https://doi.org/10.1088/1361-6595/acc9da

[43] Lin, S.-C., “Cylindrical Shock Waves Produced by Instantaneous Energy Release,” *Journal of Applied Physics*, Vol. 25, No. 1, 1954, pp. 54–57. https://doi.org/10.1063/1.1721520

[44] Jones, D. L., “The Energy Parameter B for Strong Blast Wave,” NBS-Technical note, 1962.

[45] Ansys Inc, “Ansys Fluent 2022R1 Theory Guide,” 2022.

[46] Gupta, R. N., Lee, K.-P., Thompson, R. A., and Yos, J. M., “Calculations and Curve Fits of Thermodynamic and Transport Properties for Equilibrium Air to 30000 K,” NASA-RP-1260, October 1991.

[47] Toro, E. F., “Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction,” Springer, Dordrecht; New York, 2009.

[48] Roe, P. L., “Characteristic-Based Schemes for the Euler Equations,” *Annual Review of Fluid Mechanics*, Vol. 18, 1986, pp. 337–365. https://doi.org/10.1146/annurev.fl.18.010186.002005

[49] Vlases, G. C., and Jones, D. L., “Blast Waves from an Inverse Pinch Machine,” *Physics of Fluids*, Vol. 9, 1966, pp. 478–485. https://doi.org/10.1063/1.1761701
[50] Plooster, M. N., "Shock Waves from Line Sources. Numerical Solutions and Experimental Measurements," *The Physics of Fluids*, Vol. 13, No. 11, 1970, pp. 2665-2675. https://doi.org/10.1063/1.1692848

[51] Sakurai, A., "On the Propagation and Structure of the Blast Wave, I," *Journal of the Physical Society of Japan*, Vol. 8, No. 5, 1953, pp. 662-669.
https://doi.org/10.1143/JPSJ.8.662

[52] Dewey, J. M., "The Friedlander Equations," *Blast Effects*, edited by I. Sochet, Springer International Publishing, Cham, 2018, pp. 37-55. https://doi.org/10.1007/978-3-319-70831-7_3

[53] Hilsenrath, J., Beckett, C. W., Benedict, W. S., Fano, L., Hoge, H. J., Masi, J. F., Nuttall, R. L., Touloukian, Y. S., and Woolley, H. W., "Tables of Thermal Properties of Gases," National bureau of standards, 1955.

[54] Friedman, M. P., "A Simplified Analysis of Spherical and Cylindrical Blast Waves," *Journal of Fluid Mechanics*, Vol. 11, No. 1, 1961, pp. 1-15.
https://doi.org/10.1017/S0022112061000810

[55] Osipov, A. I., and Uvarov, A. V., "Kinetic and Gasdynamic Processes in Nonequilibrium Molecular Physics," *Soviet Physics Uspekhi*, Vol. 35, No. 11, 1992, p. 903.
https://doi.org/10.1070/PU1992v035n11ABEH002275

[56] Capitelli, M., Ferreira, C. M., Gordiets, B. F., and Osipov, A. I., "Plasma Kinetics in Atmospheric Gases," Springer Berlin Heidelberg, Berlin, Heidelberg, 2000.
https://doi.org/10.1007/978-3-662-04158-1

[57] Popov, N. A., "Kinetic Processes Initiated by a Nanosecond High-Current Discharge in Hot Air," *Plasma Physics Reports*, Vol. 37, No. 9, 2011, pp. 807-815.
https://doi.org/10.1134/S1063780X1108006X

[58] Blauer, J., and Nickerson, G., "A Survey of Vibrational Relaxation Rate Data for Processes Important to CO2-N2-H2O Infrared Plume Radiation," presented at the 7th Fluid and PlasmaDynamics Conference, Palo Alto, CA, U.S.A., 1974.
https://doi.org/10.2514/6.1974-536

[59] Millikan, R. C., and White, D. R., "Systematics of Vibrational Relaxation," *The Journal of Chemical Physics*, Vol. 39, No. 12, 1963, pp. 3209-3213.
https://doi.org/10.1063/1.1734182

[60] Vincenti and Kruger, "Introduction to Physical Gas Dynamics," 1975.

[61] Cengel, Yunus. A., "Heat Transfer, a Practical Approach," McGraw-Hill, 1998.
[62] Abramowitz, M., and Stegun, I. A., “Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables,” National bureau of standards applied mathematics series, 1964.

[63] Hayes, W. D., “The Vorticity Jump across a Gasdynamic Discontinuity,” *Journal of Fluid Mechanics*, Vol. 2, No. 6, 1957, pp. 595–600.
https://doi.org/10.1017/S0022112057000403

[64] Han, Z.-Y., and Yin, X.-Z., “Chapter 3.7 - Geometrical Shock Dynamics,” *Handbook of Shock Waves*, edited by G. Ben-dor, O. Igra, and T. Elperin, Academic Press, Burlington, 2001, pp. 485–552. https://doi.org/10.1016/B978-012086430-0/50011-7

[65] Ridoux, J., Lardjane, N., Monasse, L., and Coulouvrat, F., “Comparison of Geometrical Shock Dynamics and Kinematic Models for Shock-Wave Propagation,” *Shock Waves*, Vol. 28, No. 2, 2018, pp. 401–416. https://doi.org/10.1007/s00193-017-0748-2

[66] Whitham, G. B., “Linear and Non Linear Waves,” John Wiley & Sons, Inc., 1974.

[67] Kamel, M. M., Ghoniem, A. F., Rashed, M. I., and Oppenheim, A. K. O. K., “Blast Waves in Real Gases,” *Acta Astronautica*, Vol. 4, No. 3, 1977, pp. 439–458.
https://doi.org/10.1016/0094-5765(77)90061-3

[68] Hairer, E., and Wanner, G., “Solving Ordinary Differential Equations II,” Springer Berlin Heidelberg, Berlin, Heidelberg, 1991. https://doi.org/10.1007/978-3-662-09947-6

[69] Meleshko, V. V., Gourjii, A. A., and Krasnopolskaya, T. S., “Vortex Rings: History and State of the Art,” *Journal of Mathematical Sciences*, Vol. 187, No. 6, 2012, pp. 772–808. https://doi.org/10.1007/s10958-012-1100-0

[70] Park, C., “On Convergence of Computation of Chemically Reacting Flows,” presented at the 53rd AIAA Aerospace Sciences Meeting, 1985.
https://doi.org/10.2514/6.1985-247

[71] Merk, H. J., “The Macroscopic Equations for Simultaneous Heat and Mass Transfer in Isotropic, Continuous and Closed Systems,” *Applied Scientific Research, Section A*, Vol. 8, No. 1, 1959, pp. 73–99. https://doi.org/10.1007/BF00411741

[72] McGee, H. A. jr, “Molecular Engineering,” McGraw-Hill, New York, N.Y., 1991.

[73] Raizer, Y. P., “Gas Discharge Physics,” Springer, 1991.

[74] Heijkers, S., Martini, L. M., Dilecce, G., Tosi, P., and Bogaerts, A., “Nanosecond Pulsed Discharge for CO2 Conversion: Kinetic Modeling To Elucidate the Chemistry and
Improve the Performance,” *The Journal of Physical Chemistry C*, Vol. 123, No. 19, 2019, pp. 12104–12116. https://doi.org/10.1021/acs.jpcc.9b01543

[75] Ombrello, T. M., Hammack, S. D., Carter, C. D., Busby, K., and Lefkowitz, J. K., “Scramjet Cavity Ignition Using Nanosecond-Pulsed High-Frequency Discharges,” *Combustion and Flame*, Vol. 262, 2024, p. 113335. https://doi.org/10.1016/j.combustflame.2024.113335

[76] Castela, M. L., Stepanyan, S., Fiorina, B., Coussement, A., Gicquel, O., Darabiha, N., and Laux, C., “A 3-D DNS and Experimental Study of the Effect of the Recirculating Flow Pattern inside a Reactive Kernel Produced by Nanosecond Plasma Discharges in a Methane-Air Mixture,” *Proceedings of the Combustion Institute*, Vol. 36, 2016. https://doi.org/10.1016/j.proci.2016.06.174

[77] Wang, J. M., Buchta, D. A., and Freund, J. B., “Hydrodynamic Ejection Caused by Laser-Induced Optical Breakdown,” *Journal of Fluid Mechanics*, Vol. 888, 2020, p. A16. https://doi.org/10.1017/jfm.2019.1066

[78] Lo, A., Cessou, A., Boubert, P., and Vervisch, P., “Space and Time Analysis of the Nanosecond Scale Discharges in Atmospheric Pressure Air: I. Gas Temperature and Vibrational Distribution Function of N₂ and O₂,” *Journal of Physics D: Applied Physics*, Vol. 47, No. 11, 2014, p. 115201. https://doi.org/10.1088/0022-3727/47/11/115201

[79] Blanchard, V., Minesi, N., Bechane, Y., Fiorina, B., and Laux, C. O., “Experimental and Numerical Characterization of a Lean Premixed Flame Stabilized by Nanosecond Discharges,” Vol. AIAA 2022-2255, San Diego, CA, United States, 2022. https://doi.org/10.2514/6.2022-2255

[80] Gablier, R., and Laux, C., “Hydrodynamic Recirculation Induced by Nanosecond Discharges: A Strategy to Mitigate Cumulative Heating,” In *AIAA SCITECH 2025 Forum*, Vol. AIAA 2019-0397, Orlando, FL, 2025. https://doi.org/10.2514/6.2025-0397

[81] Minesi, N., Blanchard, V., Pannier, E., Stancu, G. D., and Laux, C. O., “Plasma-Assisted Combustion with Nanosecond Discharges. Part I: Discharge Effects Characterization in the Burnt Gases of a Lean Flame,” *Plasma Sources Science and Technology*, 2022. https://doi.org/10.1088/1361-6595/ac5cd4

[82] Blanchard, V., “Thermochemical Effects of Nanosecond Plasma Discharges on Lean Flame Stabilization,” Université Paris-Saclay, PhD thesis, 2023.