Pulse-to-pulse coupling in cylindrical discharges

Several filamentary discharges can be applied to a combustible mixture, which can then ignite. The energy density of this discharge is a vital parameter, as it directly influences the local temperature rise and radical production. The goal of this article is to investigate how a previous discharge affects the energy density of a second discharge. To investigate the pulse-to-pulse coupling of filamentary discharges a one-dimensional numerical model is developed. In the developed model, the compressible Navier–Stokes equations are coupled to a plasma model. The plasma model is used to estimate the local energy density, while the compressible Navier–Stokes equations model the reactive flow. As a first step, skeletal air plasma chemistry is used, which includes fast gas heating, slow gas heating and the rapid generation of radicals. The skeletal plasma chemistry is combined with a detailed hydrogen combustion mechanism. Simulations in both air and hydrogen/air are conducted at several discharge energies and pressures. From the analysis of these results, we conclude that the main mechanism of pulse-to-pulse coupling is the reduction in molar density due to temperature rise.


Introduction
Plasma-assisted combustion, has been researched extensively in the past decades as a promising method for reliable ignition and flame stabilization in challenging conditions [1][2][3][4][5][6][7].Plasmas can be used to reduce the ignition delay by the rapid production of radicals.Non-equilibrium plasmas are considered the most interesting to activate combustion chemistry: a non-equilibrium plasma is characterized by highenergy electrons and a low gas temperature.The high-energy electrons are hypothesized to efficiently dissociate fuel and oxygen, producing initial combustion radicals [5,6], while not significantly heating the gas.As a result, the chemistry is activated directly by the electrons instead of activated thermally by plasma induced heating.
One of the key challenges in the numerical modeling of plasma-assisted combustion is the large separation in time and spatial scales involved.Nanosecond pulsed discharges are employed to allow for significant over-voltage beyond the breakdown voltage, resulting in a higher electric field.The higher electric field produces higher energy electrons, which is beneficial for chemistry activation.This inherently introduces short time scales (O(1 ns)), while ignition occurs at relatively slow timescales (O(1 ms)).Depending on conditions, like the species in the mixture, gas density and the magnitude of the electrical field, the discharge might be homogeneous or focused into thin plasma channels.Due to the large separation in time and spatial scales and the large chemical mechanisms involved fully resolving multiple discharges is not feasible.Due to their computational affordability, spatially homogeneous chemistry calculations are excellent for understanding chemical pathways, e.g.[8][9][10][11][12][13][14][15][16].However, comparing these results to experiments is challenging, not only because of uncertainties in the chemistry but also because of a nonhomogeneous plasma [17][18][19].As the chemical action performed by the plasma scales nearly linearly with its energy density [16], accurate estimation of the energy density is vital in ignition studies.
Several authors have performed 1D simulations with detailed chemical mechanisms to investigate the effect of a non-homogeneous plasma.Some of these authors have modeled plasma in a dielectric barrier setup, where the height between the electrodes is simulated [19][20][21][22], where the plasma is considered homogeneous in the direction parallel to the electrodes.While others have modeled the plasma in the radial direction and assumed a homogeneous plasma in the height and circumferential direction [23][24][25][26], more similar to thin plasma channels.
Due to the increasing computational power of highperformance computing clusters, it has recently become possible to perform two-dimensional simulations of plasmaassisted ignition.However, these simulations are still extremely time-consuming and thus limited in some way: Sharma et al [27] and Mao et al [28] studied ignition after only a single nanosecond discharge and have thus not investigated the impact of pulse-to-pulse coupling.When several discharges are investigated the flow and combustion chemistry are decoupled from the plasma simulation, e.g.Barleon et al [29,30] and Mao et al [31].Even with decoupling these simulations are still limited to ten(s) of discharges.If simulations are conducted of O (100) discharges, it is no longer feasible to simulate the plasma in detail and the plasma dimensions are prescribed, e.g.[32,33].All of these simulations are very insightful, but also computationally extremely expensive and require high-performance computing.
In this article, we will study the plasma in the radial direction similar to [25,26].We will extend upon this by considering nanosecond pulsed discharges.In these discharges, energy addition occurs significantly faster than pressure relaxation can take place.Therefore, the fluid can not be assumed to be incompressible.This implies that the gas density remains constant within a discharge and the instability demonstrated by Zhong et al [25] can not occur within a single discharge.Ideally an experimentally measured voltage profile is prescribed to activate the plasma.However, such a voltage profile can not (easily) take the shielding of the sheath and the voltage drop over dielectric material into account [20,34].This means that the actual voltage drop over the discharge gap is significantly lower than measured at a power supply.Moreover, the magnitude of the shielding depends on the breakdown of the plasma itself [34].As a result, the numerical deposited energy will be significantly overestimated and change from pulse to pulse if the experimental voltage profile is used.Alternatively a current profile could be prescribed as in Zong et al [25,26].However, we would like to have exact control of the deposited energy, as such we utilize a power deposition profile.The goal is to study the interaction of several gas discharges in air and hydrogen/air discharge, with a focus on understanding what mechanisms are important for pulse-to-pulse coupling.While not part of this article, in the future an improved understanding of pulse-to-pulse coupling can prove useful to develop improved prescribed plasma models as used in [32,33].Our model is only valid for the central region (not affected by the sheaths) of the plasma discharge, i.e. the filament or plasma channel.However, this model simplification does allow the problem to be described in 1D, which enables the use of implicit coupled solvers that are very efficient for stiff problems.As a result, our implementation is sufficiently computationally affordable to use semi-detailed chemistries and perform parameter sweeps (a simulation takes roughly one day on a single core).

Method
The numerical simulations in this article are inspired by the experimental setup of Patel et al [35][36][37], a schematic overview of their setup is shown in figure 1.In these experiments, repetitive nanosecond discharges between a spherical electrode and a flat electrode are investigated.The spherical top electrode is placed outside the quartz flow reactor and is embedded in a dielectric material, such that discharges only take place in the flow reactor.Our goal is to study the effect of pulse-to-pulse coupling on the discharge volume, the pink area in figure 1.In this study, we will consider the case in which the pulse repetition rate is sufficiently high compared to the flow timescale, such that the flow can be neglected.This allows us to assume that the discharge is homogeneous in the z-direction and is cylindrically symmetric.Thereby, the problem is reduced to a quasi 1-dimensional problem, where the rdirection is resolved and physical processes in the z-direction are modeled.As we consider nanosecond discharges in which gas heating occurs at a rate faster than the acoustic timescale, the compressible Navier-Stokes equations are solved.

Conservation equations
The conservation equations and constitutive relations for the multi-component compressible Navier-Stokes equations can be found in many textbooks, e.g.[38][39][40].Here, these conservation equations and constitutive relations are presented in  [35][36][37], the pink area is used to indicate where the discharges take place.cylindrical coordinates.The time dependent conservation of mass reads where ρ is the gas density, u is the gas velocity in the radial direction and r is the radial coordinate.In radial coordinates, the conservation of momentum is given by where τ is the viscous stress (constitutive relation provided in equation ( 15)) and p is the pressure.Here it is assumed that the electric field in the radial direction is small such that forcing due to ion wind can be neglected.The conservation of total enthalpy (where h t = chemical + sensible + kinetic is the total specific enthalpy) is given by: where q is the energy flux vector and Qp is the energy source due to the plasma.Here we would like to draw your attention to the fact that h t contains all relevant forms of energy.As such, the plasma source term Qp is all electrical energy transferred to the plasma, not only that due to Ohmic heating, i.e. sensible energy.The use of total enthalpy is also why there is no chemical heat-release source term in the conservation of enthalpy equation.The uτ term on the left-hand side is viscous dissipation.For the conditions of interest, the viscous dissipation is small, but it is included for completeness.Finally, N s equations are solved for the conservation of species, where Y i is the mass fraction of species i, U i is the drift + diffusion velocity of species i, M i is the molar mass of species i, and ωc,i and ωp,i are the source terms due to conventional and electron-induced reactions.
As charged species are included, the electric field ( ⃗ E) should be computed, as this will be a driving force in the diffusion of species.In our modeling, we assume that the electric field in the r and z-direction are independent.This assumption is valid if the applied electric field has a component only in the z-direction.As ⃗ E = −∇φ, this also implies that the electric potential (φ) can be decomposed in the r and z direction, i.e. φ(r, z, t) = φ r (r, t) + φ z (z, t) and thus E r = −∂φ r /∂r.The electric potential in the r direction is then given by the Poisson equation: where F is Faraday's constant, z i is the (signed) charge number of species i and ϵ 0 is the permittivity of free space.

Constitutive relations
To relate pressure, density and temperature (T) the ideal gas law is employed: where R is the universal gas constant and M is the mean molar mass ( X i M i , with X i the molar fraction of i).The molar density (N) of the gas is also computed from the ideal gas law: The enthalpy is related to the gas temperature via with where h 0,i is the formation enthalpy at T ref of species i and c p,i (T ′ ) is the temperature-dependent heat capacity of species i, which is evaluated using NASA polynomials [41].
The energy flux vector (in equation ( 3)) consists of three components: 1) The conductive heat transport, 2) the energy transport due to species diffusion and 3) the Dufour effect.In the radial direction, this is given by where λ is the mixture thermal conductivity and D T i is the thermal diffusion coefficient of species i. Species diffusion is modeled with the mixture averaged diffusion approximation, including drift due to the electric field and thermal diffusion.In radial coordinates the diffusive species flux reads (in equations ( 4) and (10)) where D i,m is the mixture averaged diffusion coefficient of species i and µ i is the electrical mobility of species i.The mixture averaged diffusion coefficients relate to the binary diffusion coefficients via: where D ij is the binary diffusion coefficient of species i in j, which is obtained using kinetic theory [42].The mixture thermal conductivity is computed from the species conductivities λ i using the mixing rule of Mathur et al [43] In the summation of the individual species' conductivities, the electrons are excluded.This is because the thermal conductivity of electrons is significantly larger than that of the gas.As a result, the mixture's thermal conductivity would be significantly increased when the ionization degree in the plasma becomes large.While electrons can transfer heat at a high rate, they can not, due to their small mass, transfer this energy effectively to atoms and molecules.The mobility (µ i ) is computed from the mixture averaged diffusion coefficient, using the Einstein relation: Strictly, the above equation is not valid for the electrons during the electrical discharge as the electron temperature is significantly larger than the gas temperature.However, the duration of the discharge is only short and the electric field in the radial direction is expected to be tiny.More importantly, the electric field is expected to be tiny because the externally applied electrical field only acts in the z-direction, which means that only the diffusion of ions and electrons can induce an electric field in the r-direction.In case an electric field is induced in the r-direction, the lightweight, charged and very mobile electrons will be forced (due to the electric field) to rapidly diffuse and reestablish a (close to) neutral plasma.
In other words, as long as the combination of µ e and D e,m is such that the plasma remains (close to) neutral the error induced by using the Einstein relation with the gas temperature is small.While not shown, in post-processing we validated the use of equation (14) for the electrons by checking if the plasma remains charge neutral.The forcing of species due to pressure gradients is neglected in equation (11), which can be large if X i ̸ = Y i , i.e. it can be large for the electrons.Similarly to the approximation of µ i , the duration of strong pressure gradients is short and electrons are tied to the ions.Assuming that only ∂u/∂r is nonzero the viscous stress can be computed from where µ is the dynamic viscosity, the derivation of this equation is provided in appendix A. The additional 2u/(3r) term arises due to the radial expansion of the domain with increasing r.As a first step, we have neglected the bulk (or second) viscosity in the stress tensor, which for the conditions of interest might be large.The magnitude of the bulk viscosity is condition dependent, and it is related to thermal nonequilibrium [44,45], i.e. vibrational, rotational and or electronic temperatures are not identical.It contributes to the viscous stress when heating occurs at a rate faster than the thermal relaxation timescale of the gas.For our conditions, the gas is heated several hundred degrees in less than 10 ns, which is much faster than the vibrational relaxation timescale [46][47][48][49].This implies that the bulk viscosity is expected to be non-zero for our conditions, due to vibrational non-equilibrium induced by the fast heating.As the bulk viscosity is related to thermal non-equilibrium of individual species its physical effect can be included in two ways: 1) the bulk viscosity is modeled and added to equation (15) or 2) the non-equilibrium of individual species is resolved by including internal states of molecules and atoms as species.In this study vibrational levels of nitrogen are included in our chemistry, meaning that vibrational non-equilibrium of N 2 is captured by the chemistry.As such, the effect of the bulk viscosity due to the N 2 vibrational nonequilibrium is included in the model, although not via the bulk viscosity.The vibrational levels of oxygen and the fuel are not included, but the concentration of oxygen is significantly less than that of nitrogen and the relaxation rate is faster [49], meaning the contribution of these on the bulk viscosity is significantly smaller.

