Simulation study on nitrogen vibrational kinetics in a single nanosecond pulse high voltage air discharge

High voltage nanosecond pulsed discharges (NPDs) have received great attention due to envisaged application in the fields of plasma assisted combustion, high-speed flow control and biomedical treatment.^1–4 The understanding and utilization of NPDs require fundamental knowledge of kinetic phenomena such as breakdown development,⁵ molecular energy transfer and chemical reactions,⁶ and hydrodynamic expansion^7,8 in non-equilibrium plasmas. One open question is the deposition of discharge energy in different internal degree modes of air plasma.^9,10 Basically, at the initial rising phase of the pulse, a high peak electrical field leads to gas ionization, efficient generation of electronic excited levels, and dissociation of molecule. After gas breakdown, the moderate fields favor the excitation of nitrogen vibrational levels.¹¹ The time evolution of translational temperature is associated with two-stage mechanism:^12,13 a fast heating regime related to quenching of electronic levels, and a slow heating regime related to nitrogen vibrational-translational (V-T) relaxation. As nitrogen vibrational levels serve as an energy “reservoir” in air plasma, the need of proper characterization of vibrational kinetics is unquestioned.

The progress in this area is indebted to the development of useful tools including Coherent Anti-Stokes Raman Scattering (CARS) and Spontaneous Raman Scattering (SRS).^14,15 Messina et al.¹⁶ found an increase of vibrational temperature at about 30 μs in the afterglow of an atmospheric air NPD using CARS, and attributed it to collisions with ions. However, a follow on simulation study¹⁷ interpreted it as redistribution of vibrational quanta by vibrational-vibrational (V-V) exchange. In a similar study, Pendleton et al.¹⁸ observed an increase of N₂ vibrational population (v = 1, 2) at roughly 5 μs following discharge. They concluded that energy pooling reactions such as N₂(A ³Σ) + N₂(A ³Σ) → N₂(X ¹Σ, v) + N₂(C ³Π) was the likely reason, which was proposed previously in pure nitrogen discharge by Devyatov et al.¹⁹ Later on, Montello et al.²⁰ used picosecond CARS to study N₂ vibrational levels and rotational temperature in NPD in nitrogen and air at a pressure of 100 Torr. The total vibrational quanta (v = 0-9) was found to increase significantly until about 10 μs after discharge. A phenomenological model of electronic-to-vibrational (E-V) energy transfer was provided to explain the results. More recently, Lo et al.²¹ reported similar phenomena after a NPD in atmospheric air, and also its spatial dependence by SRS. The above behavior was explained by a convective transport induced behind the shock wave following fast gas heating.²² 

A growing number of kinetic modeling is developed to address the experimental data for example, the puzzling results of increasing vibrational quanta after the discharge pulse.^6,20,23 However, a complete description of vibrational kinetics is rather challenging, in particular of time-resolved measurements. The present model, apart from the usually considered nitrogen V-V, V-T, and vibrational de/excitation by collision with plasma electrons (e-V), also takes into account the E-V and chemical-vibrational (Chem-V) processes by a statistical approach of state-specific reaction rates. The purpose of the work is not to compare the model and previous ones to show whether an advantage in more accuracy exists. Instead, it is designed to study the role of E-V and other processes in producing additional vibrational quanta, which is yet not well established by previous models. The vibrational kinetics is then coupled self-consistently to the gas thermal balance equation, electron energy distribution function (EEDF) solver, and plasma chemical kinetics. The modeling results were compared with the in burn (∼100 ns) and afterglow phase (∼10 ms) of a pin-to-pin NPD in 100 Torr flowing dry-air.²⁰ Good agreement with experimental results was obtained for time-dependent vibrational and translational temperatures. The modeling result showed a non-negligible rising of N_q in the early afterglow when E-V and Chem-V mechanism were included. The quantitative discrepancy is discussed in light of zero-dimensional modeling and real experimental conditions.

A quasi-zero-dimensional multi-temperature (translational temperature of neutral gas, ions and electrons, vibrational temperature of N₂ and O₂) model ZDPlasKin²⁴ is used to describe the air plasma kinetics in NPD. Such model is less computationally expensive than higher dimensional models, and retains the rather uniform character of the diffuse plasma filament.

The number density n_A for gas species A is described by the kinetic equation

where the source term on the right side is a sum of corresponding plasma physical/chemical reaction sources S_Aj for all processes involving species A from j = 1 to j = j_max. For example, this source term for reaction

can be expressed using reaction rate R_j and rate constant k_j

The inverse reaction process is defined manually if it exists. The list of species and reactions are converted automatically to a system of ordinary differential equations, and solved numerically using ZDPlasKin module. The basic model includes 55 air plasma components and 784 reactions. Table I shows the components of air plasma, and Table II shows the excited N₂ and O₂ effective states clustered around actual molecular levels (due to availability of reaction rates). Table III shows a classification of all reaction processes for nitrogen-oxygen mixtures.^25–34 

The vibrational kinetics is calculated using a state-to-state master equation model, with all reactions for nitrogen listed in Table IV. The basic model of N₂ vibrational kinetics includes e-V, V-V and V-T processes involving first 8 vibrational levels N₂(X¹, v0-8), while O₂ vibrational kinetics includes e-V, and V-T processes involving first 4 vibrational levels O₂(X³, v0-4). The master equation for N₂ vibrational kinetics in the basic model can be written as follows:

The rate constants k_eV for N₂ e-V processes are obtained by solving electron Boltzmann equation with corresponding cross sections σ_e-V. σ_e-V for excitation of ground level are adopted from Lxcat database,²⁶ and excitation of excited levels are constructed using a scaling relation σ_e-V(w, v) = σ_e-V(0, v − w)/(1 + 0.05v) proposed by Gordiets et al.²⁷ Because the vibrational-vibrational exchanges between N₂ and O₂ (V-V’) are much slower than V-V processes in N₂, V-V’ processes are not considered. The multi-quanta transitions for V-V and V-T have limited effect at gas temperature lower than 1000 K, therefore are also not considered. The forced harmonic oscillator free rotator (FHO-FR) model is used for N₂ V-V rates,³⁵ with $k 1 , 0 0 , 1$ taken from a direct experimental measurement.²⁸ The simple HO model is used for N₂ V-T rates, with k_1,0 taken from Ref. 25 for relaxation by N₂, N, and from Ref. 29 for O by fitting of experimental data over a wide range of temperatures.

In an extended model of vibrational kinetics, higher vibrational levels v = 9-12 are also added, to check whether downward V-V processes influence the population of lower vibrational levels v = 0-8. In order to describe the increase of N_q in the afterglow adequately, E-V and Chem-V processes which may produce additional vibrational quanta are added. However, the rate constant for E-V and Chem-V is rather scarce for a kinetic modeling. A statistical approach from information theory is used to calculate the state-specific reaction rate. For E-V processes such as $O ( 1 D ) + N 2 →O ( 3 P ) + N 2 ( v ̄ = 2 . 2 )$⁠, rate constant k_v is supposed to follow a Poisson distribution³⁶ 

where $v ̄ =2.2$ is the measured average vibrational level in the products (i.e. about 30% of released energy goes into vibrational excitation and the remaining goes to gas heating).³⁷ Above theory was found to agree with the experiment result of E-V process for CO vibrational excitation.³⁸ If $v ̄$ is unknown from any measurement, such as in reaction N₂(A ³Σ, B ³Π, C ³Π, a¹ ′Σ) + O₂ → N₂(X, v) + O + O, the model assumes similar fraction of released energy into vibrational excitation, and $v ̄$ is calculated to satisfy

It should be noted that above E-V processes are the dominant mechanisms for quenching of electronic states in air discharge, while the energy pooling processes dominates in pure nitrogen.^20,30 The rate constant for Chem-V process such as $N ( 4 S ) +NO→O+N2 ( v ̄ = 3 . 4 )$ is also calculated by a similar method,³⁹ with fraction of the released energy transferred into vibrations known from measurements.⁴⁰ The surface-vibrational (Surf-V) process is included in the 0D model by the characteristic diffusion time $τ diffusion − 1 ≃D/ ( d / 2 ) 2$ to describe particle diffusion and reaction at the wall (D diffusion coefficient, and d diameter of discharge filament). After incorporating E-V, Chem-V and Surf-V processes in the extended model, the complete air plasma kinetics include 59 species and 1045 reactions.

While the reduced electric field E/N and electron density n_e are two governing factors determining the energy partition in non-equilibrium plasma, they were input from the measured gap voltage (deduce the voltage drop on cathode layer) and current density by parametric equations.⁴¹ The EEDF is calculated by a steady-state two-term expansion solver Bolsig+⁴² to obtain rate coefficients of electron impact excitation, and dissociation rates. This is justified by much shorter EEDF relaxation time compared with the characteristic change time for applied voltage.⁴³ Because the shape of EEDF strongly influences the populations of excited states and radicals in the discharge,¹¹ Coulomb collisions between electrons are invoked when electron mole fraction is above 10⁻⁵. The second kind collisions (super-elastic collisions) from vibrational excited states and first level of electronically excited states are included for an accurate description of the tail of EEDF.^44,45 The density of 25 species, namely N₂(v0-12), N₂(A³), O₂(v0-4), O₂(a¹), N, O, NO, O₃, N₂O, and 169 electron impact processes (shown in Table V) are used for Bolsig+.

The variation of gas temperature T_g considers the thermal energy balance between power release in plasma physical/chemical reactions and heat conduction:

In Eq. (7), N_g is the gas number density, C_p is the heat capacity under constant pressure, λ_g is the thermal conductivity, and W_R is the released power transferred to gas heating. In 0D simulation of spatially averaged gas temperature, the heat conduction term can be written as³⁰ 

where the wall temperature T_wall ≈ 300 K is assumed to be roughly the environment temperature, and R ≈ 3.5 mm is about the same as the radius of discharge electrode.²⁰ The temperature-dependent heat capacity and conductivity are adopted from Ref. 30. The gas heating mechanisms^29,46 are listed in Table VI (the detail information about dominant E-V and Chem-V processes are marked in bold type).

The ion temperature differs from the gas temperature due to ion drift in the electrical field, and may influence the reaction rate involving ions:⁴⁷ 

Taking ion N⁺ as an example, the reduced mobility μ_iN = 8.0 × 10¹⁹ (V ⋅ cm ⋅ s)⁻¹ is retrieved from measurements.⁴⁸ 

The simulation in this work considers the experiment of a NPD in 100 Torr dry-air (78.5%N₂-21.5%O₂) carried out in Ref. 20 where time resolved N₂ vibrational and rotational temperature were measured in a 100 ns pulse and subsequent 10 ms afterglow. The discharge was in a point-to-point geometry with a gap distance of 1 cm, produced with a peak voltage of 9 kV and peak current of 48 A. The typical calculation time for complete air plasma kinetics using a modern workstation is around 30 hours, with 0.5 ns time step in discharge and 50 ns time step in the afterglow. During the 100 ns in burn phase, Bolsig+ calculates non-Maxwellian EEDF automatically and incorporates values of electron transport and rate coefficients into ZDPlasKin self-consistently. Between 100 ns to 10 ms when the discharge pulse is turned off, a Maxwellian EEDF is assumed to be mainly determined by super-elastic collisions with N₂ vibrational excited states.

Fig. 1 presents the simulated EEDF f_e during different time of discharge. At time t = 1 ns when applied voltage pulse rises to peak, the reduced electric field inside of the plasma is estimated to be 240 Td (1Td = 10⁻¹⁷ V/cm²). As a result, the tail of EEDF is rather high, which favors ionization, dissociation and excitation of electronic states. The electron density rises to about 10¹⁴ cm⁻³ by avalanche ionization, and breakdown occurs. The applied voltage quickly drops due to an increase of gas conductivity. Then, the main part of discharge is relatively stable at electric fields of about 120 Td. The shape of EEDF shows a remarked drop at electron energy ε_e around 2 eV during this wide time range, because large amount of deposited energy goes into N₂ vibrational mode.⁴⁹ At the end of discharge pulse, the reduced electric field is around 60 Td with an average electron temperature T_e around 1 eV (⁠$T e = 2 3 ∫ ε e f e dε$⁠).

Fig. 2(a) presents the simulated vibrational temperature of the first level for N₂ and O₂, and gas temperature. The dashed line in a light grey background depicts the in burn regime, and the solid line depicts the afterglow. Because the vibrational distribution function (VDF) of N₂ and O₂ shows a “bimodal” structure different from Boltzmann distribution,⁵⁰ two different vibrational temperatures are defined

It can be observed that T_v01 for N₂ is quite higher than T_v01 for O₂ in the whole time range, while the latter is close to gas temperature. This result is attributed to the low electron-impact vibrational excitation cross sections of O₂ compared to N₂, and the strong V-T relaxation of O₂ by O atoms. Fig. 2(b) and 2(c) compares the simulated and experimental results for N₂ vibrational temperature and gas temperature, respectively. Both results show quantitative agreements during in burn regime and in the afterglow. The above agreements prove the validity of proposed vibrational kinetics and gas heating mechanisms. During the in burn regime t ≤ 10⁻⁷ s, the rapid rise of T_v01 for N₂ is mainly caused by e-V excitation process; between 10⁻⁷ to 10⁻⁴ s in the afterglow, the rise of T_v01 is mainly due to downward V-V process; between 10⁻⁴ to 10⁻² s in the afterglow, the decrease of T_v01 is mainly due to V-T relaxation by O atoms.

For the variation of gas temperature, a two-stage energy thermalization mechanism is confirmed in the afterglow.^12,13 A rapid gas heating rate about 10 k/μs between 10⁻⁷ to 10⁻⁶ s mainly corresponds to quenching of electronic states N₂(A³, B³, a^‘1, C³) by O₂ and quenching of O(¹D) by N₂. A slow gas heating rate about 0.5 k/μs between 10⁻⁴ to 10⁻³ s mainly corresponds to N₂ V-T relaxation by O atoms, which accompanies the decrease of T_v01 for N₂ during this time range. The decrease of gas temperature at time later than 10⁻³ s is caused by heat conduction. In our quasi-0D model, the isobaric approximation was used to describe the temporal evolution of gas temperature both in plasma and in afterglow. The above approximation was generally true for discharge in flow conditions or even in static conditions when energy deposition is not too high.³⁰ However, the pressure may change in NPD due to gas dynamic expansion induced by fast gas heating. According to a discharge model coupled with compressible Navier–Stokes equations for the same conditions,^6,23 the pressure at the discharge center line may increase to about 150 Torr at 1 μs, followed by a decrease to 90 Torr at 5 μs, and remain 100 Torr after 10 μs. At the peripheral area, the pressure change is around 10%. Therefore, the isobaric approximation is valid in most areas for a long time scale, which can be justified in Fig. 4(c). Although at time around 1μs the isobaric approximation was not valid to a certain extent, the variation of gas temperature would not induce considerable uncertainty because the heating terms in this time period is limited (just between the two thermlization stage). Therefore, the heavy particle kinetics involving temperature-dependent rate constant will not be influenced significantly. The effects of isobaric approximation on vibrational kinetics will be discussed later separately.

Further information on time resolved VDF for N₂ can be found in Fig. 3 and 4. The relative population of N₂ vibrational state n_v/[N₂] for levels v = 0-1 are predicted quite well with the model, as can be seen in Fig. 3 over the whole time scales. The increase of N₂ vibrational populations v = 0-1 in the initial phase of afterglow is mainly caused by downward V-V transition. However, the model under-predicts n_v/[N₂] for levels v = 2-4 for time later than 2 × 10⁻⁶ s. The detail of this discrepancy can be found in Fig. 4. In Fig. 4(a), the model over-estimates the vibrational populations v = 5-8 at the end of discharge pulse. This was possibly due to lack of reliable data for corresponding electron-impact cross sections (an inaccuracy of 30% was estimated for excitation of higher vibrational levels). In Fig. 4(b), the model under-predicts the vibrational populations v = 2-8 in time later than 10 μs, but not in the situation at 1μs (a slight over-prediction is observed by the model). This suggests that additional vibrational loading mechanism, other than E-V and Chem-V energy transfer, exists between 1 μs to 10 μs, but is not covered in the 0D model.

Fig. 5 compares the simulated and experimental vibrational quanta in burn and in afterglow. The vibrational quanta can be accurately predicted for e-V excitation in burn and V-T relaxation in afterglow by the model. However, the rise of vibrational quanta observed in the experiment is not reproduced by the model. Note that, incorporation of E-V mechanism gives a non-negligible rise (about 2%) of vibrational quanta during time scale 10⁻⁷ ∼ 10⁻⁶ s. It can be seen that the model over-estimates the vibrational quanta at the end of discharge pulse, i.e. it under-estimates the restored energy in electronic excited state. As a result, the E-V energy transfer is under-estimated in the current model (the use of effective electronic states may be another reason for this under-estimation). However, the E-V mechanism alone could not explain the increased amount of vibrational quanta and its time scale. The time evolution of number density for main electronic states and chemical species in air plasma are demonstrated in Fig. 6. The stored energy in electronic states is calculated from number density in Fig. 6 and energy levels in Table II, and estimated to be no more than 60 meV per molecule. For E-V reactions like N₂(A ³Σ) + O₂ → N₂(X, v) + O + O, a large fraction of N₂ excitation energy (6.17 eV) is consumed by O₂ dissociation (5.12 eV), then the remaining released energy is partitioned among gas heating (assumed about 70%) and N₂ vibrational excitation. The available energy from electronic levels into vibrational excitation is so small, and in contradiction to the measured 60% rise of N_q. Besides, the quenching processes of N₂(A³, B³, a^‘1, C³) by O₂ and O(¹D) by N₂ is so fast that the number density for all electronic excited states drops to at least 1/10 of their peak values at 1 μs. This disagrees with the fact that N_q increased significantly between 1 μs to 10 μs. Therefore, the rise of N_q at time later than 1 μs is unlikely to occur by E-V mechanism.

The Chem-V mechanism may also produce vibrational excited N₂, but the mole fraction of reactants N(⁴S, ²D) and NO is so low (less than 10⁻³) to give an observable effect. The electronic excited nitrogen atom N(²D) is quickly quenched by oxygen molecule O₂, and has little chance to react with NO to produce vibrational excited N₂. The ground state nitrogen atom N(⁴S) may have sufficient time τ (about 20 μs) to react with NO, and its contribution to vibrational loading can be estimated as $k [ N ] [ NO ] v ̄ τ/ [ N 2 ] ≃3×1 0 − 3$ (species in bracket denotes its number density, k is rate constant). In fact, during the time when Chem -V processes are believed to occur, the temporal variation of Nq remains almost constant. The contribution to rise of N_q from higher vibrational levels v > 8 is also ruled out, because our extended model has already incorporated levels v = 9-12 with e-V cross sections from Biagi database.⁵¹ Above result confirms the previous conclusion,²⁰ where N₂ e-V cross sections for v > 8 were estimated using the results of Huo et al.⁵² 

We note that the authors²⁰ observed a compression wave due to rapid gas heating at t = 1 ∼ 10 μs and an expansion of discharge filament diameter from 2 mm to 4 mm. In their recent model,²³ one radial coordinate was incorporated in the afterglow to study this effect and the random transverse motion of the filament jitter in experiment. The result of 1D model showed that in the time range between 1 ∼ 5 μs, the vibrational quanta of N₂ increased by 20%. The above explanation reminds us to look at the spatiotemporal data in a similar experiment.^21,22 The variation of N_q in a NPD in atmospheric air by Lo et al.^21,22 is calculated based on measured data and plotted in Fig. 7. While the N_q on-axis decreases monotonically with time in the afterglow, the N_q off-axis at a radial position of r = 0.72 mm increases at time t = 1 ∼ 10 μs. The reason was explained by convective transport associated with the compression wave. Because the vibrational temperature on-axis is higher than that off-axis (also presented in Fig. 15 for a NPD in nitrogen²⁰), the gas-dynamic expansion of the channel would sweep the vibrational “hot” N₂(v) from the axial central line to the periphery,⁵³ leading to an increase of vibrational temperature and quanta observed. This hydrodynamic expansion not covered in current 0D model, together with the under-estimated E-V mechanism, may explain the quantitative discrepancy between modeling and experimental results in Fig. 5.

We presented a modeling study on N₂ vibrational kinetics during a NPD in dry-air (78.5% N₂/21.5% O₂) and the subsequent afterglow. The model is based on self-consistent solution to time-dependent gas thermal balance equation, coupled to electron, vibrational and chemical kinetics. The E-V and Chem-V mechanisms are included in the vibrational kinetics by state-specific reaction rates using a statistical approach. Good agreement between modeling and experiment has been obtained. The quantitative discrepancy in rise of vibrational quanta in the afterglow is attributed to the under-estimated E-V energy transfer, and the hydrodynamic expansion not covered in this 0D model. This work will be improved in the future by considering the following issues: a higher-dimensional model in predicting the electrical field using Poisson equation; incorporating hydrodynamic expansion phenomena associated with rapid gas heating; a detailed analysis of the energy partition in E-V energy transfer process and a complete database for accurate cross sections/rate constants.

The authors thank Igor V Adamovich from Ohio State University for helpful discussions on experimental and modeling results. This work is supported by the National Natural Science Foundation of China (Grant No. 11505015).