Spatio-temporal analysis of power deposition and vibrational excitation in pulsed N2 microwave discharges from 1D fluid modelling and experiments

Vibrational excitation of N2 beyond thermodynamic equilibrium enhances the reactivity of this molecule and the production of radicals. Experimentally measured temporal and spatial profiles of gas and vibrational temperature show that strong vibrational non-equilibrium is found in a pulsed microwave discharges at moderate pressure (25 mbar) in pure N2 outside the plasma core and as an effect of power pulsing. A one dimensional radial time-resolved self-consistent fluid model has been developed to study the mechanism of formation of vibrationally excited N2. In addition to the temperature maps, time-resolved measurements of spontaneous optical emission, electron density and electron temperature are used to validate the model and the choice of input power density. The model reveals two regions in the plasma: a core where chemistry is dominated by power deposition and where vibrational excitation starts within the first ∼10 µs and an outer region reliant on radial transport, where vibrational excitation is activated slowly during the whole length of the pulse (200 µs). The two regions are separated by a sharp gradient in the estimated deposited power density, which is revealed to be wider than the emission intensity profile used to estimate the plasma size. The low concentration of excited species outside the core prevents the gas from heating and the reduced quenching rates prevent the destruction of vibrationally excited N2, thereby maintaining the observed high non-equilibrium.

Vibrational excitation of N 2 beyond thermodynamic equilibrium enhances the reactivity of this molecule and the production of radicals.Experimentally measured temporal and spatial profiles of gas and vibrational temperature show that strong vibrational non-equilibrium is found in a pulsed microwave discharges at moderate pressure (25 mbar) in pure N 2 outside the plasma core and as an effect of power pulsing.A one dimensional radial time-resolved self-consistent fluid model has been developed to study the mechanism of formation of vibrationally excited N 2 .In addition to the temperature maps, time-resolved measurements of spontaneous optical emission, electron density and electron temperature are used to validate the model and the choice of input power density.The model reveals two regions in the plasma: a core where chemistry is dominated by power deposition and where vibrational excitation starts within the first ∼10 µs and an outer region reliant on radial transport, where vibrational excitation is activated slowly during the whole length of the pulse (200 µs).The two regions are separated by a sharp gradient in the estimated deposited power density, which is revealed to be wider than the emission intensity profile used to estimate the plasma size.The low concentration of excited species outside the core prevents the gas from heating and the reduced quenching rates prevent the destruction of vibrationally excited N 2 , thereby maintaining the observed high non-equilibrium.