Plasma modeling
The model for the z-direction of our problem is inspired by the studies of Popov [23,24] and Zhong et al [25,26].The goal of this model is to provide spatial-temporal energy ( Qp , in equation ( 3)) and species ( ωp,i , in equation ( 4)) source terms.In this model, these source terms are dependent on the local gas composition and gas temperature.In the models of [23][24][25][26], a user-defined plasma current is required.As a result, the deposited energy can not be directly controlled.To have direct control over the energy deposition, our implementation is slightly modified, such that a temporal power deposition function P(t) can be used.
Assuming that the plasma behaves like an Ohmic resistor, the power deposition relates to electrode voltage (φ elec (t)) and current (I(t)) via The current due to the electron flux in the discharge gap can be computed from Note the use of the superscript I in the electrical mobility, as the Einstein relation (equation ( 14)) is not used here.Because In the center the electric potential crosses the discharge gap, which has a distance dgap.Outside the center line, the electric potential also encounters the dielectric material surrounding the electrode lengthening its path with d dielec .
the electric field in the z-direction is large, the Einstein relation is inaccurate.Strictly, the electron mobility then depends on the reduced electric field (E/N).However, in our modeling we assume that Fz e µ I e N is constant.As a first step, this is a reasonable approximation as Bolsig [50] calculations show that µ I e N is nearly constant for a sufficiently large reduced electric field (E/N > 10 Td).
The electric field in the z direction is approximated as where l(r) is the effective path length of the electric field lines.φ elec (t) will evolve temporarily but not spatial and l(r) evolves spatially but not temporarily.The effective path length is introduced to introduce a spatial dependence of the electric field at t = 0, such that the plasma is most intense in the center of the domain.l(r) is based on the physical model illustrated in figure 2. As illustrated, it is assumed that the applied electric field inside the discharge only has a component in the z-direction.Moreover, the electric field attaches tangent to both electrodes.Near the top electrode, a dielectric material is present (with relative permittivity ϵ r ), which avoids discharge from occurring outside the discharge gap.The effective path length of the electric field can then be reconstructed as two straight lines, with where d g is the height of the discharge gap and R is the radius of the electrode.
Combining equations ( 16)-( 18) allows us to derive an analytical expression for φ elec : where the previously made approximation (Fz e µ I e N is constant) is used.Placing these variables, as a constant, outside the integral allows us to provide the analytical expression of equation (21).
Knowing the voltage, the current density (J(r, t)) can be computed as which can then be used to compute the local plasma power (needed in equation ( 3)) Finally, plasma reactions are introduced to the model to include the chemical effect of the plasma.The method used is similar to the approach by Castela et al [51] (and derived works [30,52,53]), but it is extended to include the influence of E/N.In our implementation, these plasma reactions do not directly depend on the electron concentration.Instead, they depend on the plasma power Qp , i.e. the reaction rate is defined as an energy branching.For example, a plasma reaction in this framework might read whose rate is determined by the rate coefficient k p,j .This rate coefficient is computed from an Arrhenius-like expression in E/N: where A j , n j and [E/N] a,j are the (Arrhenius) constants of reaction j.The Arrhenius-like expression is used to model the dependency of the energy branching on the reduced electrical field, i.e. mean electron energy.Energy branching fractions are obtained if equation ( 24) is multiplied by the enthalpy change of the reaction and the mole fraction of the colliding partner.
The rate of progress of reaction j and the plasma source term of species i (ω p,i ) are then computed per usual.While (R.1) does, intentionally, not depend on the electron concentration the rate of progress does via Qp .As Qp depends on the local concentration of e, the reaction rate itself does depend on the electron concentration (see equations ( 22) nad ( 23)).Formulating the rate coefficient of plasma reactions like equation (24) guarantees consistency between the deposited energy and the change in formation enthalpy of the mixture.

Numerical methods
The resulting set of equations is challenging to solve and require dedicated solvers.We have extended our in-house low-mach flame solver (Chem1d, e.g.[54][55][56]), to the entire mach regime.To achieve accurate solutions while maintaining stability, the MUSCLE method is employed for finite volume interpolation, with the van Albada limiter [57].To ensure consistent face properties, density, velocity, temperature and species mass fractions are interpolated to the face.From these properties the mixture enthalpy and pressure are reconstructed at the face.Using these face interpolated values, numerical convective fluxes at the cell faces are estimated using the SLAU2 flux-splitting method [58,59], designed for all-speed flows.Application of these methods provided sufficient stability such that no preconditioning of the conservation equations is utilized.At r = 0, a symmetry boundary condition is implemented by forcing the velocity to zero and all other variables have zero gradients.On the right boundary r = R max the NSCBC procedure [60] is employed to provide a well-possed outflow boundary.The effect of chemistry is taken into account at the boundary [61], while the effect of diffusive transport is neglected.No corrections have been applied to the NSCBC procedure for the fact that the domain is in cylindrical coordinates.From our testing, we concluded that this was not required as more than 99.9% of the acoustic energy is absorbed by the boundary.

Model limitation
The plasma is assumed to be uniform in the z-direction.However, near the bottom and top electrode this assumption is unlikely to hold.Close to the electrodes a plasma-sheath is formed, which can have a large reduced electric field.As a result, a significant portion of the energy is deposited near the electrodes [20].As this is energy is deposited close to the wall a significant portion of it is directly lost (via heat loss to the wall), such that the largest temperature rise and radical production is halfway between the electrodes [21].The region at some distance of the electrodes is expected to be uniform in the z-direction and-due to the lack of heat loss to the wall-is expected to be most relevant for ignition.This uniform central region of the plasma is exactly the region that we intend to model in this article.As a result of our modeling assumption, the temporal power deposition P(t) is the energy deposited exclusively in the filament, not the entire plasma.For a DBD-discharge, where the deposited energy remains rather constant from pulse-to-pulse [20,62], this implies that the fraction of energy going into the filament is assumed constant in our model.

Chemistry
In this article plasma discharges in air and hydrogen/air are studied.The combustion chemistry is modeled using the detailed hydrogen mechanism from NUIGMech 1.1 [63].This hydrogen combustion mechanism is extended with the ozone chemistry from HP-mech [64].The air plasma mechanism is the fast gas heating model of Popov et al [7,65] and Mintoussov et al [66], from which we have constructed a skeletal mechanism.The reduction of the detailed mechanism to the skeletal mechanism is explained in the remainder of this section.The final chemical model includes fuel oxidation, ion chemistry, vibrational non-equilibrium and electronically excited chemistry.Thermodynamic properties of electronic and vibrationally excited species are obtained with the consistent method outlined in Hazenberg et al [67].As the focus of this article lies on the impact of chemical reactions on the pulse-to-pulse coupling of nanosecond discharges plasma discharges, it is assumed that hydrogen does not react due the plasma itself.However, it does react due to the oxidation chemistry in NUIGMech 1.1 [63].This implies that the plasma chemistry outlined below is mainly designed for air discharges.This assumption is more reasonable than it seems at first: In our model, this reduced electron mobility is the only parameter that links the electron scattering processes to the local energy deposition (reaction rates are derived from this).We have performed Bolsig calculations with H 2 /O 2 /N 2 mixtures, from which we conclude that the reduced electron mobility does not change significantly with the addition of hydrogen.Moreover, the electronically excited state chemistry of H 2 is not extensive, electronically excited hydrogen either relaxes to the ground state (generating heat) or dissociates, see the reaction mechanisms of, e.g.[14,15,68].The dissociative effect of electrons on hydrogen is neglected, but this is only a relatively small fraction of the energy (based on Bolsig calculations, less than 5%), while the relaxation to ground state is implicitly included as (in our model) all energy not spend on chemical reaction will be transferred as sensible energy.The chemical mechanism is made available through the supplementary materials.

Skeletal ion chemistry
Repetitive nanosecond discharges are studied in this article.Using the plasma model outlined in section 2.3 the spatial distribution of the energy deposition is estimated.The spatial energy loading depends strongly on the local electron concentration from before the discharge.Therefore, the main goal of our ion chemistry is to predict the spatial distribution of electrons.At the same time, we would like to achieve this with the smallest possible set of species and reactions to reduce the computational load.To this end, the following ions are included N 2 This choice is rationalized as follows; the N 2 + and O 2 + ions are generated directly by the plasma via where one electron has been placed between brackets as it is, numerically, implicitly included (see section 2.3).At sufficiently high pressure (P ∼ O(1 bar)) [65,69], both of these ions are quickly converted to N 4 + and O 4 + via the third body reactions: Reactions and rates involved with N 2 ( 3 ), k is in cm 3 (s mol) −1 and Ea in cal mol −1 .The reaction rates are based on Popov et al [65].
1.300 × 10 13 0.000 0.0 1.300 × 10 13 0.000 0.0 Dissociative attachment (R.6) requires the O 2 bond to be broken and thus requires a sufficiently large temperature.Moreover, as near atmospheric pressure is considered here, the third body process (R.7) is expected to be fast.The O − 3 ionincluded in Popov et al [7,65]-is not included in our mechanism, because the formation of O − 3 requires O 3 .While at room temperature a plasma can produce a significant amount of O 3 , at elevated temperatures the produced O 3 is expected to quickly dissociate.As we are considering a plasma in which a significant temperature rise will occur, O 3 can not accumulate over multiple discharges.As a result, the rate of progress of is expected to be low.Therefore, only the O − 2 ion is considered here.Compared to the work of Popov et al [7,65], the ion chemistry is extended with the important detachment reactions by O 2 and N 2 [70]: where the rate is obtained from [48].As a first step, it is assumed that H 2 , H 2 O or any combustion intermediates do not affect the ion chemistry.

Skeletal electronically excited chemistry
Electronically excited states are often considered to be important in plasma-assisted combustion.However, as a first step-and to significantly reduce the computational resources required-we have excluded all electronically excited states except a combined N 2 triplet state (N 2 ( 3 )).This combined triplet state consists of the following electronic states: The rationale behind this reduction is as follows: First, all triplet states quench, at ambient pressure, faster than typical combustion timescales (O(100 ns) [5,51,65]).Second, A 3 Σ + u quenches slowly against ground state N 2 (X 1 Σ + g ) [49,71,72].Third, when Fourth, all triplet states quench fast against O 2 and H 2 , which dissociates O 2 and H 2 [4,7,65].Fifth, A 3 Σ + u does quench rapidly against itself, this produces one N 2 (X 1 Σ + g ) and [49,65,73,74].To summarize, the N 2 ( 3 ) states are strongly coupled, and their main chemical effect is the dissociation of fuel and oxidizer.As this chemical effect occurs way faster than typical combustion timescales, it is reasonable to simplify the entire N 2 ( 3 ) chemistry to a single species.The resulting triplet chemistry is summarized in table 1, the branching over N 2 (v = x) will be explained in the next subsection.Electronically excited O 2 , is not included as the singlet states (O 2 (a 1 ∆ g ), O 2 (b 1 Σ + g )) have been shown to have a small effect on the oxidation of hydrogen [5,75].

Vibrational non-equilibrium
It is assumed that O 2 remains in vibrational equilibrium, while the lowest 9 vibrational levels of N 2 are included in the reaction mechanism.As the N 2 vibrational energy levels are an important energy reservoir in air discharges, while O 2 vibrational levels relax quickly [11,49,51,65].The thermodynamic data of the lowest 9 vibrational levels of N 2 are obtained with the method outlined in Hazenberg et al [67].The nitrogen vibrational kinetics include V − V and V − T processes.For both the V − V processes and the V − T processes the scaling relation given by Guerra et al [49] is used.V − T relaxation by O, N 2 , O 2 and H 2 O is taken into account, where the O 2 relaxation rate is obtained from the N 2 rate and the conversion provided by Guerra et al [49].The rates of Table 2. Rate parameters of electron processes, the computation of the rate coefficient is given in equation (24).The units of A j is mol/J, n j is dimensionless and [E/N] a,j is in Td.Note reactions starting at N 2 (v = 0) are included for every vibrational level of N 2 , see section 3.4.

Reaction
) is based on the fast gas heating work of Popov et al [24,66].We have matched the heat release of Popov [24] by distributing the chemical energy over two vibrational levels.

Electron chemistry
The chemistry activation by the electrons is significantly reduced.It is assumed that the rate constants are not dependent on the gas composition.Note, the reaction rate itself still depends on the species concentration and the power density, see section 2.3.These rate constants are obtained from Bolsig [50] calculation using the IST-Lisbon [76] crosssections for O 2 and N 2 .In the Bolsig calculations, all processes are included to get the correct energy transfer.The reaction rates are then obtained by normalizing the rate by this energy transfer.This implies that excluded (or reduced) processes are assumed to instantaneously relax to thermodynamic equilibrium.An overview of all included processes is provided in table 2. Some simulations are conducted in dry air.In these simulations, when the gas temperature is low, the vibrationaltranslation relaxation is longer than the time between discharges.As such, significant vibrational non-equilibrium can accumulate over multiple pulses.To account for this exciton from v > 0 is included in the chemistry.For vibrational excitation, it is assumed that the rates of and are equal.This is reasonable if the energy of vibrational level i, j and i + j can be approximated by a harmonic oscillator.For excitation to an electronic state or ionization the rates of and are also equal, which relies on the assumption that the threshold energy of excitation and ionization is much larger than the energy of vibrational excitation.

Results
The model that is outlined above requires several parameters.These parameters are based on the experimental setup of Patel et al [36,37].An overview of the key parameters is provided in table 3. The discharge gap and radius of the top electrode sphere are directly taken from the setup dimensions.The relative permittivity is set to 5, which is close to the value of quartz.This value is chosen as it provides a small fall-off for the reduced electric field outside the discharge center.This falloff ensures that the first discharge occurs in the center of the domain.After the first discharge, the temperature and residual electrons ensure that consecutive discharges also occur in the center.The reduced electron mobility (Fz e µ I e N) is based on Bolsig [50] calculations in dry air at a representative reduced electric field.It is assumed that the power deposition (P(t)) is approximated by a Gaussian, which has a width of σ p = 2 ns.In all our simulations the pulse repetition rate is set to 3 kHz.Finally, simulations are conducted at three energy depositions ∆E p = 437.5,875, and 1750 µJ and two pressures P 0 = 500 and 750 mbar.The deposited energy is related to P(t) via with where t 0 is the start time of the discharge (t 0 = (N p − 1)/f pulse , with N p the pulse number and f p the pulse frequency), t off = ln(1000)σ p is an offset with respect to t 0 to ensure the full  energy is deposited in the first pulse, t end the end of the discharge and C is a scaling constant of the Gaussian, which can be computed from ∆E p .The energy depositions are a quarter and a half of the energy deposition in the experiments of Patel et al [36], while the pressures are in the same range.The simulations are conducted on a domain of r = 0 to R max = 1 cm, with 2000 finite volume cells.R max is sufficiently large such that the entire discharge occurs within the domain.Moreover, it is sufficiently far away, such that the boundary curvature is not too large to affect the performance of the Navier Stokes characteristic boundary conditions.The 2000 grid cells results in a cell size of 5 µm, which is sufficiently small to resolve the plasma.Finally, in all simulations it is assumed that the initial mole fraction of electrons has a small positive value (X e,0 = 1 × 10 −14 ), to avoid a division by zero in equation (21).

Single discharge in air
We will first discuss the results of a single discharge in air.This discharge has an energy deposition of 875.0 µJ, where the gas initially is at 300 K and 500 mbar.In figure 3 the temporal evolution of the reduced electric field and the (prescribed) power deposition is shown.The initial reduced electric field is zero as the power deposition is zero.Within the first 5 ns, the reduced electric field quickly rises, shows a rapid drop, and then decreases nearly linearly.The short peak in the reduced electric field is the result of our modeling assumptions, let us elaborate: As the electron concentration is very low (X e = O(1 × 10 −12 ), see figure 4), the gas conductivity is extremely low.To deposit energy with such a low electron concentration the applied voltage (computed with equation ( 21), note the division by the electron molar fraction) must be very high.As a result, the reduced electric field quickly increases to a peak value of approximately 300 Td at r = 0.0 mm.Such a high electric field is sufficient to rapidly accumulate electrons and the voltage quickly drops to a lower and linearly decreasing value.The short peak in the reduced electric field is not expected to impact the results as the deposited energy is extremely low between 0 and 5 ns, see right frame of figure 3 for reference.Outside the domain center the reduced electric field does not increase to that extent, which is due to the artificial lengthening of the electric field lines, see equation (19).
The increased reduced electric field results in the increase of electrons in the discharge gap.As the rate of electron production increases with increasing electric field, most electrons will be produced in the center.This can also be seen in figure 4, where the electron mole fraction at r = 0.0 and 2.5 mm is shown.The fraction at 5 mm is not included, as the reduced electric field at that position is insufficient to significantly increase the electron concentration.Even when the reduced  electric field is reduced by only O(50 Td) the difference is substantial, while the electron mole fraction increases to nearly 1 × 10 −4 at 15 ns in the center, it remains below 1 × 10 −10 at 2.5 mm.As a result of this fall-off in the reduced electric field, the discharge will deposit most of the energy in the center of the domain.
In figure 5 the spatial distribution of the time-integrated power deposition is shown for three different cases.The first case is the same as above (875.0µJ and 500 mbar), the second is at half the deposited energy (437.5 µJ and 500 mbar) and the third at increased pressure (875.0 µJ and 750 mbar).For the first case (875.0 µJ and 500 mbar), the peak energy density is 138 kJ m −3 and it falls off to very low values within 1.5 mm.This is to be expected based on the results of figure 4, as the electron concentration is only marginally increased at 2.5 mm.When the deposited energy is halved to 437.5 µJ (case two) the entire curve is reduced in height.Upon careful inspection, it can be seen that the peak value is not halved (it decreases from 138 kJ m −3 to 75 kJ m −3 ), which implies that the width of the discharge at lower energy deposition is slightly smaller.For the case where the pressure is increased to 750 mbar (while the energy deposition is maintained at 875.0 µJ), the peak value of the energy deposition is very slightly increased.As a result, the energy loading per molecule reduces with pressure in these simulations.Finally, the shape of all these deposited energy curves is comparable to a Gaussian distribution.
In our chemical model the energy of the discharge can activate several chemical processes: 1) Vibrational excitation of N 2 , 2) Electronic excitation of N 2 , 3) Dissociation of O 2 and 4) ionization of oxygen and nitrogen.Of these, ionization is a minor contribution to the total energy transfer in the plasma (< 5%), which is to be expected at these reduced electric field e.g.[3,20,65].The fraction of energy in mode 1 to 3 is computed and shown in figure 6.In the left frame, the spatial distribution of these energy fractions is shown at t ≈ 13.8 ns.It can be seen that the energy branching over the different chemical modes is practically independent of radius.The energy branching does fall off to pure sensible energy (all fractions go to zero) above r = 1.5 mm.But in this region, almost no energy is deposited, see figure 5.The fraction of energy in N 2 ( 3 ) does start decreasing at a smaller radius than the energy in N 2 (v) or O.This can be explained from the larger activation energy for N 2 ( 3 ) (see table 2), which requires a larger reduced electric field to overcome.When comparing the solid line (discharge at 875 µJ) to the dashed line (discharge at 437.5 µJ) it can be seen that the distribution over the different energy modes is rather independent of discharge energy.
In the right frame of figure 6, the temporal evolution of the energy branching is shown.Initially, all energy is in sensible form (all fractions are zero).When the discharge starts, all fractions immediately increase.Roughly at the moment of peak energy deposition (t = 13.8 ns), the fraction of energy in N 2 ( 3 ) reaches its maximum at approximately 25%.At this instance, a plateau in O 2 dissociation is observed at approximately 15%.In these conditions, the N 2 ( 3 ) state has a lifetime (in this model controlled by the decay of N 2 (A 3 Σ + u )) in the order of 100 ns [5,51,65].This can also be seen from the fraction of energy in the N 2 ( 3 ) chemical mode, which starts decreasing directly after the peak and becomes zero within 500 ns.Between 10 ns and 20 ns, a fraction of the energy is used for vibrational excitation via This process is extremely fast (it occurs within the discharge time) and thus limits the peak value of N 2 ( 3 ).It is also responsible for the rapid decline of N 2 ( 3 ) between 13 and 20 ns.Note that reaction (R.15) is second order in N 2 ( 3 ), meaning it is very sensitive to the concentration of N 2 ( 3 ).
After 20 ns the reaction (R.15) is no longer the dominant mechanism for removal of N 2 ( 3 ).The dominant mechanism is then dissociative quenching against O 2 forming O. Roughly half of the remaining energy in N 2 ( 3 ) energy is spent on this mechanism, as the fraction of energy in O increases from 15 to 25% during this time.The other half of the N 2 ( 3 ) energy is not used in a chemical mode and is thus spent on gas heating.After 1 µs 30% of the energy is in vibrational excitation, 25% in O 2 dissociation, and thus 45% of the energy in sensible form.Compared to the empirical model of Castela et al [51], these numbers are rather different: In Castela et al 45% is spent on vibrational excitation and 35% on O 2 dissociation.The main reason for the difference is due to the skeleton chemistry of N 2 ( 3 ).In our model quenching of N 2 ( 3 ) by O 2 (which induces O 2 dissociation) does not result in vibrational excitation of N 2 , which in practice might be the case [7,65].It should be noted that the energy branching in an air plasma is still up to debate; in the more recent work of Barleon et al [30], only 20% is spent on vibrational excitation.
In figure 7 the centerline temperature (T c ) is shown for the same three cases as in figure 5.During the discharge (peak power deposition at 13.8 ns), the temperature quickly rises by at least 50 K in all three cases.The temperature rise during the plasma is largest for the 875.0 µJ and 500 mbar case and smallest for the 437.5 µJ and 500 mbar, which matches expectations based on the deposited energy (figure 5).In all cases, the temperature remains increasing after the discharge till about 300 ns.This increase is due to the quenching of N 2 ( 3 ) (see also figure 6), which is slightly faster when pressure is increased.Based on the energy fractions of figure 6, this second increase in temperature is expected to be roughly half of the direct plasma heating, as ±50% of the N 2 ( 3 ) energy is spent on O 2 dissociation.
Between 0.3 and 2 µs, the temperature is reduced, which is due to the pressure relaxing back to ambient.The discharge itself takes place in the nanosecond timescale, as a result, the energy is deposited in an effectively constant volume gas.As a result, the pressure has increased by the same ratio as that of temperature.The pressure field then relaxes within the acoustic timescale, in this case approximately 2 µs.This relaxation produces a pressure wave.The pressure wave leaving the discharge area has both a peak and a valley.As such, the pressure in the center initially drops below the ambient pressure before restoring to the ambient value.This is due to the cylindrical domain (gas expands) and the outward moment of the peak pulling air with it, creating a reduced pressure weak.For both the 875 and 437.5 µJ at 500 mbar case the pressure wave is shown in figure 8, where the acoustic pressure is shown P ′ = P − P 0 .As a result of this pressure wave shape, T c first undershoots before stabilizing after 3 µs, see figure 7.After approximately 20 µs the temperature starts to increase again.This increase in temperature is due to the formation of ozone (O 3 ) via which is an exothermic reaction.Finally, when comparing the peak pressure of the 875 and 437.5 µJ case, it can be noted that the peak is reduced by nearly a factor of three.This is interesting to note, as the energy density in the center is a bit more than half when comparing the 875 and 437.5 µJ cases, see figure 5.Moreover, the fraction of energy in each of the chemical modes very similar between the two cases (see figure 6), meaning one would except that the temperature rise and thereby pressure rise differ by a factor of two.However, the relaxation of the N 2 ( 3 ) occurs faster than the pressure relaxation, but the time scales still partially overlaps.In other words, the relaxation of N 2 ( 3 ) does neither occur at constant pressure nor constant volume, i.e. compressibility should be taken into account for the N 2 ( 3 ) relaxation.

Single discharge in hydrogen-air
Simulations are also conducted in a hydrogen/air mixture with a fuel-equivalence ratio (ϕ) of 0.7.The fraction of energy in the different chemical modes for the 875 µJ at 500 mbar and the 437.5 µJ at 500 mbar cases are shown in figure 9.In this graph, the fraction of energy in H 2 dissociation instead of the fraction of energy in N 2 ( 3 ) is shown, as the shape of N 2 ( 3 ) is very similar to figure 6.Moreover, until 1 µs the shape of O 2 dissociation and N 2 (v) is also very similar to that in dry air  (figure 6).However, the absolute fraction of energy in O 2 dissociation and N 2 (v) is significantly reduced compared to dry air.This can easily be explained based on the lower mole fractions of O 2 and N 2 , which implies that production rates from electron collisions are reduced.When comparing the 875 µJ (solid) and the 437.5 µJ (dashed) cases, it can be noted that the fraction of O is slightly higher between 1 and 100 µs.This is due the temperature dependency (centerline temperature is lower for the 473.5 µJ case) of the rate at which O attacks H 2 , as a result the peak value of H is also much lower.After 10 µs, large differences in the fraction of energy in O 2 dissociation are present when comparing figures 6 and 9.When hydrogen is present O starts attacking it via The increase in H radical concentration can also be seen in figure 9, where the fraction of energy in the H 2 dissociation increases slightly between 10 and 100 µs.As the temperature is still relatively low (see figure 10) the oxidation capabilities of this mixture are limited [16].At these temperatures, hydroperoxyl (HO 2 ) is an important intermediate in the oxidation of H 2 , which is mainly formed via the third body reaction: Figure 10.centerline temperature for an H 2 /air mixture at ϕ = 0.7 for three different conditions: 1) 875 µJ at 500 mbar, 2) 437.5 µJ at 500 mbar and 3) 875 µJ at 750 mbar.
When the temperature is below 450 K the H, OH and HO 2 radicals mainly accelerate the recombination of O into O 2 [16].As a result, the chemical heat release due to this first discharge is very limited.In the centerline temperature graph (figure 10) this can be observed: In the 875.0 µJ and 500 mbar case, the temperature rises an additional O(30 K) due to the formation of water vapor between 10 and 100 µs.This can be compared to a temperature rise of O(10 K) in air, see figure 7.In the reduced energy case (437.5 µJ) and the increased pressure case (750 mbar), the energy density is lower.As a result, the temperature in the center of the discharge is also lower, making the hydrogen oxidation even less efficient.

Multiple discharges
In this study, the plasma is operated in pulsed mode at a frequency of 3 KHz.After the first discharge, ion chemistry quickly removes free electrons from the discharge gap.For the case of 875 µJ and 500 mbar in air, the temporal evolution of all ions at the centerline is shown in figure 11.In the left graph, the negatively charged species are shown.The electrons produced by the discharge attach against the abundantly available oxygen and form O − 2 within a couple of nanoseconds.The mole fraction of O − 2 peaks roughly 40 ns after the discharge and then decays.After some time a balance is formed between electron attachment, detachment and recombination resulting in the first-order decay of charged species after 2 µs.In the right graph of figure 11 + .This series of ion conversion tions are very typical for N 2 /O 2 type discharge [77].
The ion chemistry is rather sensitive to the temperature and local gas composition.As a result, the ions do not decay everywhere at the same rate.In figure 12 the temperature and ion distribution is shown just before the second discharge (t = 1/3 kHz ≈ 0.333 ms).Near r = 0 mm, the electrons and ions all have mole fractions in the order of 1 × 10 −10 .The mole fraction of electrons has a nearly flat profile up to 1.5 mm and then falls off quickly.When comparing this to the temperature shown in the right graph, it can be concluded that this corresponds to the edge of the discharge.Everywhere in the domain, the electrons quickly attach to the available oxygen, which results in a higher mole fraction of O − 2 than e.A change in the balance between O 2 + and O 4 + due to temperature can also be seen.In the center, the mole fraction of O 2 + increases and becomes almost equal to that of O 4 + , while at the edge the O 4 + is very much dominant.Outside r = 0.5 mm, the gas remains charge neutral by balancing the concentration of O − 2 and O 4 + .While not directly shown, the gas has almost exactly zero net charge everywhere and at every time in the domain.Finally, the temperature profile approximates a Gaussian with a maximum of 416 K in the center.
In figure 13 the spatial distribution of ions and temperature are shown for H 2 /air at ϕ = 0.7 for the same discharge settings as in figure 12.When H 2 is added to air, a notable difference can be observed: The mole fraction of e is significantly reduced, while that of O − 2 (and positive ions) is increased.This can be explained as fellows: In the mixture containing H 2 , the lifetime of O is reduced (compare figures 6 and 9  9).The free electrons keep recombining with the positive ions, reducing their concentration, while the O − 2 ion concentration remain relatively stable.As a result, the electron mole fraction before the second discharge is expected to be lower compared to dry air when H 2 is added.Finally, when comparing the temperature plots in figures 12 and 13, it can be noted that the width of the profiles is very similar.However, the higher temperature in the H 2 /air case compared to the air case is due to the chemical heat release which occurs when H 2 is present.
The spatial energy distribution of the first and second pulse for both ϕ = 0.0 and ϕ = 0.7 are shown in figure 14.Not easily read from the graph, but the first discharge in both air and air/H 2 are practically identical.In both the ϕ = 0.0 and ϕ = 0.7 cases, the second discharge has a significantly smaller radius.As a result, the energy deposition in the center is significantly increased.For the air case, the peak value is increased from 138 to 812 kJ m −3 , while for ϕ = 0.7, it is increased from 138 to 935 kJ m −3 .Based on the residual electrons from the previous discharge (see figures 12 and 13) one might expect that the second discharge in air/H 2 is wider.Because the electron concentration is lower in the center than at the sides compared to the air discharge.However, the temperature for the air/H 2 case is higher in the center, which increases the reduced electric field in the center compared to the edges.
The reduced electric field at peak energy deposition (13.8 ns after the start of the pulse), is shown in figure 15.In the left frame, the result for air is shown and in the right frame, the ϕ = 0.7 case is shown.For the first discharge, the spatial reduced electric field is completely determined by the electric field effective path length (equation (19)).Due to this effective path length, the reduced electric field reduces away from center of the domain.The reduced electric field is higher for ϕ = 0.7 than in the ϕ = 0.0 case, as the ionization of hydrogen is not included in the model.Note, the ionization of hydrogen requires higher electron energy (and thus higher E/N) than the ionization of oxygen, as such one would expect this effect to be present even if the ionization of hydrogen was included in   the model.In the second discharge, the temperature field (see figures 12 and 13) results in a lower number density in the center of the domain.The resulting reduced electric field is now a combination of the effective path length and the temperature field, which results in an enhancement in the center of the domain.This increases reduced electric field in the center of the domain, which results in an increased ionization rate.The increased ionization rate then increases the electron concentration.Finally, the increased electron concentration increases the energy deposition, resulting in the observed contraction of the discharge.

Pulse-to-pulse coupling
To further analyze the pulse-to-pulse coupling a reduced order model is introduced.During a nanosecond discharge, the production rate of electrons is so fast that transport can be neglected.As such, the discharge is effectively a 0D process in which cells are only coupled via the electrode voltage computation (equation ( 21)).Moreover, for an exponential growth of the electrons to occur the production of electrons must be much faster than attachment or recombination reactions.When we additionally assume that changes in the gas composition and density are negligible, a reduced model can be obtained, see appendix B: where k eff,prod is an effective electron production rate coefficient.Note that the N term remains constant even if large amounts of energy are deposited as the energy is added sufficiently fast that it can be considered constant volume heat addition.Under the assumption of constant electrode voltage, the solution to this ODE reads The electron mole fraction increases exponentially, which is only valid if the production rate is much larger than any loss process.The reduced electric field at which this occurs is the breakdown field, which in air is approximately 120 Td [70].In appendix B, we also provide the solution when loss processes are considered and show that this is equivalent to equation (28) as long as X e has not approached its quasi-steady state value.The effective electron production rate is given by where N p,e are the number of plasma reactions involving electron production, the rate coefficients are those of equation (24) and v j,i is the stoichiometric coefficient of species i in electronproducing reaction j.As indicted by example reaction (R.1) v j,e = 0, as the electron concentration is in the Qp term.In other words, this effective rate is the sum of the electron-producing reactions per equation (24), such that it only needs to be multiplied by Qp to obtain the electron production rate.In figure 16 this effective electron production rate k eff,prod in O 2 , N 2 and air as a function of E/N are shown.From this figure it can be seen that the effective electron production rate becomes nearly linear above 150 Td, which is the case for the center of the domain, see figure 15.This implies that the electron mole fraction after a pulse relates like ln when the electric field is sufficiently large.As a result, the increase of electron concentration during the discharge is extremely sensitive to the local reduced electric field.
Using the analytical solution of the electron mole fraction (equation ( 28)), it can be shown that the local energy density (∆e p = ´Q p dt) scales as or using the scaling relation of equation ( 30) The number density depends on temperature via the ideal gas law and that pressure has relaxed before the second discharge.As such, the local energy density is expected to scale as This implies that the temperature increase from a previous discharge will significantly intensify the discharge in that area during the next discharge.
Based on the assumptions required to arrive at equation (28), it is clear that the final electron concentration depends linearly on the initial electron concentration.This means that accurate prediction of the spatial distribution of electrons, due to recombination or transport between discharges, does impact the discharge size in consecutive discharges.However, the influence of the reduced electric field on the final electron concentration is much larger.Typically, X e,end /X e,0 = O(1 × 10 6 ) (see figure 4), such that exp((E/N) 2 ) (and thus exp((ET) 2 )) is not small and is thus very sensitive to small changes in E/N.At the same time, we should realize that the uncertainty in ion chemistry is rather large, meaning that the spatially accurate prediction of electrons is more challenging than that of temperature.
In other words, the radius for which half of the energy is deposited at a radius smaller than it.In table 4 the discharge radius is given for the first discharge in air and H 2 /air at two pressures and three energy depositions.From the table, it can be seen that the discharge size grows slightly in size for increasing discharge power.This can be explained by electron-ion recombination reactions: As the deposition energy increases the maximum electron concentration in the center of the discharge is expected to increase.The increased electron concentration also implies an increased concentration of positively charged ions.As a result, the rate of progress of third-body recombination reactions like is significantly increased in the center of the domain.This slows down the power deposition if the electron concentration becomes large, which slightly increases the discharge size.This effect is not included in our analytical analysis of equations ( 28), ( 31) and ( 33), but as can be seen it only has a minor effect.When comparing the discharge size of the ϕ = 0.7 to the ϕ = 0.0 cases, it can be seen that they are very similar.This is not surprising as our ion chemistry does not include the influence of hydrogen ionization.As such, the mixture acts like a slightly more difficult to ionize air mixture.This implies that the reduced electric field is higher in the entire domain (see figure 15), but as the temperature is uniform in the first discharge the shape of the E/N profile remains similar.As the E/N profile determines the discharge shape, the shape of a discharge in ϕ = 0.7 is expected to be similar to that in ϕ = 0.0.Finally, the size of our discharges is of a similar order of magnitude as those measured by Patel et al [37]; a few hundred of micrometer.
In figure 17 the ratio of discharge size of the second to first discharge is shown as a function of the temperature ratio.The temperature ratio is defined as the ratio of temperatures in the center of the domain just before the discharge (in all cases T c,1 = 300 K).For reference, the spatial distribution of the temperature at the inception of the second discharge can be seen in figures 12 and 13.A low ratio of discharge sizes, indicates that the second discharge is more contracted with respect to the first discharge.In the graph data for ϕ = 0.0 and ϕ = 0.7 is shown at both 500 and 750 mbar, each data point represents the result of a single simulation.The ratio in temperatures is the result of the discharge energy and possible chemical heat release.The discharge energies are identical to those in table 4. For all presented cases the same trend is observed at the lowest temperature ratio, i.e. discharge power, the contraction of the plasma is weakest, while at the highest the contraction is strongest.Increasing the pressure reduces the temperature ratio, which can easily be explained as more gas needs to be heated.However, the contraction of the discharge is not significantly altered by the modification of pressure.
In all cases, the discharge contraction does not decrease linearly by the discharge energy (or temperature ratio).This can be explained via the previously discussed electron-ion recombination reactions.As the discharge energy becomes higher-and as a result, the second discharge radius becomes smaller-the electron concentration increases significantly.
The electron-ion recombination rate increases mostly in the center of the domain slowing down further ionization and thereby widening the discharge.
For the cases in which hydrogen is present, the chemical heat release increases the temperature ratio.At the same time, the hydrogen reduces the discharge contraction for the lowest energy setting.This can be explained via ion chemistry, which is based on that of air, meaning, the ϕ = 0.7 cases are chemically similar to a lower-pressure air case.However, the presence of hydrogen does imply that the ratio of temperatures becomes higher for the same discharge energy.For the higher energy cases, this results in a stronger contraction of the plasma.This implies that the presence of chemical heat release can enhance the discharge contraction.

Conclusion
In nanosecond discharges the energy deposition is sufficiently fast that it occurs within acoustic and diffusive transport time scales.This implies that the distribution of the energy in the discharge is exclusively the result of electron chemistry and initial conditions.The classical explanation of the heatinginduced plasma instability can, therefore, not be applied.However, in repetitively pulsed discharges a similar instability is present.In this article, we have investigated the coupling between discharges in a nanosecond type discharge by detailed numerical modeling.This investigation is performed for both an air mixture and an air-fuel mixture.We conclude that two effects are at play.First, a variant of the classical heating-induced effect, i.e. gas heating from previous discharges increases the reduced electric field which (locally) enhances the electron production rate.As a result, the volume of the discharge tends to reduce when multiple pulses are applied.Second, residual electrons of the previous discharge might not be distributed homogeneously.This is the result of varying electron attachment and detachment rates with both temperature and composition.From our analysis, we conclude that the temperature-induced effect is significantly stronger than that of residual electrons.In the future, it would be interesting to perform this study with a 2D or 3D numerical domain to confirm our observations.
Accurately predicting the effect of pulse-to-pulse coupling is crucial to further our understanding of plasma-assisted ignition.As we have demonstrated, the plasma size is extremely sensitive to the temperature field produced by the previous discharge.The size of the plasma determines the energy density, which impacts the temperature increase and radical production due to a discharge.As both the temperature and radical concentration have a significant impact on the rate of oxidation (or chemical heat release), a positive feedback loop is present.This implies that hot spots produced by one discharge are expected to receive a higher energy loading in consecutive discharges, accelerating the oxidation chemistry in the spots at which it is already fast.

Figure 1 .
Figure 1.Schematic overview of the experimental configuration of Patel et al[35][36][37], the pink area is used to indicate where the discharges take place.

Figure 2 .
Figure 2. Illustration of the modeled electric field.The bottom electrode has a potential of φz = 0 and the top electrode is modeled by φ elec .In the center the electric potential crosses the discharge gap, which has a distance dgap.Outside the center line, the electric potential also encounters the dielectric material surrounding the electrode lengthening its path with d dielec .

Figure 3 .
Figure3.The temporal evolution of the reduced electric field (left), at r = 0.0, 2.5 and 5.0 mm.The temporal evolution of the total deposited power inside the domain (right).Both for a discharge at 500 mbar, 300 K and an energy deposition of 875.0 µJ.

Figure 4 .
Figure 4. Mole fraction of electrons at r = 0.0 and 2.5 mm, for a discharge at 500 mbar, 300 K and an energy deposition of 875.0 µJ.

Figure 5 .
Figure5.The spatial distribution of the time-integrated power deposition for the first discharge in air for three cases: 1) 875 µJ at 500 mbar, 2) 437.5 µJ at 500 mbar and 3) 875 µJ at 750 mbar.

Figure 6 .
Figure 6.Fraction of energy in different chemical modes (vibrational excitation, N 2 ( 3 ) excitation and O 2 dissociation).In the left frame the spatial distribution is shown at 13.8 ns (peak energy of deposition), in the right frame the temporal evolution of the centerline is shown.The dotted line in the right frame is for reference and indicates 13.8 ns.The solid lines are at 875 µJ and the dashed lines at 437.5 µJ, all lines are at 500 mbar.

Figure 8 .
Figure 8.The acoustic pressure P ′ = P − P 0 at three different time instances for the 875 µJ at 500 mbar (left) and the 437.5 µJ at 500 mbar (right) case in air.

Figure 9 .
Figure 9.The fraction of energy in different chemical modes for an H 2 /air mixture at ϕ = 0.7.The discharge has an energy of 875 µJ and the gas has an initial pressure of 500 mbar.The solid lines are at 875 µJ and the dashed lines at 437.5 µJ, all lines are at 500 mbar.

Figure 11 .
Figure 11.Ion mole fraction at the centerline of the domain after a discharge of 875 µJ at 500 mbar in air.Left frame shows the negatively charged species and right frame the positively charged.

Figure 12 .
Figure 12.Distribution of the ions in air just before the second discharge t = 1/3 KHz ≈ 0.333 ms (left frame) and temperature (right frame).The ambient pressure is 500 mbar and the discharge has an energy deposition of 875 µJ.

Figure 13 .
Figure13.Distribution of the ions in H 2 /air at ϕ = 0.7 just before the second discharge t = 0.333 ms (left frame) and temperature (right frame).The ambient pressure is 500 mbar and the discharge has an energy deposition of 875 µJ.

Figure 14 .
Figure 14.First and second pulse spatial distribution of energy deposition in air and air/H 2 at ϕ = 0.7, both are at 500 mbar and 875 µJ.

Figure 15 .
Figure 15.The reduced electric field at peak energy deposition for the first and second discharge in air (left) and air/H 2 (right).

Figure 16 .N φ 2 elec
Figure 16.The effective electron production rate in O 2 , N 2 and air.

Figure 17 .
Figure 17.The ratio of the second to first discharge size as function of the centerline temperature ratio.Results are for both ϕ = 0.0 and ϕ = 0.7 at two pressure and three different discharge energies.Note, each point is an individual simulation, lines are only shown to indicate trends.

Table 3 .
Main simulation parameters used.