Introduction
In recent years vibrational excitation of nitrogen molecules has been shown to increase the production of reactive radicals and to enhance the reactivity of the molecule [1].In particular, vibrational activation has been studied in the context of NO x production and has been shown to improve NO x yield in various types of discharges [2][3][4][5].The presence of oxygen in the mixture has been indicated as the main obstacle to obtaining the needed high degrees of vibrational non-equilibrium in N 2 -O 2 [6], suggesting that activation of N 2 prior to the reactive zone could be beneficial [3].
Recent experimental studies have shown the possibility of obtaining non-equilibrium conditions in nitrogen either by reducing power input and pressure or by exploiting the different time scales of vibrational excitation and quenching [7][8][9].The latter relies on variation of the exposure of the particles in the plasma to the power source by changing the input gas flow rate and/or providing the power in pulses of varying duration and offers more possibility of control, since pulsing frequency and duration can be externally imposed.
Non-equilibrium as a function of distance from the plasma core in continuous wave (CW) mode has been investigated experimentally, showing that values of the ratio between vibrational and gas temperature (T v /T g ) up to 4 can be maintained in the region outside the plasma core [7].This behaviour had been attributed to the low gas heating outside the core due to the low pressure and low power input.Using the same operational mode with N 2 -O 2 (50%-50%) feed gas, Kelly et al [10] have found through simulations that vibrational nonequilibrium of nitrogen molecules can be obtained by tailoring the mass flow rate and therefore the residence time of the gas in the hot plasma region.This was identified as one of the reasons for the enhanced NO x yield.
Another way to exploit the different timescales of chemical reactions, vibrational excitation, and gas heating is through the modulation of the input power in time by using a pulsed power input [6,8,11,12].The end of the non-equilibrium time window measured in a 25 mbar pure nitrogen discharge by van de Steeg et al [11] was attributed to the raise in N atoms number density causing fast quenching of vibrations.In the works by Bahnamiri et al [6,12] the matching between residence time and pulse duration was proposed as a way to achieve optimal NO x yield through vibrational non-equilibrium; this was maintained during the whole pulse due to the vibrationaltranslational (V-T) collision rates being ∼100 times lower than other processes favouring vibrational energy loading into N 2 .
In the work by Kelly et al [10], 0D modelling was employed to identify the main reaction paths [4,10].The complexity of energy transfers between the different degrees of freedom and the addition of transport phenomena calls for validation of the models used to predict reactor efficiency.In particular, power deposition profiles and transport characteristic times are subject to approximations in 0D models and chemical kinetics data has not been extensively validated for the intermediate pressure regimes that can be maintained in microwave (MW) discharges.
The value of deposited power density (P d ) is generally used as input in a typical model of MW discharges; the value of the reduced electric field (E/N) needed to solve the electron Boltzmann equation (EBE) is then estimated from it assuming Joule heating of the electrons [8,10].This input parameter is calculated using the measured absorbed power and the plasma volume estimated from optical emission intensity profiles, assuming that there exists a linear dependence between absorbed power and emission intensity.This assumption has been previously tackled for CO 2 MW discharges by Viegas et al [13].In that work it was shown that the emission profile is generally narrower than the applied power density profile.Groen et al [14] estimated the reduced electric field radial profile in a CO 2 MW discharge using measured electron density and gas temperature, and the principle of impedance matching.They showed a non-concave electric field and reduced electric field profile, which, depending on the profile of electron temperature and density, might result in a non-concave profile of deposited power density.The large uncertainties linked to the estimation of P d lead to large variations in model results in 0D and one dimensional (1D) models [8,[15][16][17].
Modelling of MW discharges is also affected by the treatment of both radial and axial mass and energy transport.Transport in the radial direction in 0D models is often treated by means of characteristic diffusion times that are determined by assuming that charged particle number density profiles are described by a zeroth order Bessel function [18], which is true for systems dominated by diffusion.Similarly, a parabolic profile is assumed for gas temperature to determine losses of heat due to conduction to the reactor walls [19].These assumptions were recently checked by Viegas et al [20] for O 2 DC glow discharges up to 13 mbar, revealing their validity in that pressure range.Modelling of CW CO 2 discharges at higher pressures (above 100 mbar) with a 1D radial fluid model has shown non-parabolic T g profiles [16] and non-Bessel electron density profiles [13].1D fluid modelling of the same kind of discharge in pure nitrogen by Volynets et al [21] has shown parabolic profiles for T g up to 53 mbar.However, the non-uniform electric field in a MW discharge might affect the balance between transport phenomena and chemical kinetics, leading to radial profiles that are not solely determined by diffusion.Moreover, the presence of a vortex flow employed in MW discharges to stabilise and centre the plasma might affect particle balance and gas cooling at the edge of the discharge due to axial convective transport, which is negligible in glow discharges.Transport in the axial direction is usually included by considering a plug flow [4, 8, 15, 16], which, given the complex flow patterns in an MW discharge [22], might be underestimating the gas residence time within the hot plasma zone.
Uncertainties are present also on the chemistry set to be used for intermediate pressure (10-1000 mbar) discharges.Extensive work has been dedicated to the modelling of lowpressure discharges in pure nitrogen [23], resulting in a selfconsistent chemistry set that was recently employed to simulate CW experiments for p > 50 mbar [17].Given the role of N atoms in the heating of the gas in N 2 discharges [10,11], rate coefficients for dissociation mechanisms are often carefully treated in literature.0D modelling of a pulsed MW discharge at 25 mbar by Kelly et al [10], based on the chemistry set employed in a previous work by Van Alphen et al [8], includes dissociation by means of the Fridman-Macheret approximation, with the only recombination mechanism leading to the formation of N 2 in the electronic and vibrational ground state.More precise derivation of state specific dissociation for high temperature air plasmas including also multiquanta processes is done by means of the forced harmonic oscillator (FHO) theory used by Macheret and Adamovich [24] and later by da Silva et al [25,26].Low pressure models in pure N 2 include additional dissociation mechanisms involving electronically excited states of N 2 and vibrational-vibrational (V-V) collisions and recombination to the excited states N 2 (B), but neglect multi-quanta collisions with gas molecules, which become relevant as pressure and temperature are increased.Lastly, the work by Volynets et al [21] suggests the importance of excitation of N 2 to predissociative states as one of the main sources of atomic nitrogen.Moreover, in that work, a new set of reactions for dissociation is proposed, which neglects collisions between electronically excited and vibrationally excited N 2 .
Finally, due to the role of electronically excited species in ionisation, dissociation and heating [4,17], their creation and destruction mechanisms are also under investigation in literature, with particular attention to state-specific impact excitation from vibrational levels [27][28][29].
In this work, a 1D radially resolved model is validated against experimental measurements on a pulsed MW N 2 discharge.Time resolved measurements of electron density in the plasma core are used to validate the ionisation-recombination mechanisms included in the model.These measurements combined with temporal measurements of electron temperature and optical emission intensity are used to confirm the correctness of the power density input profile.Raman scattering measurements of vibrational and gas temperature with spatial and temporal resolution are used to validate the chemistry set employed in the model.Results obtained by implementing the dissociation schemes used by Kelly et al [10] and by Volynets et al [21] are compared.The relative importance of axial and radial transport within the timescale of the pulse is estimated to highlight model shortcomings.Finally, a detailed investigation on the reason underlying the mismatch between emission intensity and power deposition is presented.

Experiments
The MW (frequency of 2.45 GHz) 200 µs long power pulses used to create the discharge are provided by a solid state MW power source (Ampleon) at a frequency of 30 Hz inside a quartz tube of 27 mm of inner diameter and 30 mm of outer diameter [11].The discharge is created in pure N 2 at a pressure of 25 mbar.Due to practical difficulties in igniting the plasma, the power in between pulses is kept low (∼17% of the pulse power), but not zero.Matching of the plasma impedance to assure minimum reflected power at the beginning of the pulse is achieved with a three stub tuner.The fraction of absorbed power decreases with time during the pulse due to changes in plasma properties, specifically electron density.Figure 1 shows the forward, reflected and absorbed power during the pulse.A delay generator (SRS-DG645) controls the delay between the power source trigger and the firing of the laser.This allows time-resolved measurements.The diagnostics is provided by a 30 Hz frequency doubled Nd:YAG laser (λ = 532 nm) of 400 mJ focused on the centre of the wave guide.The light scattered by the plasma is collected perpendicularly to the beam (see figure 2).The collection optics for Raman spectroscopy consists in an f = 75 mm lens focusing the scattered light into a fibre bundle; stray light and Rayleigh light are filtered out using a 532 nm long pass filter positioned in front of the fibre array.For Thomson scattering measurements, light is collimated by a f = 100 mm lens onto a Bragg grating, which reflects the scattered Rayleigh light [9]; the light transmitted by the grating is focused into the fibre bundle by an f = 100 mm lens.Further reduction of stray light is obtained by a closing aperture positioned in front of the fibres.The light collected by the fibres is projected into a custom built Littrow spectrometer, equipped with an intensified charge-coupled device (iCCD) camera.The setup can be moved perpendicularly to the laser beam to collect light from radial positions different from the core, specifically at 3, 6 and 9 mm from the tube axis.
Thomson scattering is employed to obtain electron density n e , which is proportional to the intensity of the scattered signal, and electron temperature T e , which is proportional to its width [11,30,31].An example of a typical Thomson spectrum is shown in figure 3. N 2 rotational Raman spectrum at room temperature and 50 mbar of pressure was used as calibration, necessary to obtain an absolute value for electron density.
The collection optics for Raman scattering is positioned so that light is collected from a 1 cm long region in the axial direction.In the pulsed regime, plasma properties are not uniform, meaning that within the 1 cm region temperatures vary considerably.While the signal for the rotational spectra is intense enough to allow a spatially resolved analysis in the axial direction (2 mm of resolution), the vibrational spectra are analysed without this distinction, therefore providing an average value  of T v for the region of collection.An example of a typical vibrational Raman spectrum is shown in figure 3.
Finally, the broadband plasma emission is imaged with an ICCD camera and Abel-inverted to obtain the emission intensity radial profile.The latter is used to estimate the power density input.

Model
The 1D radial fluid model employed in this work is described in the work by Vialetto et al [16] and used to simulate CO 2 MW discharges in CW mode.Similarly to what is described in that work, in the present case the model solves a system of coupled mass balance equations for 61 heavy species: 46 vibrational levels of the N 2 electronic ground state, electronically excited states of N 2 (N 2 (A 3 , nitrogen atoms in their ground electronic state (N( 4 S)) and their electronically excited states (N( 2 P) and N( 2 D)), ionic species (N + , N + 2 (X 2 Σ + g ), N + 2 (B 2 Σ + u ), N + 3 and N + 4 ).Electron number density is calculated by imposing local quasi-neutrality.The gas temperature equation is also solved self-consistently.
The original model has been further developed to allow: • time-dependent mass balance and gas temperature solutions; • a time-dependent equation for the electron mean energy density and therefore self-consistent temporal evolution of electron mean energy, allowing employment of the local mean energy approximation [32]; • solution of the two-term stationary homogeneous EBE with Bolsig+ [33] at different nodal points and with updates in time.
These three additions and the input parameters required by the model are discussed in the following along with a general description of the code.The explicit dependence on time (t) and radial position (r) of each quantity will be introduced once and then not reiterated.

Mass balance equations
In 1D cylindrical coordinates, the time-dependent mass balance equation solved in the model for the ith species is: (1) where i is the species index, y i is its mass fraction, J i is its mass flux, described through the drift-diffusion approximation, m i and ω i are the mass and the production rate, and Ω i is its loss rate due to axial transport.For a precise definition of all the variables the authors refer to the work by Vialetto et al [16].Differently from that work, y i (r, t) is solved both as a function of time and space.Equation (1) is solved for all species except for electrons and N 2 (X, v = 0), whose number densities are imposed from charge-neutrality and conservation of mass flux densities, respectively, as is described in the following.
The last term on the right hand side of equation (1) describes either the loss or the gain of ρ i due to axial convective transport.This model assumes that the gas injected in the system consists exclusively of N 2 in the ground electronic and vibrational state.Thus the only species for which where L is the length of the plasma, which will be discussed in detail later on, and v axial (r, t) is the speed of the inflow gas and is calculated from the conservation of mass as: where Φ is the volumetric input flow rate, ρ env is the gas mass density of the environment (pure nitrogen gas at 300 K) and S tube is the section of the tube.Notice that equation (3) can be derived from: ˆ2π by assuming that Φ ρ env is uniform over the tube section and that mass conservation holds for each elementary ring of area π . The definition of velocity v axial as in equation (3) assumes that the flow velocity never changes sign for different radial positions, while previous works on CO 2 discharges have indicated the possible presence of a backflow and a recirculation of previously heated gas back to the core [22].This phenomenon would lead to a redistribution of energy, either in the form of internal energy or translational energy to the outer layers and would reduce losses in the core.
Unlike N 2 (X, v = 0) all other species are exiting the reactive region.Ω i (r, t) changes accordingly: This same approach was introduced by Viegas et al [15] to take into account axial transport in 0D models.
The mixture diffusion coefficient used to calculate the diffusive component of the mass flux J i (r, t) is calculated as [34,35]: where the binary diffusion coefficients D ij for species i in species j are calculated as indicated by Guerra et al [23] (equation (37) in the paper and adapted from Hirschfelder et al [36]) and depend on Lennard-Jones binary interaction potential parameters [36].In the state-to-state formulation used in this work, electronically and vibrationally excited states are considered as chemically distinct species, with their own binary diffusion coefficient.This approach is reminiscent of what is done in the work by Pintassilgo et al [37] where the electronically excited states have a distinct diffusion coefficient in N 2 and also of the work about the modelling of hydrogen internal energy transport by Bruno et al [38], where diffusive velocities are defined for each excited state.
and mass conservation is guaranteed by imposing: Finally, isobaric conditions are maintained by enforcing: where i runs over all heavy particles and the translational temperature of all heavy particles is assumed to be equal to T g .Note that equation ( 9) is obtained from Dalton's law of partial pressures by neglecting the electron contribution to the total pressure.This is a valid assumption in the conditions investigated in this work.The electron density is obtained from the quasi neutrality condition.

Mean electron energy density equation
The time evolution of electron mean energy density (n ϵ = n e ϵ, with ϵ being the electron mean energy) in the radial direction is described by the electron energy balance equation [35]: where q e is the elementary charge, E amb is the radial component of the ambipolar electric field, J e is the radial component of the electron mass flow and S ϵ is a term containing both sources and sinks for the electron mean energy density due to external power input and collisions.The mean energy density flux (J nϵ ) is defined as: where D nϵ is the energy diffusion coefficient and µ nϵ is the energy mobility; both are calculated by Bolsig+.The term S ϵ is calculated as: in which the power density is indicated by P d , the second term takes into account the power lost due to inelastic or elastic collisions, which is calculated by Bolsig+ and the last term takes into account the mean energy losses due to axial transport.

Gas temperature equation
The temporal and spatial evolution of the gas temperature T g (r, t) is described by the energy balance equation: In the left hand side of equation (13) c p is the mixture averaged molar heat capacity at constant pressure, calculated with data from Burcat and Ruscic [39] and n m is the molar density of the mixture.In the first term on the right hand side λ(T g ) is the thermal conductivity of the mixture [40].Both c p and λ are calculated considering the contributions of N and N 2 , without distinguishing their excited states.Given the limited concentration of atomic nitrogen, these quantities are mainly determined by N 2 properties.
The second term on the right hand side of the gas temperature equation takes into account heat losses due to the mass flux of the ith species with enthalpy h i [39].The term Q in is the net contribution of chemical reactions to gas heating.This term is calculated as: where Q V-V and Q V-T are due to vibrational kinetics; their form is detailed in the works by Pintassilgo and Guerra [41,42].The last term in Q in takes into account elastic collisions with electrons and has already been introduced in equation ( 12) as a loss term for the electron energy density.Q chem contains heat exchanges due to chemical reactions, calculated as: where the sum runs over all processes p, the product runs over all reactants in reaction p, n i,p is the number density of said reactants, k p is the rate coefficient of the reaction p and ∆H p is the energy released or absorbed by the processes.
A detailed list of all the exothermic reactions considered can be found in a previous work by the authors of this paper [17].
Lastly, term Q axial takes into account gas cooling due to convective transport and Q rad is the heat loss due to radiation, both are describes in previous publications [16,17].
Note that, differently from the work by Vialetto et al, in equation ( 13) the input power density does not appear as a source of heat.Instead it appears in the electron energy density equation as an electron heating term.The conversion to translational energy happens thanks to exothermic reactions involving excited species.

Boundary conditions
Owing to the cylindrical symmetry of the system, homogeneus Neumann boundary conditions are considered at r = 0 for all the equations solved in the model.
At the tube wall, a flux for neutral and charged particles has been used in accordance with the work of Hagelaar et al [43]: where γ i is the fraction of particles of species i that is adsorbed at the wall.This value varies depending on the species and is indicated in a previous work by the same authors of this paper [17].The boundary condition for the mean electron energy density is formulated in a similar way [44,45]: where the following relation has been used: with ϵ w being the energy carried away per electron, which is equal to 2T e = 4 3 ϵ [44].The heat flux at the quartz tube wall is assumed to be [16,46]: where R in = 13.5 mm is the inner radius of the tube, R out = 15 mm is the outer radius of the tube, γ w is the thermal conductivity of quartz (1.4 W m −1 K −1 ) and T g,env is the temperature of the environment (300 K).
3.5.Input parameters 3.5.1.Power density from optical emission.The deposited power density used as input in equation ( 12) is usually estimated assuming it overlaps exactly with the optical emission intensity radial profile (I(r, t)).In the conditions investigated in this work this quantity appeared to reach a maximum at r > 0 mm and decrease at positions closer to the axis of the tube.To obtain optimal fit of the experimental values, the optical emission intensity was assumed to follow a decaying exponential profile for r < r 0,I , where r 0,I represents the position of the peak, and a step-function profile for r > r 0,I : The parameter δ p,exp measures the characteristic length of decay of the emission intensity within the core of the plasma, δ p describes the position of the edge of the step function with respect to r 0,I and σ p measures the characteristic length of decay of the emission intensity after the peak is reached.This shape was chosen to fit the experimentally determined profile more accurately.Based on the results of the emission intensity fit, the power density profile was modelled as: The peak of power density is set at r 0 = 1.1r 0,I , while the characteristic decay of the power density within the plasma core was set to be δ p .Besides the presence of an exponentially decaying power density profile in the core of the plasma, the concave shape of P d used in this work is different from the convectional Gaussian profile used in previous works [13,16].This choice was made to maintain consistency with the fit function used for the experimental emission intensity profiles.Furthermore, the pressures investigated in those studies were above or equal to 100 mbar, at least four times higher than the conditions investigated in this work.In a work by Lobaev et al [47] dealing with modelling of an MW chemical vapor deposition reactor it was shown that the radial power density profile tends to match the plasma optical emission profile with increasing pressure.Moreover it showed a step-like power density radial profile at ∼80 mbar.The value of P max (t) in equation ( 21) is found by imposing that: where P abs (t) is the experimentally measured absorbed power, L is the plasma length discussed in the next section and R in is the inner radius of the quartz tube.The intensity of spontaneous optical emission from a pure nitrogen plasma is mostly due to the second (N 2 (C) → N 2 (B)) and the first ((N 2 (B) → N 2 (A)) positive system [48], second positive system (SPS) and first positive system (FPS), respectively.The former is located in the wavelength region between 260 and 550 nm, while the latter between 478 and 2300 nm and their relative detected intensity depends on the radiative lifetimes of the excited states, the emitter number density and the camera quantum efficiency.
The quantum efficiency of the intensifier used in this work peaks in the wavelength region between 400 and 600 nm, sharply decreases for λ < 400 nm and slowly decreases for λ > 600 nm.Both emission bands overlap with the intensifier operating range.The wavelength dependence of the intensity of the peaks in each band is related to the population of vibrational levels of N 2 (A,B,C), which is not trivial to predict [12,49]; it is therefore assumed that the shape of the bands does not vary relevantly within the conditions analysed in this work.This allows us to consider the quantum efficiency as constant, meaning the only factor determining the emission profile would be the emitters' density.
Two main effects might cause a hollow density profile for the emitters: (i) Sharp increase in gas temperature and therefore decrease in number density in the core, as pressure remains constant.(ii) Lower density of emitters due to lower power density and therefore lower electron density, electron temperature or electronic excitation.
The temperature effect (1) can be discarded, as an increase in gas temperature is also linked to an increase in electronic excitation and therefore population of emitters, so that the entity of the temperature gradient necessary to cause a visible minimum of emission intensity in the centre of the plasma should be higher than the one measured.
Effect (2) is relevant if impact with electrons is the main creation mechanism for electronically excited N 2 .This is indeed the case, as it will appear clear in the results section.Regions where power density is lower can result in a lower degree of electronic excitation and therefore a lower emission intensity.Other reactions, like associative ionisation or dissociation, or transport phenomena can relevantly affect the emitters density profile, so that the emission profile is different from the applied power density input.
The optimal parameters for the power density profile were checked by comparing experimental and simulated optical emission, considering the FPS and the SPS as the main constituents of the broadband emission from the N 2 plasma.

Effective axial transport.
Plasma length (L) is used in the model to calculate a characteristic time for heat or mass losses due to convective transport in the axial direction, which is not resolved by the fluid equations.This parameter is also used to determine the peak power density value through normalisation.Note that overestimating the plasma length would lead to low power density input (see equation (22)), but also slower mixing with cold gas.On the other hand, shorter plasma would lead to higher power input, but faster losses.If all processes (mainly species production and gas heating) rely in the same way on transport and power density, overestimation or underestimation of L would play a small role in determining plasma parameters.However, the different timescales of gas heating, electronic excitation and vibrational excitation make axial transport and residence time play a relatively different role in the dynamics of the aforementioned quantities.
The plasma length (L) effectively determines the characteristic time of axial transport which controls both losses of particles and cooling of the gas.Its implementation follows the reasoning of Wang et al [50].The idea is that the cold pure N 2 gas entering the discharge from upstream is mixing with the hot gas within the discharge, diluting its properties.The infinitesimal change of any property φ (dφ) in the time dt is proportional to the ratio between the infinitesimal volume of 'cold' gas coming from upstream during that time (dV) and the total plasma volume (V): In this model: where r is any radius and R = r + dr.The background is considered to be pure N 2 gas at room temperature.This description assumes that any property φ is described by a step function in the axial direction, so that a clear separation between background and plasma can be drawn.In reality composition and temperature change smoothly along z and an effective length needs to be used instead.
This effective length is estimated using optical emission intensity profiles.This is an approximation as there is no single plasma property directly correlated to deposited power density.Nevertheless sharp gradients in emission intensity are correspond approximately to sharp gradients in other plasma properties.The emission profile along the axial direction is niether Gaussian nor with a clear analytic shape.The plasma length was therefore evaluated as the distance between the two diametrically opposite points at 1/2 of the peak intensity value.Note that this effective length has been estimated for each radial position.It becomes smaller with increasing radial position.
The same effective length introduced in this section was used to impose the normalisation of the power density profile, so that the total absorbed power is equal to the total forward power from the MWs.While in a time resolved measurements this could be an accurate picture, in a steady-state discharge other mechanisms like diffusion and convective axial transport might be playing a role in smearing the density of emitters over regions where the power density is low.This behaviour would make our estimate of peak power density lower than its real value.
Nevertheless in this work different regions of the plasma will be identified based on the relative efficiency of transport or power deposition in driving chemistry.

Chemistry set
The chemistry set used in this paper is the same introduced in a recent work by the same authors [17]; reaction numbers from that work are used in the present paper.The changes to such chemistry set are discussed in the following.
The set of monoquantum V-V (V-V i ) collisions of the type: has been expanded from exclusively i = 1, 2 to i = 1, .., 20.
Differently from that previous work, where gas temperature exceeded 2000 K, in this work the rate coefficients of V-V i and V-T collisions have been calculated according to the SSH (Schwartz, Slavsky, Herzfeld) theory, as described in recent works by Guerra et al [23] and Davies et al [51].
State-specific dissociation of N 2 through collisions of vibrationally excited nitrogen and gas molecules has been added to the previous chemistry set.The state specific rate coefficients used in this work are calculate according to the FHO theory [26].The corresponding recombination reactions are obtained through the principle of detailed balance.Figure 4 shows the total dissociation rate obtained by averaging the state-specific coefficients from FHO theory (used in this work) on a Boltzmann vibrational distribution function (VDF) compared to the experimental fit tabulated by Park et al [52].The rate coefficients obtained using the Fridman-Macheret approximation as in the work by Van Alphen et al [8] and Kelly et al [10], and the rate coefficients proposed by Guerra et al [23] (which were part of the chemistry set previously used [17]) are also averaged in the same way and compared in the same figure.It is clear from the figure that for a Boltzmann VDF the processes indicated as relevant for dissociation by Guerra et al [27] are dominating over dissociation of N 2 (v) due to impacts with gas molecules (as calculated in this work) for T < 4000 K. Above such temperature, these processes become relevant and necessary to match the experimental fit reported by Park et al [52].The dissociation rate calculated from Fridman-Macheret state-specific rate coefficients provides a higher average dissociation rate for T > 2000 K.The rate coefficient (k I9 (T g )) for the dissociation of N + 4 through impact with molecular nitrogen (N 2 ) (reaction (I9) in [17]): has been modified to prevent exponential growth with increasing gas temperature.This reaction was originally introduced by Gordiets et al [53] for low-pressure (2.7 mbar) and lowtemperature (T g < 1200 K) conditions.In this work the rate coefficient increases exponentially up to 1500 K and then maintains a constant value:

Electron kinetics
The complete set of reactions included in the electron kinetics is discussed in a previous paper by the authors of this work [17].In the following notable additions are discussed.The state-specific excitation of N 2 (X, v) to electronically excited states has been discussed as a relevant mechanism for nitrogen excitation in the plasma afterglow by Guerra et al [27] and was also introduced by Dyatko et al [28].Considering the high degree of vibrational excitation shown by experiments, the present model also includes electronic excitation from vibrationally excited states of nitrogen.The cross section for such processes has been calculated using the semi-empirical formula indicated by Loureiro et al [54,55], with the Franck-Condon coefficients from different papers [48,56,57].The role of these processes in providing electronic excitation will be discussed in the results section.
A recent publication by Volynets et al [21] has indicated the relevance of electron impact dissociation of vibrationally excited species through the excitation of dissociative states of N 2 (B,C and a).In this work those processes are added to the electron kinetics by adapting the cross sections for electronic excitation using the semi-empirical formula mentioned in the paragraph above.Note that in the work by Volynets et al [21] the Gryzinski approximation is used instead, as indicated by Cacciatore et al [58].These cross sections differ slightly from the ones calculated by Loureiro et al [55].
Electron-impact processes rate coefficients are calculated based on the electron energy distribution function (EEDF) obtained by the solution of the EBE.Due to its dependence on gas composition the EEDF is calculated at different radial positions; a detailed discussion on this can be found in the next section.The EBE is solved for different values of reduced electric field (E/N) so that the rate coefficient for every reaction can be tabulated as a function of electron mean energy (ϵ).
Differently from the work by Kelly et al [10], in this work superelastic collisions with electrons causing de-excitation of N(P) and N(D) to N(S) are included.This is relevant, since only recombination of N(S) is added to the model: if cascading back to the electronic ground state is hindered by the absence of de-excitation due to collisions with electrons, the recombination rate would decrease significantly.
Moreover, superelastic collisions with all excited atoms or molecules are included in the model, as they relevantly change the shape of the the electron energy distribution function, causing the rate constants for electron impact collisions to undergo variations.An example of this is shown in the work by Guerra et al [23].In this work the effect of the absence of superelastic collisions is not investigated.

Numerical implementation
Equations (1), (10) and ( 13) have been discretised on a 1D spatial grid using the finite volume method, as already applied in previous works [16,17].The resulting matrix has been inverted using the tridiagonal matrix algorithm.The spatial grid has been divided into 200 cells with a width ∼0.06 mm each.This value was found to guarantee smooth profiles for all the quantities described by the model.The fluxes are discretised on the grid using the central difference scheme.Time integration is performed explicitly.
Due to the large difference in electron and heavy particle dynamics characteristic times, two different time steps have been defined.A time step ∆t e (typically ∼10 −12 s) has been used to advance equations ( 10) and (1) for charged species.This choice has been made to ensure a correct evaluation of y e (and subsequently n e ) via equation ( 7), which is needed to obtain ϵ = n ϵ /n e .A different time step (∆t ∼ 10 −10 s) was introduced for the heat balance equation and the mass balance equations of neutral species.In both cases ∆t has been chosen so to satisfy the Courant-Friedrichs-Lewy condition, that is: where r p is the rate of a given process p; these include chemical reactions, diffusion and radial convection (for equation (1)), heating, radial convection enthalpy losses and conduction for equation (13), and electron heating (from equation ( 12)), diffusive and convective radial transport for equation (10).As discussed before, the EEDF depends on gas composition, which changes both in time and in space.In order to have a radially accurate evaluation of the electron distribution function, the Boltzmann equation has been solved using BOLSIG+ [33] at designated nodal points on the 1D grid.The positions have been chosen in order to better resolve the regions with sharper gradients, that is the region between 3 and 6 mm from the centre of the tube, as revealed by measurements of optical emission.The designated nodal point do not change with time.The value of all rate coefficients and transport parameters on all other cells in the grid has been calculated via linear interpolation.The Boltzmann solver is called whenever at least one of the following conditions is satisfied: • The electron mean energy falls outside the range covered by the previous iteration.• The maximum relative change in the mass fraction of any species at any of the radial positions exceeds 0.1.
Whenever one of the above conditions is satisfied the EEDF is updated at every position in order to prevent sharp discontinuities in the calculated rate coefficients which are otherwise detected.
The reduced electric field values used as input by the Boltzmann solver are bound to be in the interval [0.8E/N(ϵ min ), 1.2E/N(ϵ max )], where ϵ min (ϵ max ) is the absolute minimum (maximum) value of the electron mean energy in the domain.The actual values of E/N are shown in the results section.
A total of five pulses (frequency of 30 Hz) are simulated in this work, since a periodic steady-state is observed after and including the third pulse.The total simulated physical time is 135 ms (the low input power phase of the last pulse is stopped after 3 ms).Using 16 CPUs the total wall time of the simulation is ∼2.9 × 10 5 s (∼3 days and 8 h).The simulation of the first 200 µs of the pulse is the most computationally expensive, requiring ∼30 s of real time for 1 µs of physical time.After the initial power pulse, the computational cost decreases to ∼1.5 s for 1 µs of physical time.

Results
Figures 5(a) and (b) show the experimentally detected emission intensity and the one obtained from the self-consistent simulations, that is the sum of the emission intensity from the SPS and the FPS described in the methods section.The two profiles compare well, which corroborates the choice of input power density profile, which is shown in figure 5(d).
The decrease in power density in the core of the plasma as a function of time (estimated from experiments according to equations ( 20) and ( 21)) is sharper than the decrease in emission intensity both measured and simulated in that region, making it clear that the two are not in a linear proportionality relation.In fact, the population of the two triplet states involved in the transitions for the second and first positive system (N 2 (B) and N 2 (C)) relies on electron impact processes, whose rates are proportional to reduced electric field, electron temperature and electron density.The latter, as shown in figure 5(c) does not relevantly decrease in the core.Moreover, destruction mechanisms and transport of the aforementioned excited species play a role in shaping the emission profile obtained by the simulations.
The input power density profile chosen in this model, therefore, provides a good agreement with the experiments within our framework, that is with the chosen chemical and vibrational kinetic set, electron kinetics, heating mechanisms and transport parameters.Cross checking of other quantities like temperatures and electron density, illustrated in what follows, will strengthen the validity of the model.
Figure 6 shows the temporal evolution of electron temperature and density, both measured and simulated.Note that, in order to avoid numerical instability in the model, the input power density increase starts at −1 µs and P d peaks at 0 µs, anticipating the rise of T e with respect to the experimental   measurement; this artefact does not influence the remainder of the discharge neither does it invalidate the obtained trends.Moreover, the discontinuities in the simulated temporal evolution of T e are caused by the fact that the EBE is not solved at every time step.After the 200 µs pulse peak of power density shifts back to the axis of the tube (r = 0 mm), explaining the increase in T e immediately after t = 200 µs in the core; similarly to what happens during the pulse, the discontinuities afterwards are due to the discrete updates of the EBE, which during the post discharge are less frequent.
Electrons are the only species to receive energy from the applied electric field.As such the electron temperature has an initial sharp increase around 0 µs and reaches ∼2.5 eV in less than 1 µs.This is shown both by experiments and simulations in figure 6(a).At this point in time, when electron density is still low and power density has its maximum value, the electron mean energy density increases sharply due to the increase in deposited power density P d (r, t), which is the only relevant source term in equation ( 10), while the peak in electron density is reached only after ∼5 µs (see figure 6(b)) at r = 0 mm.This initial rise in electron density in the core is caused by electronimpact ionisation of N 2 , as shown in figure 7(a).Soon after this initial rise, associative ionisation (blue full lines in figure 7) of excited N 2 through reactions: becomes the dominant process in all regions of the discharge, while associative ionisation involving excited atomic nitrogen (see full orange lines in figure 7) becomes a relevant contribution outside the core after ∼100 µs.The relevance of this processes for ionisation effectively lowers the ionisation potential for N 2 and decouples the temporal evolution of electron density from the temporal evolution of the applied power density.
Electron-ion recombination (dashed lines in figure 7) mainly consists in dissociative recombination of either N + 4 and N + 2 .The former is dominant at all positions at the start of the discharge; as the rate coefficient for reaction (I9) increases exponentially with temperature during the time of the pulse, N + 4 is dissociated to form N + 2 , which becomes the dominant ion and eventually recombines through reactions (I11) and (I12).
The gradual loss of electrons from the core of the plasma during the duration of the pulse was also observed by Kelly et al [10] for a similar discharge and by Shkurenkov et al [59] for a pulsed nanosecond discharge at 133 mbar.Both gas heating and ionisation-recombination kinetics are responsible for this behaviour.
The rising gas temperature (figure 12) is partially responsible for the decreases in the absolute number density of electrons in the core.After ∼50 µs the molar fraction of electrons starts decreasing, indicating that the increasing gas temperature is not the only contributor to the low electron density reached at the end of the pulse.
The other cause of this behaviour the plummeting in electron temperature following the increased electron density, which in turn favours dissociative recombination over associative ionisation.This decrease in electron temperature also lowers the rate coefficient for electron-impact ionisation of one order of magnitude with respect to the initial value.
Radial transport contributes relevantly to the increase in electron density in the region of the plasma for r > 7 mm.This shows that electron density is mainly determined by chemistry and by the competition between ionisation and recombination mechanisms, rather than transport.
After the first 5 µs the electron temperature is maintained almost constant through the duration of the pulse at all radial positions.The relative change of T e from the value from 5 µs to 200 µs is of ∼10% at 0, 3 and 9 mm and ∼20% at 6 mm.The almost constant T e in the core is a symptom of electron energy density decreasing at a similar rate of electron density.The temporal evolution of the electron energy density is mostly determined by power density input (source) and inelastic collisions (sink).Excluding the outer position (9 mm), the reduced rate of energy losses due to inelastic collisions ( P inel ngas (ϵ) equation ( 10)) is almost constant.Therefore reduction of the energy sink term in equation ( 10) can only happen if heavy particle density decreases (in other words T g increases) or electron density decreases.The timescale of both these processes is longer than the timescale of the decrease in power density input, causing a net decrease in n ϵ .Similar measurements in CH 4 performed by Butterworth et al [60] show a much smoother decrease in electron temperature accompanied by a smooth increase in electron density; this could be linked to the presence of slower ionisation mechanisms in CH 4 or slower gas heating.
The electron density profile is almost flat in the core (r < 4 mm) of the plasma (figure 5(c)) despite the input power density being hollow.This is due to the fact that the main ionisation mechanism involves electronically excited molecules.The efficiency of formation of these species depends on E/N.Though the model does not determine the reduced electric field self-consistently, it is possible to estimate it through interpolation knowing the electron temperature T e ; both of these quantities and the electron density profile are shown in figure 8. Due to the increase in gas temperature in the core, E/N presents an almost flat profile in this region, maintaining the production rate of electronically excited species high enough to provide ionisation.More details on this will be given in the final section, where formation of electronically excited species is analysed.

Vibrational excitation
In this section the experimental and simulated vibrational temperatures are compared.The experimental values of T v,01 rely exclusively on the intensity of the first two peaks of the Q-branch of the Stokes Raman emission, whereas T v,15 is based on the intensities of the peaks from 1 to 5 and calculated with a Boltzmann fit, making it a more reliable estimate.Moreover, the vibrational temperature based on the population of levels from 1 to 5 is representative of the mean energy stored in vibrations of N 2 ; for these reasons only T v,15 is presented in this section.Figure 9 shows the evolution in time of the vibrational temperature T v,15 for different radial positions.Note that the simulated value of T v at the beginning of the pulse (0 µs) is ∼1000 K higher than the measured one.This could be due to the fact that, when the power is turned on, the plasma moves significantly making the first measurement subject to more uncertainties.Moreover, long-term vibrational cooling is underestimated by the 1D model as it fails to take into account turbulent effects.The experimental vibrational temperature for radial positions r ⩽ 6 mm is well reproduced by the simulation.In particular the different timescales for vibrational heating of the core and the outer positions are well matched by the computed T v .
In the first tens of microseconds (t < 20 µs) of the discharge, vibrational energy loading due to impact of electrons on heavy molecules dominates.The rate coefficients for these processes depend on electron mean energy and density, which in turn depend on how much power is given as input to the system.It has been shown above that electron temperature and density increase sharply in this period of time due to the peak in power density input, while transport phenomena are negligible.
After 20 µs the electron density in the core decreases along with the electron mean energy, leading to slower e-V processes and to slowing down of the vibrational temperature temporal evolution; moreover, as the vibrational temperature increases so does the relative vibrational energy transferred to electronic excitation (V-E), which becomes the main loss mechanism for vibrational energy.
According to the simulation the vibrational temperature at r = 0 mm reaches a maximum at around 50 µs, while this happens at 100 µs at r = 3 mm.T v in outer regions slowly increases until the end of the pulse.This change in trend is caused by transport becoming the leading mechanism for vibrational heating of the outer layers of the discharge.Figure 10 shows the relative contribution of axial and radial transport and chemical reactions to the losses and gains of vibrational energy: between 4 and 6 mm from the centre diffusive radial transport is in competition with chemical kinetics.Note that this is the same region where the sharp gradient in optical emission and therefore input power density is located.
In the region closer to the wall (9 mm) the vibrational temperature is overestimated by ∼1000 K.According to figure 10 in this region the vibrational heating balance is dominated by diffusive transport and axial convective transport, making it extremely sensitive to assumptions on the flow pattern and gas dynamics.
As a result of different mechanisms dominating vibrational heating and cooling, the shape of the VDF as a function of time depends on the radial position investigated (see figure 11).The dominance of electron-impact excitation over V-V collisions in the core causes the evolution of the VDF towards a Boltzmann-like distribution with a depleted tail due to dissociation (figure 11(a)).Differently, at 9 mm (figure 11(b)) the VDF shows an evident overpopulation of higher vibrational levels (with respect to a Boltzmann distribution).This is a sign of both slow electron impact processes and V-T quenching.

Gas temperature
Figure 12 shows a comparison between the temporal evolution of T g measured in experiments and the one obtained from the simulations.Figure 12 shows that both experiments  and simulations highlight different timescales for gas heating depending on the radial position.Namely, heating in the core is faster than at outer positions.This behaviour is explained by the faster production of electronically excited species in the core due to the higher initial power density.The gas heating of the core is generally in agreement with the experimental measurements within experimental error.Gas temperature at 3 mm is overestimated by ∼250 K for the first half of the discharge and then underestimated up to ∼500 K during the second half, revealing a different heating dynamics with respect to the one modelled.The measurements at 6 mm reveal a gas temperature reaching ∼1000 K at the end of the pulse, while simulated results increase up to 750 K at the same position.
We show in figure 13 the relative contributions of transport and chemistry to the total gas heating or cooling as a function of both radial position and time.For r < 4 mm radial diffusive transport represents ∼75% of the gas cooling, while for the outer region (r > 8 mm) it is the main heating mechanism.Heat balance in the region between ∼4 and ∼8 mm is governed by an interplay of transport and chemistry.Exothermic reactions and axial transport are the main heating and cooling mechanism, respectively, in all other cases.
Experimental measurements at 3 mm show a fast gas heating starting at ∼60 µs, which is not reproduced by the simulation.In fact, the gas heating predicted by the model in this region is mostly determined by exothermic chemical reactions, while conduction is the main cooling process.This causes the temporal evolution of T g to proceed like in the core.This dynamics changes for the region at 6 mm, where conduction is a heating mechanism and cooling only comes from axial convective transport.This dynamics is much much more similar to the one showed by the experiments at 3 mm; the position of the crossing between these two different dynamics is placed where the power density input peaks.Uncertainty in this quantity can justify the mismatch between model and experiments.
The rather low gas temperature predicted by the model at 6 mm can be explained with an underestimation of the exothermic reactions heating rate.In fact, despite being the most important cooling mechanism for r > 4 mm, the absolute value of the cooling rate due to axial convection at 6 mm is ∼2 orders of magnitude lower than the heating rate in the same region (figure 14).The slow heating of the region around 6 mm is therefore most likely due to an underestimation of exothermic reaction rate.The comparison of different heating mechanisms shown in figure 14 shows that quenching of the state N 2 (B) is dominating over the whole plasma volume, meaning that underestimation of N 2 (B) population might be the cause of the wrong estimated heating rate.
Another way to increase heating in the 6 mm region could be adding an effective turbulent conductivity coefficient to the already computed thermal conductivity λ.This approach would smooth out the peak of gas temperature causing the heating of the core to become slower and therefore would deteriorate the comparison with experiments in that region.The turbulent effective conductivity used in works simulating CO 2 [13,16,22] was specifically calculated for those discharges; the implementation of such term to simulate an N 2 discharge would introduce an additional approximation.
Gas heating depends on electron density more loosely than T v as molecules need to be excited first and then quenched to heat the gas.The space resolved measurements by Gatti et al [7] showed that the high non-equilibrium outside the core is caused by the lower gas temperature in this region as opposed to the high vibrational temperature.This effect is not due to lower V-T quenching rate, as this is not the main heating mechanism outside the plasma core, but rather to the lower electronic quenching rate.
The quenching of N 2 (B) was identified as the first dominant mechanism for heating of pure N 2 nanosecond discharges [59,61] at 133 mbar.In those works, this mechanism is soon overtaken by quenching of atomic nitrogen metastables, which in those conditions are quickly formed due to the high electron density (∼10 20 m −3 ) reached thanks to the high initial applied reduced electric field (∼350 Td, as opposed to the ∼125 Td of this work).The higher contribution of N 2 (B) quenching over pooling reactions (N * 2 + N * 2 → N * 2 + N 2 in figure 14) was also highlighted by Pintassilgo and Guerra [42] at 2.6 mbar, where the heating rate of such reaction was shown to be dominant under the hypothesis that all the available energy (1.18 eV) is used for heating.Finally, quenching of N 2 (B) is also identified as the main heat source at atmospheric pressure in a constricted glow discharge [62].
Furthermore, energy transfer from vibrations to heat through V-T quenching with molecular and atomic nitrogen and non-resonant V-V collisions provides ∼20% of the total heat at the later stages of the power pulse, with V-T (N 2 ) being the dominant process; it is interesting to notice that the heating rate of these processes peaks in the centre, where T g is higher and therefore oppose the tendency of reaction (N7) [17] to broaden the gas temperature profile.The role of V-V collisions in gas heating was highlighted by Pintassilgo and Guerra [42] and by Volynets et al [21] for glow discharges at pressures between 2.6 and 53 mbar.In our work, the relative net contribution of V-V process to gas heating is about 2% during the later stages of the discharge.The work of Pintassilgo and Guerra [42] shows a much more relevant contribution of these processes, but on a longer timescale (∼10 ms).Moreover, in that work the total heating rate is two orders of magnitude lower than the one shown in this work, which leads to a gas temperature T g < 400 K and to the formation of an evident plateau in the VDF, populating the levels involved in V-V collisions.If only the exothermic component of the contribution to the heat equation provided by V-V collisions (instead of the net contribution, represented in figure 14) is considered, such processes make up ∼20% of the total heating, which is in accord with the estimated values of Volynets et al [21].In none of these studies V-T collisions with N 2 is mentioned as a relevant contribution to gas heating, differently from what is seen in this work.The main difference in our case is the higher electron density (∼1 order of magnitude).At these values of n e the VDF is mainly shaped by e-V collisions, which allow the tail of the VDF, where V-T collisions are most effective, to be more populated [63].
Previous modelling of a similar discharge by Kelly et al [10] showed a greater contribution to heating from V-T collisions with atomic nitrogen (N) than the one predicted in this work.This was due to the significantly higher peak power density applied in that simulation (four times higher than the one estimated in this work).The higher power density would enhance electron-impact dissociation, which is the leading cause of N formation and lead to a higher dissociated fraction, therefore enhancing V-T (N) rates.Due to the dominance of electron impact collisions, the difference between the chemistry set used by Kelly et al [10] and the one used in this work in terms of dissociation rates is minimal.Nevertheless, the set of reactions used by the former is used to show the impact of these changes on the final results.Moreover, results obtained by implementing scenario B from the work by Volynets et al [21] are calculated and compared.The most significant changes can be observed for the gas temperature temporal profiles, which are compared in figure 15.Compared to the model employed by Kelly et al the gas temperature calculated in this work differs negligibly, due to the limited effect of N on gas heating and the fact that the set of exothermic reactions that provide heating is the same.
The model suggested by Volynets et al underestimates the gas temperature by ∼1000 K for t > 50 µs.These authors excluded from the chemistry set reaction (N8) and lowered the enthalpy contribution of reaction (N7), which are the main sources of heat in the conditions investigated in this work (orange line in figure 14).This choice was made in that work to match the degree of dissociation observed in experiments and lower the otherwise too high gas temperature.That model was developed to simulate a glow discharge, where axial transport and cooling can be neglected, due to low flow rate.In the case of our MW discharge, axial cooling is the main cooling mechanism of the edge of the plasma and is responsible for ∼30% of the cooling in the core (figure 13).

Electronic excitation of N 2
As discussed in the previous sections, electronically excited species in pure N 2 plasmas are responsible for ionisation and gas heating.In turn, electron creation influences the vibrational heating rate, making electronically excited species vital also for this mechanism.In this section, the formation and destruction of the main excited species is studied to elucidate the relative role of chemistry and transport.
Figures 16-18 show the formation and destruction of electronically excited species; similarly to the case of vibrational excitation and gas heating, it is possible to draw a separation between the region at r < 5 mm (core) and the outer region where radial diffusive transport provides a relevant contribution to either creation or destruction of excited states.
In the core the three triplet states are strongly coupled by electron impact excitation and de-excitation and heavy particle impact de-excitation.In particular the population of the states N 2 (A) and (B) is mostly determined by the balance of the reaction: N 2 (B) radiative de-excitation to N 2 (A) accounts only for ∼7% of its losses in the core and its rate is always slower than the one of reaction (N7) (figure 17).Moreover, it is overrun by  heavy particle impact de-excitation to the ground state (N8) in the outer regions.The emission from the FPS (N 2 (B)→N 2 (A) + hν) will therefore be suppressed at r > 5 mm.The excited state N 2 (C) is for ∼80% formed through collisions involving state (A).Radiative de-excitation is the main (∼99%) cause of loss of N 2 (C), meaning that the emission intensity profile for the SPS should match the profile of formation of N 2 (C).
The radial profiles of the number density of these three states is similar due to the fast coupling reactions.If these reactions are ignored, from figures 16-18 one can see that the reactions pumping energy into the triples states are: N 2 (v) + e → N 2 (A, B, C) + e (31) while depopulation happens mainly due to: N 2 (A, B, C) + e → N 2 (v) + e (32) N 2 (B) + N 2 → N 2 + N 2 (33) and diffusive transport.Thus the actual shape of the number density profile of molecules excited to the low-lying triplet states is determined by the balance between the electronimpact excitation reactions causing energy pumping into the metastable N 2 (A), which depend on and are enhanced by the relatively high vibrational temperature, and the corresponding superelastic collisions.Impact between N 2 (B) and ground state N 2 is responsible for ∼15% of the losses of energy from the triplet system and tends to push the population of N 2 (B) towards the thermodynamic equilibrium prediction.Those collisions are particularly effective at low temperature.Diffusive transport of both N 2 (A) and N 2 (B) contributes to ∼9% to the loss of triplet states in the region between ∼2 and ∼5 mm, which smooths the peak that would otherwise grow due to the higher power deposition.
Figure 19 shows the relative amount of electron energy used for different processes at 200 µs.The low rate of electronimpact electronic excitation in the outer region is due to the decreasing reduced electric field (see figure 8).This causes the slow heating rate in this region, as shown in previous sections due to the absence of N 2 (B).

Conclusions
A 1D fluid model with radial resolution and time resolution has been developed to simulate a pulsed MW plasma discharge in pure N 2 .The validity of the model was checked through a thorough comparison with temporally resolved measurements of electron temperature and density in the plasma core, and spatially and temporally resolved measurements of gas and vibrational temperature, and spontaneous broadband optical emission.The latter has been used to validate the choice of power density input profile within the framework (chemistry set and transport parameters) of this work; qualitative and quantitative agreement between simulation and experiments has been reached within ∼15%.Simulated and measured electron parameters agree qualitatively and quantitatively within ∼20%.Measurements in the core of gas and vibrational temperature agree with simulations qualitatively and quantitatively within 15%, with very few exceptions.Between 3 and 6 mm, agreement is within 20%, also accounting for experimental error.Vibrational temperature in the outer region is overestimated of about 50%.The source of uncertainty in this region is identified to be in the treatment of diffusion and/or axial transport.
The optical emission intensity radial profile appears narrower than the imposed power density profile, making it an unreliable means to determine the main input parameter for modelling of MW plasmas.This is due to the lower fraction of electron energy directed to electronic excitation in the outer region as the decreasing gas temperature with radial position in this region causes a sharp decrease in reduced electric field.The radial decrease in power density (∼50%) from the peak towards the core is greater than the decrease in optical emission intensity (∼25%).On the one hand, this is due to diffusive transport causing electron and excited species losses from the peak region.On the other hand, the dominance of associative ionisation for electron formation in the core also contributes to the decoupling of electron number density from the applied power.
The development of a more precise collisional radiative model is advisable; the inclusion of escape factors or the formulation of a self-consistent model for radiative deexcitation and self-absorption of photons are among possible improvements.
Mass and energy diffusive radial transport decouples the plasma behaviour from the applied power evolution in time and space.It has been shown that radial and axial transport are in competition with chemical reactions for the formation and destruction of electrons, vibrationally excited species and the heating of the gas at r ⩾ 4 mm.This is important as it changes the residence time of molecules transiting those regions and their specific energy input.
Further decoupling of plasma temporal behaviour from power density is caused by associative ionisation being the main ionisation mechanism in the plasma core.The slow decay of electron density in this region is caused by the reduction of associative ionisation sustained by the long lived metastable N 2 (A) which is formed at the start of the discharge.The excited states involved in associative ionisation are mostly formed by electron impact.
Finally, the value of the reduced electric field (E/N), rather than the power density, determines energy partition of energy between electronic excitation and vibrations.This affects the observed emission intensity and the heating of the plasma.The low E/N at r > 4 mm and therefore slow electronic excitation explains the low temperature kept in this region and hence the long life of vibrationally excited species that are diffusing there from the core.This explains the non-equilibrium previously observed in CW discharges in pure N 2 .
Axial convective transport always represents the main loss mechanism for gas heating and vibrational heating, while it is only a minor sink of electronically excited molecules and electrons.As such, any underestimation or overestimation of the particle velocity will result in a significant decrease of vibrational and gas heating.

Figure 1 .
Figure 1.Measured forward (blue line) and reflected (red line) MW power, and resulting absorbed (green line) power.

Figure 2 .
Figure 2. Schematics of the experimental setup with the optics used for Thomson scattering measurements.Reprinted with permission from [9].Copyright (2022) American Chemical Society [9].

Figure 3 .
Figure 3. Examples of experimentally measured spectra (black line) compared to their best fit (red line).The two examples are chosen to show the highest signal-to-noise ratio recorded in the measurements.(a) Thomson scattering spectrum recorded in the core at 200 µs; the resulting electron temperature and density are, respectively (0.76 ± 0.08) eV and (1.9 ± 0.4) × 10 18 m −3 .(b) Vibrational Raman scattering spectrum recorded in the core at 200 µs; the resulting vibrational temperatures are T v,10 = (4.5 ± 0.6) × 10 3 K and T v,15 = (6.8± 0.4) × 10 3 K.

Figure 4 .
Figure 4. Dissociation rate coefficient as a function of gas temperature obtained by averaging different state-specific dissociation reaction sets (see text), compared to the dissociation rate obtained by fitting experimental data [52].

Figure 5 .
Figure 5. (a) Experimental optical emission intensity profile normalised by its maximum value.(b) Simulated optical emission (from SPS and FPS) intensity profile normalised by its maximum value.(c) Simulated electron density profile.(d) Power density input profile.

Figure 6 .
Figure 6.Simulated and experimental electron temperature (a) and electron density (b) temporal evolution.Simulated results are presented at four different distances from the axis of the tube, while experimental values are only obtained at 0 mm.The shaded area includes values calculated at ±0.5 mm around the selected position.

Figure 8 .
Figure 8. Radial profiles of reduced electric field E/N and normalised electron temperature Te and density ne at three different time: 0 µs (a), 100 µs (b) and 200 µs (c).

Figure 9 .
Figure 9. Simulated and experimental temporal evolution of the vibrational temperature T v,15 at four radial positions.The shaded area includes values calculated at ±0.5 mm around the selected position.

Figure 10 .
Figure 10.Radial profiles of main vibrational activation (solid lines) and deactivation (dashed lines) mechanisms at 0, 100 and 200 µs.The abbreviations V-E and V-D indicate energy gain or losses due to electronic excitation or dissociation from vibrational levels, respectively.

Figure 11 .
Figure 11.Simulated vibrational distribution functions of the ground state of N 2 (a) in the core and (b) at 9 mm from the centre of the tube.Different colour indicates a different time.

Figure 12 .
Figure 12.Simulated and experimental temporal evolution of the gas temperature Tg at four radial positions.The shaded area includes values calculated at ±0.5 mm around the selected position.

Figure 13 .
Figure 13.Composition of sources and sinks of translational energy as a function of radial position at different times.

Figure 15 .
Figure 15.Temporal evolution of the gas temperature in the plasma core (r = 0 mm) obtained with three different chemistry sets (see text).

Figure 16 .
Figure 16.Radial profiles of the main creation (right) and destruction (left) mechanisms for N 2 (A) at 200 µs.

Figure 17 .
Figure 17.Radial profiles of the main creation (right) and destruction (left) mechanisms for N 2 (B) at 200 µs.

Figure 18 .
Figure 18.Radial profiles of the main creation (right) and destruction (left) mechanisms for N 2 (C) at 200 µs.

Figure 19 .
Figure 19.Radial profile of the relative contribution of different processes to the losses of electron energy at 200 µs.
In this work, the Lennard-Jones coefficients used to calculate the D ij 's of electronically excited species and vibrationally excited species are considered equal to the ones used for the ground state for lack of more precise data.The absence of a turbulent term, differ- ently from what has been done byViegas et al [13]and Vialetto et al[16]for CO 2 is mostly related to the lack of data regarding N 2 MW discharges; its impact on the simulation results is briefly discussed in the Results section.Finally, quasi-neutrality is guaranteed by imposing: