Charged particle kinetics and gas heating in CO2 microwave plasma contraction: comparisons of simulations and experiments

This work investigates kinetics and transport of CO 2 microwave plasmas through simulation results from a 1D radial fluid model and experiments. Simulation results are validated against spatially resolved measurements of neutral species mole fractions, gas temperature, electron number density and temperature obtained by means of Thomson and Raman scattering diagnostics, yielding good agreement. As such, the model is used to complement experiments and assess the main chemical reactions, mass and energy transport in diffuse and contracted plasma regimes. From model results, it is found that, as pressure is raised, the inhomogeneous gas heating induces significant gradients in neutral and charged species mole fractions profiles. Moreover, the transition from diffuse to contracted plasma is accompanied by a change in the dominant charged species, which favours electron–ion recombination over dissociative attachment. Associative ionization rates increase in the plasma core from diffuse to contracted regime. These processes contribute to the increase in the peak electron number density with pressure, that determines radial plasma contraction.


Introduction
In the past few years, a lot of attention has been devoted to the study of CO 2 dissociation by means of microwave (MW) discharges [1][2][3][4][5].Most of this work has been motivated by high energy efficiencies, up to 80%, that have been reported in experimental investigations with subsonic vortex-confined MW reactors by former Soviet Union scientists [6][7][8][9].These high energy efficiencies have been mainly explained by preferential mechanisms leading to the onset of non-equilibrium between vibrational and translational degrees of freedom of the CO 2 molecules, favouring CO 2 dissociation at low gas temperatures.Despite numerous efforts by researchers around the world, these high energy efficiencies have never been reproduced.In fact, in plasma-assisted CO 2 conversion there is an increasingly complicated interplay of many effects, namely transport of reactive flows, complex chemistry with several species and electromagnetic fields.The interplay of these effects is not well known, but it is yet important for optimization of the reactor performances [10].
Advanced laser scattering diagnostics have been employed to obtain in situ measurements of rotational and vibrational temperatures in CO 2 MW discharges [11][12][13].It has been observed that non-equilibrium between different modes of vibration of the CO 2 molecule can be obtained by modulation of the MW power.Nevertheless, these studies suggest that, in both pulsed and continuous-wave (CW) experiments, CO 2 dissociation does not preferentially exploits non-equilibrium effects, but rather the high gas temperature in the plasma core, ranging from 3000 K to 7000 K [14][15][16].In these conditions, the complex chemistry of CO 2 MW plasmas is typically studied by means of 0D models [17][18][19][20][21].In particular, the study by Pietanza et al [22] suggests that CO and O can be described by a thermodynamic approach, while other minority species present large deviations from thermodynamic equilibrium.Similarly, in the work by Viegas et al [23], a chemical kinetics scheme for the study of CO 2 dissociation in MW discharges at intermediate to high pressures (from 60 mbar to 300 mbar) has been proposed.In [23], it has been observed that mass transport is fundamental to describe the composition in the core of a vortex-stabilized MW reactor [14].It is generally found that 0D models can qualitatively explain the behaviour of the discharge core, but discrepancies between calculations and recent experimental measurements are still present.For example, calculations by Pietanza et al [22] at 100 mbar and 250 mbar predict a neutral species composition that is largely different than the one obtained by means of Raman scattering [10,24], with C-atom mole fractions higher than in the measurements.Moreover, for the same MW input power, both models in [22,23] predict an increase in the electron temperature with increasing pressure, that is in contrast with recent experimental measurements obtained via Thomson scattering at DIFFER [25].These discrepancies suggest that further studies are needed to accurately describe the main electron-related properties of the discharge, such as electron temperature and electron number density.Indeed, chemical kinetics schemes employed in these models have not been extensively validated against experimental measurements of discharge parameters, such as neutral and charged species mole fractions, electron and gas temperatures.This is mainly due to the lack of data or to strong assumptions that are implicit in 0D model calculations that limit the applicability of these models to spatially homogeneous conditions.
In order to describe the spatial inhomogeneities of discharge parameters in CO 2 discharges, 1D or 2D fluid models are typically employed.For example, Kotov and Koelman [26] developed a 1D steady-state axial plug flow model for a MW plasma reactor.Nevertheless, the chemistry set used in [26] has been originally employed for studies of dielectric-barrier discharges (DBDs) [19] and not for MW discharges, that typically exhibit very different power densities and gas temperatures.A similar reaction scheme has been used by Ponduri et al [27] and by Wang et al [28] in a 1D axial fluid model to study CO 2 dissociation in a geometrically symmetric DBD and in a gliding arc (GA) discharge.In particular, in [28], it has been found that CO 2 dissociation in GA discharges occurs mainly through collisions with oxygen atoms and exploiting the non-equilibrium vibrational distribution function (VDF) of CO 2 molecules.However, a detailed validation of that model has not been performed due to the lack of experimental data for GA plasmas in literature.CO 2 dissociation in MW plasmas has been studied by Sun et al [29] by means of a 1D axial fluid model.They found that there is little difference between gas, vibrational and electron temperatures in the plasma, thus CO 2 dissociation in MW discharges is mostly driven by thermal chemistry.The same conclusions about thermal dissociation in CO 2 MW plasmas have been drawn by Wolf et al [30], where a 2D fluid model of the reactor used at DIFFER has been employed.However, the model developed by Wolf et al [30] considers the plasma simply as a heat source.Hence, it neglects the role of charged particle kinetics on CO 2 dissociation and on the radial contraction of the plasma taking place between moderate and atmospheric pressures.
The discharge contraction plays an important role in CO 2 conversion.In fact, the optimal condition for CO 2 conversion is achieved concurrently with this phenomenon [15,16].The mechanisms triggering contraction in molecular plasmas are not fully understood and spatially-resolved (fluid) models are important tools to shed light on them.Specifically, based on recent studies on atomic plasmas, two different mechanisms have been identified as responsible for plasma contraction: (i) electron-electron collisions, or (ii) non-uniform gas heating along the tube radius.The role of electron-electron collisions has been investigated by Golubovskii et al [31] and by Petrov and Ferreira [32], with the conclusion that these collisions have an important influence on the electron energy distribution function (EEDF) and, in turn, on increasing electron impact ionization rate coefficients in the plasma core and not on the edges.The non-uniform gas heating along the discharge radius has been considered by Martinez et al [33], by Ridenti et al [34], and by Moisan and Pellettier [35] as the main mechanism driving contraction in Ar, Ne and N 2 MW plasmas.In particular, in [33], it has been hypothesized that radial gradients of gas temperature have an influence on the reduced electric field and on decreasing charged particle losses through dissociative recombination in the core where atomic ions emerge, but not on the edges.Nevertheless, a satisfactory interpretation of the contraction phenomenon in CO 2 MW plasmas is still lacking since, in particular, plasma chemistry models employed so far have not been extensively validated against experiments.Schematic of the combined MW plasma and laser scattering setup [24].Reprinted with permission from [25].Copyright (2022) American Chemical Society.
In this work, a 1D radial fluid model is developed to study the contraction of CO 2 MW discharges.Moreover, results of the model are validated against spatially-resolved measurements of the main neutral species mole fractions, electron number density, gas and electron temperatures obtained by means of laser scattering diagnostics.The present study aims at extending the detailed investigation by Viegas et al [23], where a 0D plasma chemistry model has been developed for the study of contraction in CO 2 MW plasmas.In this work, the radial profiles of the aforementioned quantities and their relation with discharge contraction are investigated.
The paper is organized as follows.Section 2 summarizes the experimental conditions considered in this work.Section 3 presents a 1D radial fluid model coupled with a Monte Carlo flux (MCF) code for electron kinetics [36,37], with particular emphasis on the assumptions and input parameters that have been used in the model.Results of the model are compared with experimental measurements of the main neutral species mole fractions and gas temperature in subsection 4.1, and with electron number density and temperature in subsection 4.2, for a MW discharge operated at 700 W and 1000 W input MW power, 10 slm gas flow rate, and pressures from 100 mbar to 400 mbar.Insights into CO 2 MW plasma contraction are illustrated in subsection 4.3.Subsection 4.4 investigates the effect of uncertainties in parameters, such as the input MW power density and chemical rate coefficients, on the model results.Conclusions and perspectives for future work are illustrated in section 5.

Experimental conditions
Figure 1 shows a schematic of the setup in which experimental data are acquired.The experimental setup has been described in detail by Wolf et al [14,15], Viegas et al [23], and van de Steeg et al [24] and it consists of a vortex-stabilized MW reactor, with a 2.45 GHz MW source combined with laser scattering diagnostics.A pure CO 2 gas is injected tangentially into a quartz tube of 27 mm inner diameter that is placed inside a MW cavity consisting of rectangular waveguides.Gas flow rates (Φ) of 10 slm and pressures (p g ) in the range 100-400 mbar are considered in this work.MWs are supplied to the system by a magnetron operating at 700 W or 1000 W input power.The impedance of the system is optimized to minimize reflected power.The laser diagnostics allows one for simultaneous and spatially resolved measurements of rotational (gas) and electron temperatures as well as mole fractions of CO 2 , CO, O 2 , and O, and electron number densities [24,25].
Figure 2 shows the intensity distribution of the 777 nm spectral line emission of atomic oxygen (O(3p 5 P) → O(3s 5 S)) recorded for a CO 2 MW plasma at 100 mbar and 250 mbar, a gas flow rate of 10 slm and MW input power of 1000 W.
In particular, a broad intensity distribution in the radial and axial directions has been recorded at 100 mbar, that is typical of a diffuse plasma regime, whereas, at 250 mbar, the intensity distribution is narrower in the radial direction and more elongated in the axial direction, that is typical of a contracted plasma regime [15].The relation between intensity distribution and electron density profile in CO 2 MW plasmas has been investigated by Viegas et al [38].
In this work, two different sets of input parameters for the model are considered.The two sets of data are acquired from the same experimental setup that is shown in figure 1.The first set is summarized in table 1 and it is used for comparison of model results with measurements of neutral species mole fractions and gas temperature obtained with spontaneous Raman scattering.The application of Raman scattering is described in [24] and these conditions correspond to a total input MW power of 700 W. Overview of input parameters used in the fluid model for a CO 2 MW discharge at a total input MW power of 700 W. The values of the parameters have been obtained by the experiment described in [24].Overview of input parameters used in the fluid model for a CO 2 MW discharge at a total input MW power of 1000 W. The values of the parameters have been obtained by the experiment described in [25].Specifically, the peak power density at the plasma center (P in (0)) has been obtained by normalization of the emission intensity profile to the total input MW power [14].The assumptions underlying this derivation are discussed by Viegas et al [38].The width of power density profile (W ) and the length of the plasma (L) are defined as the full-width at half maximum (FWHM) of the spontaneous emission intensity profiles of the atomic oxygen triplet lines in radial and axial direction, respectively [14].
The second set of input parameters is summarized in table 2 and is used for the comparison of model results with measurements of electron number density and electron temperature obtained with Thomson scattering.The use of Thomson scattering is described in the work by van de Steeg et al [25] and this set corresponds to a total input MW power of 1000 W.
The reason behind the use of two different sets of input parameters, one for 700 W (table 1) and the other for 1000 W (table 2), is twofold.First, at 1000 W, the Raman signature in the spectra gets perturbed by the presence of C 2 emission [24,39], while the same emission has not been recorded for the 700 W conditions. Second, at 1000 W, the Thomson signal can still be fitted over the C 2 emission.Moreover, the 1000 W case is preferred over the 700 W for Thomson scattering measurements also because of the higher electron number density peak that is obtained as a result of the higher deposited MW power density.As discussed in subsection 4.1, by validating neutral species mole fractions and gas temperature for the first set of conditions, we consider the model to accurately reproduce these quantities also for the second set of conditions.

Model
A novel 1D radial fluid model for CO 2 MW plasmas has been developed.Due to interactions with other neutrals and charged species, CO 2 molecules dissociate into the plasma giving rise to a multi-component mixture of species.In this respect, the , and e), and an equation for the gas temperature (T g ).The fluid model is coupled to a MCF solver for the electron Boltzmann equation [36,37] that provides the reduced electric field (E/n) and electron properties of the discharge, such as electron impact rate coefficients, electron transport parameters, electron temperature (T e ).The coupling is discussed in the following subsections.One of the main assumptions considered in this work is that only variations of discharge parameters in the radial direction of the discharge vessel are taken into account.Moreover, the system is considered at steady-state, as in experiments.The importance of radial mass and energy transport for the description of CO 2 dissociation in the DIFFER reactor has been highlighted by Wolf et al [30], where a 2D fluid model has been used to model CO 2 reactive gas flows.Furthermore, note that axial transport is not completely neglected in the present model and the effect of axial gas flow is taken into account as a source/loss term for the mass and energy balance equations.The experimental parameters in tables 1 and 2 are used as input for the model.The calculated species mole fractions, gas and electron temperature are validated against spatially-resolved measurements performed at DIFFER.The model is also used to complement experiments, by providing calculations of charged species mole fractions that are not directly measured.

Mass balance equations
The mass conservation equation involves the total number density n and mass density ρ of the mixture that are defined as where m i , n i , and ρ i are the mass, number density and mass density of the ith species, respectively [40].From equations ( 1) and ( 2), the species mole fraction x i and the species mass fraction y i are defined as Moreover, in a multi-component reactive mixture, the mixture mass can be defined as In 1D cylindrical coordinates, the mass balance equation for the ith species in the mixture is expressed as with J i the diffusive mass flux of species i, ω i the production rate of species i due to chemical reactions, and Ω i accounts for the transport of species due to convective axial flow.The species production rate is calculated from with ν i j the stoichiometric coefficient, and R j the reaction rate for the jth chemical process.The effect of axial gas flow is accounted for as a source term in the right-hand side of equation ( 5) representing the inflow of CO 2 as where L is the plasma length that is estimated from measurements the emissivity profile of the 777 nm spectral line emission of atomic oxygen (O(3p 5 P) → O(3s 5 S)), and v axial is the axial speed of the inflow gas defined from the conservation of the mass flow rate as where S tube is the cross sectional area of the tube and ṁ = ρ env Φ is the input mass flow rate, with Φ the input gas flow rate and ρ env the environment mass density that is composed by CO 2 at 300 K. Similarly, the outflow for the ith species, other than CO 2 , at the end of the reactor is given by Note that equations ( 7)-( 9) are defined such that the continuity equation for the mixture as a whole is satisfied.Indeed, the same approach has been adopted by Viegas et al [23,38] in order to take into account gas flow effects in 0D and 1D model calculations.
The mass flux J i for the ith species [41,42] is dT g dr , (10) where q i is the charge of species i, μ i is the mobility of species i in the mixture, E amb is the ambipolar electric field in the radial direction due to radial charge separation, D i is an effective diffusion coefficient for species i in a multi-component mixture.Note that three contributions are included for the radial diffusion, that are (i) gradients of mole fractions, (ii) gradients in mixture mass, and (iii) thermophoretic forces, where the expressions for the thermal diffusion ratio k T i of the ith species are taken from the work by Fristom and Monchick [43].The mass diffusion of species i in a multi-component mixture is described by means of a diffusion coefficient in which D M i is a mixture-averaged diffusion coefficient for species i, and D T eff is an effective coefficient accounting for turbulent mixing [44].The expression for D M i is calculated according to [45] as where D i, j is the binary diffusion coefficient for the (i, j) pair of species, derived from first order collision integrals [46].In particular, collision integrals for neutral-neutral and electron-neutral interactions are taken from the work by Laricchiuta et al [47], while electron-electron and electron-ion collision integrals are tabulated by Mason et al [48].In case of an electron-heavy species interaction, the collision integrals are computed at the reduced temperature given by T * = (m e T g + m i T e )/(m e + m i ), with m e and m i the electron and heavy species mass, respectively, T e and T g the electron and the gas temperature, respectively [49].Note that equation (12), also used by Kubečka et al [50] and by Baeva et al [42], is a generalized Fick expression that is strictly valid under the trace-species approximation [51].A self-consistent approach to describe mass diffusion in a multicomponent mixture is based on the solution of momentum balance equations, as defined in [52,53].This approach has been considered in [24], where the Stefan-Maxwell equations have been solved numerically together with mass balance equations for a CO 2 -CO-O 2 -O-C mixture.Nevertheless, we note that, for the conditions under consideration in this work, deviations between neutral species mole fractions obtained using equation ( 12) and the ones obtained by a self-consistent model based on the solution of the Stefan-Maxwell equations do not exceed a few percent.The expression for the effective turbulent diffusion coefficient D T eff is taken from [30,54] as where ν T is the turbulent kinematic viscosity, and Sc T is the turbulent Schmidt number.In particular, supported by computational fluid dynamics simulations of the vortex flow including a heat source representing the plasma [30], it is assumed that ν T has a quadratic profile as a function of the radius: where R in is the tube inner radius, and the value of the coefficient ν T (0) is taken from the work by Wolf et al [30].Note that ν T (0) exhibits a dependence on both gas pressure and gas flow rate, as described in [30].A value of Sc T = 0.7 is considered in equation (13), as in works considering comparable gas flow rates [55].
The ambipolar diffusion theory, as reported by Oskam [56], for multi-ionic plasmas is used to describe the charged species radial transport.The same approach has been adopted by Hassouni et al [41] and by Wang et al [28].According to Oskam's theory, the ambipolar field is where the summation runs over all charged species including electrons, and the mobilities of different charged species μ i are calculated using the Einstein relation, from the gas temperature T g and the corresponding mixture-averaged diffusion coefficient at every radial position [54].Note that Oskam's model has been derived assuming that the total current density and the total charge density obtained by summing over all species are both equal to zero [56].These assumptions are typically violated when the plasma contains several ion species [42].However, as it will be shown in subsection 4.2, O + 2 ion is the dominant ion for a diffuse plasma, over the whole radius, and for a contracted plasma, further than 1 mm from the center of the tube.Hence, we expect that the current ambipolar diffusion model is still accurate, for the conditions under considerations.
The electron mass fraction is calculated from quasineutrality as y e = m e q e i =e where the summation runs over all charged species, except electrons [40].Note that, due to the usage of effective diffusion coefficients in equation (11), it is not guaranteed that diffusive mass fluxes sum up to zero.In other words, the continuity equation for the mixture as a whole does not hold.Hence, the mass fraction of CO 2 is obtained by enforcing the mass conservation as where the summation runs over all charged and neutral species, except CO 2 .Furthermore, the species mole fractions are derived from the respective mass fractions according to the relation In section 4, the species molar fractions calculated from the model are compared with the ones derived from experimental measurements.

Gas temperature equation
In the model, neutrals and ions are characterized by a translational temperature T g .Moreover, supported by detailed experimental observations [10,12], it is assumed that, at steady-state, the translational, rotational and vibrational temperatures of the molecules have all the same value T g , that is calculated from the following energy equation: where λ is the heat conductivity of the mixture, J i is the mass diffusive flux of species i (equation ( 10)), h i is the enthalpy of species i, P in (r) is the MW power density that is locally absorbed by the electrons, Q flow is the heat loss due to convection by the axial flow, and Q rad is the energy lost from the plasma by radiation.
The heat conductivity in equation ( 19) is calculated as the sum of a mixture-averaged conductivity λ M g and an effective turbulent conductivity λ T eff , as proposed by Gleizes et al [44]: A semi-empirical formula is used for calculations of λ M g : which is only dependent on the species molar fractions x i , and on the pure species heat conductivities λ i .According to Mathur et al [57], this expression gives results with errors in the range of a few percent.Moreover, in 1D model calculations, equation ( 21) is preferred to other semi-empirical formulas, such as the one proposed by Wilke [58], since it does not require computations of single species gas viscosities.The pure species quantities λ i are cast in a polynomial form, where the coefficients of the polynomials are tabulated in the NASA library [59].The effective turbulent component of the heat conductivity is written as where Pr T is the turbulent Prandtl number, and C p,mix is the mixture-averaged specific heat.A value of 0.85 is used for Pr T , as in simulations for a similar setup [30].Moreover, the expression of C p,mix is defined as where C p,i is the specific heat of species i [51].The pure species quantities, h i and C p,i , in equations ( 19) and ( 23) are taken from a NASA technical report [60].
The loss of heat due to convective cooling by the axial gas flow is accounted for in equation (19), by a term depending on the axial flow speed (equation ( 8)) and on the plasma length where the subscript 'env' denotes the environment species composition and temperature, that is considered as pure CO 2 at 300 K.A similar expression has been taken into account by Wang et al [28].The radiation transfer term in equation ( 19) is calculated as where N is the net emission coefficient of radiation, for a given temperature.The values of N are taken from the work by Aubrecht and Bartlova [61], where N is computed for a thermal CO 2 plasma, at various pressures and plasma thickness.
The local absorbed MW power density in equation ( 19) is estimated from measurements of the total plasma emissivity as The same expression for P in (r) has been adopted by Wolf et al [30] and the validity of equation ( 26) has been discussed in details by Viegas et al [38].We note that results of the model are highly dependent on the choice of the input peak MW power density (P in (0)).Hence, accurate estimations or measurements of P in (0) are sought after.In subsection 4.4, the effect of the uncertainties on experimental estimation is P in (0) on the model results are presented.The power density loss is also calculated from MCF as where the expressions for the elastic (Q e,el ) and inelastic (Q e,inel ) electron power losses are taken from [62] as where ν ei is the microscopic collision frequency for an electron of energy which results in a transfer of momentum to the ith species, ν e j is the microscopic collision frequency for inelastic collisions of electrons with the jth species, V j is the corresponding energy threshold, ν inv e j is the microscopic collision frequency of the superelastic inverse process of each individual excitation process, and the symbol • represents the average of the quantity within the brackets over the EEDF.Additionally, in the model, the electron temperature is defined as k B T e = 2 /3, where is the electron mean energy, although the EEDF deviates from a Maxwell-Boltzmann distribution.The local value of E/n is determined through an iterative procedure until the condition P in (r) ≈ P calc (r) is met, as described in subsection 3.4.Note that equation ( 27) is an approximate form of the electron energy balance equation that is valid when electrons are close to local field equilibrium and the electron energy transport terms are negligible [63].This assumption is valid for the conditions under consideration as the energy exchanges through inelastic collisions of electrons are the dominant contribution in the electron energy equation [23].
According to Dalton's law, the static pressure p g is defined as the sum of the partial pressures of the different species, thus leading to Equation ( 30) is used as local constraint to determine the total gas number density.

Boundary conditions
In the center of the tube, due to axial symmetry, a homogeneous Neumann boundary condition for the neutral species mass fractions and gas temperature has been considered.A flux boundary condition has been used for neutrals and ions at the tube wall (i.e. at r = R in ) [64,65]: where n is the normal vector pointing towards the wall, H(sgn(q i )μ i E amb ) is the Heaviside step-function that is set to one if the drift velocity is directed towards the wall and to zero otherwise, v th,i is the thermal velocity of the ith species, i.e.
and γ i is the fraction of particles that is adsorbed at the surface that is set to 10 −3 for neutral species, based on similar values for adsorption of O atoms in silica-like walls [66], and it is assumed to be equal to 1 for the ions [41].Moreover, recombination of atomic species into molecular ones the surface and secondary electron emission from the surface are not taken into account in the model.Note that the aforementioned assumptions on adsorption probabilities do not have an important influence on the model results for the conditions under consideration [67].In fact, the species concentration in the MW reactor is mainly determined by the gas phase chemistry and not by surface reactions.At the quartz wall boundary, the heat flux (J g ) is assumed to be continuous [40], such that where R in = 13.5 mm is the inner tube radius of the reactor, the thickness of the quartz tube is 1.5 mm, such that the outer radius is R out = 15 mm, λ w is the quartz tube heat conductivity, that is assumed to have a constant value of 1.4 W m −1 K −1 , and T g,env is the environment temperature for the outside part of the wall and it is set to 300 K.

Numerical solutions of the fluid equations and coupling with Monte Carlo flux
The finite volume method has been used to discretize the fluid equations ( 5) and ( 19) [68,69].As such, all the variables and sources have been defined at cell centers (i.e. the nodal points of the grid), whereas fluxes are defined at cell interfaces.The exponential scheme by Scharfetter and Gummel [70] has been employed for the discretization coefficients.The tridiagonal matrix algorithm has been employed for matrix inversion.
The fluid model is coupled to an MCF code, that provides calculations of the EEDFs, T e , electron rate and transport coefficients and electron power losses (equations ( 28) and ( 29)) in each cell node, for a given E/n profile and composition.The composition and E/n are calculated through an iterative procedure.The principles of the iterative solution procedure are similar to the ones adopted by Hassouni et al [67] and the schematic of the fluid-MCF model implementation is shown in figure 3. The simulation starts by setting different input parameters, such as the profile of MW deposited power (MWPD inp (r)), the pressure of the gas (p g ) and the gas flow rate (Φ), and an initial guess for profiles of reduced electric field (E/n), species mass fractions (y i ), gas temperature (T g ), and reduced angular frequency of the field (ω/n).A first set of MCF simulations is performed, in order to obtain estimations of the electron parameters averaged over the EEDF, under the aforementioned initial condition.In order to reduce the CPU time associated to MCF simulations, MCF runs are performed only for nodal points labelled by odd-indexes in the grid.Hence, a linear interpolation of MCF outputs from the odd-index to the even-index nodal points is performed.The accuracy of the calculations is expected to be affected no more than a few percent by this procedure, due to the high number of cells that are used in the grid leading to a cell size lower than 0.1 mm.The fluid model calculation starts by computing transport properties, thermodynamic parameters and source terms that are used for the numerical solution of the mass and energy balance equations (equations ( 5) and ( 19)).A first check is performed to ensure that the species composition and gas temperature reach numerical steady state.The numerical steady state is reached when the residuals of y i and T g are lower than a pre-defined tolerance ( ), that is set to 10 −6 .Typically, around 10 000 iterations are required for this condition to be fulfilled.Once the steady-state is reached, the MW power losses are computed from equation (27).The discharge maintaining E/n is determined by the condition that |P calc (r) − P in (r)| < , where is the tolerance for the power balance between the MW power deposited and the power losses by the plasma, that is set to 10 −2 .In this case, the E/n profile is adjusted iteratively in the code by multiplying the E/n found from the previous iteration by a term proportional to (P calc (r) − P in (r)), in order to fulfil the aforementioned condition.Typically, around ten iterations are required to satisfy the power balance relation.Note that, if E/n is updated, also MCF simulations are repeated, taking into account the updated composition calculated from the fluid model.
All the calculations presented in section 4 have been performed with a uniform grid of 150 cells, from the center of the tube to the inner tube radius.The total number of cells has been chosen in order to guarantee uniformity in the calculated quantities and source terms within the control volumes, as required by the finite volume approach [69].A relaxation factor of 0.99 is applied on discretization coefficients and chemical source terms of mass balance equations [71].As initial condition, the electron and O + 2 ion number densities are set to be equal to 10 2 cm −3 in the whole domain.Moreover, a background CO 2 gas at 300 K is considered as initial condition of the simulation.The total CPU time for a fluid-MCF simulation is up to 15 h on a 2.6 GHz Intel i7 single core processor, mostly due to the iterations required for the MW power deposited to be equal to the power lost.

Reaction scheme
The reactions scheme used in this work is mostly taken from the work by Viegas et al [23], that is an improvement of previous works by Kozák and Bogaerts [18] and by Koelman et al [19] for application to CO 2 MW plasmas at intermediate pressures.In this subsection, the main features of that scheme are summarized, together with the changes in the chemistry set that have been adopted for this work.

Thermal chemistry.
The reactions scheme for neutralneutral collisions has been taken from the GRI-MECH 3.0 database [72].In particular, table 3 presents forward reaction rate coefficients (k f,i ) for CO 2 , CO, O 2 , O and C.Moreover, in van de Steeg et al [10] it has been shown that the use of this dataset of reactions gives the best agreement between calculations and experiments of total (volume averaged) CO production and input MW power, where k 0,i is the low-pressure limit rate coefficient, whose expression is taken from [72].The rate coefficients for the reverse reactions (k r,i ) have been obtained considering that, at equilibrium, the rate-of-progress of all reactions is zero, which yields the following relationship between the forward and backward rate coefficients where K eq,i is the equilibrium constant for the ith reaction [73].
The equilibrium constant is related to the change in Gibbs free energy across a reaction ΔG i by where Δν i is the difference between the stoichiometric coefficients of products and reactants.Additionally, according to [23], it is assumed that, at steady state and for the high T g calculated in the plasma core, electronically excited species have a negligible contribution, and all vibrational/rotational states of the ground electronic state of the molecules are in thermal equilibrium, as measured by van de Steeg et al [10].
3.5.2.Electron impact reactions.Electron impact reactions used in the fluid model can be found in table 4, where all neutral and charged collision partner species have been considered in their ground electronic state.For a collision of an electron with a neutral atom/molecule, the reaction rate coefficient is given by where is the electron kinetic energy, k j and σ j ( ) are the electron impact rate coefficient for the jth collisional process and the corresponding cross section, respectively, and f 0 ( ) is the EEDF that is calculated from MCF. Normalization of the EEDF is such that Due to a lack of cross sections for electron-ion interactions (E20-25), rate coefficients for those reactions have been computed from expressions that are uniquely dependent on the electron temperature (T e ).
In addition to the reactions reported in table 4, the following processes have been included for accurate EEDFs calculations in MCF.Note that those processes are not included in the fluid model reactions scheme, but only in MCF.In particular, electron impact cross sections with CO 2 molecules have been taken from the Biagi database of LXCat [74].Specifically, the elastic momentum transfer for the ground electronic state of CO 2 has been modified by taking into account a Boltzmann population at T g of the excited vibrational levels of the ground electronic state and rotational levels of the ground vibrational state, as described in [37].Additionally, cross sections for Table 3. List of neutral thermal chemistry reactions.Rate coefficients are provided for the forward reactions.Rate coefficients for the reverse reactions are calculated assuming thermodynamic equilibrium.Units for rate coefficients for two-body reactions are cm 3 s −1 and for three-body reactions are cm 6 s −1 , n is the total number density in cm −3 , P r is the reduced pressure as defined in [72], and T g is the gas temperature in K. electron impact dissociation of CO 2 have been taken from the work by Polak and Slovetsky [81], as the calculated rate coefficients from these cross sections are comparable with experimental measurements [37,84].According to the recommendation by Morillo-Candas et al [84], Polak and Slovetsky's cross sections have not been included in MCF simulations, but they have been used only for the calculations of the corresponding rate coefficients.The use of Polak and Slovetsky's cross sections within the Biagi's set of electron impact cross sections has been found to provide good agreement between calculated and measured total dissociation rate coefficients, within a few percent, for E/n up to 110 Td [37].The total ionization cross section from the Biagi database has been replaced by the ones taken from the Itikawa database of LXCat [75] that reproduce the same total ionisation cross section, but include also the dissociative contribution to the formation of different ions.Fridman's approximation [9] has been used to obtain electron impact cross sections for stepwise excitations from excited vibrational levels, from ν = 1 to ν = 21 for the asymmetric stretching vibrational mode of the ground electronic state.The same approximation has been adopted by Pietanza et al [22] and by Berthelot and Bogaerts [85].Superelastic collisions cross sections from vibrationally excited states of the ground electronic state of CO 2 have been calculated using the formula of Klein-Rosseland [86].

#
Electron impact cross sections with CO molecules have been taken from Biagi's code Magboltz v11. 10 [76].This set replaces the IST-Lisbon dataset [87] that has been used in the work by Viegas et al [23].Anisotropic scattering for rotational collisions has been included according to the dipole-Born anisotropic scattering model described in [88].As for the CO 2 cross sections set, the total ionization cross section from Magboltz v11.10 has been replaced with different dissociative ionization cross sections that are included in the Itikawa database of LXCat [75].Cross sections from the work by Laporta et al [89] for electron impact vibrational excitations from vibrational levels of the ground electronic state have been included, for levels from ν = 1 to ν = 80.These cross sections have not been included in our previous work [23], but they are present in the work by Pietanza et al [22].The population of rotational and vibrational levels has been considered as a Boltzmann distribution at T g .Moreover, superelastic collisions from those levels have been included using the micro-reversibility principle [86].Additionally, note that electron impact dissociation of CO molecules in the ground electronic state is not included in the Biagi's set.Therefore, the cross section for this process has been taken from the work by Cosby [90].The same cross section has been included in the recent IST-Lisbon set [87].
As in [23], electron impact cross sections with O atoms and C atoms are taken from the IST-Lisbon [77] and the B-spline atomic R-matrix databases [78], respectively.
Electron impact cross sections with O 2 molecules are taken from the Biagi database of LXCat [74].These cross sections replace the ones from the IST-Lisbon database [77] considered in [23].Indeed, Biagi's cross sections have been preferred to the IST-Lisbon ones due to (i) consistency with the MCF method used in this work, and (ii) the use of an elastic momentum transfer cross section, instead of the effective cross section that is proposed in [77].Ionization cross sections have been taken from the Itikawa dataset [75], in replacement of the total ionization cross section that is included in the Biagi's set.A Boltzmann distribution for the first 11 vibrational levels (i.e.ν = 0-10) of the ground electronic state of O 2 has been considered.Moreover, stepwise vibrational excitations/de-excitations between vibrational states characterized by quantum numbers from ν = 0 to ν = 10 have been included.Cross sections for these processes have been taken from the work by Laporta et al [91].Electron impact de-excitation cross sections have been calculated using the Klein-Rosseland formula [86].
In this work it is also assumed that the plasma is weakly ionized, such that self-collisions of electrons (i.e.electronelectron collisions) are negligible.This assumption is discussed in more detail in section 4.

Ion-ion and ion-neutral reactions. Ionic conversion (I1-20), electron detachment (I21-24), ion-ion collisions
Table 4. List of charged particle reactions: electron-impact ionisation, attachment and dissociation; electron-ion recombination.Units for rate coefficients for two-body reactions are cm 3 s −1 and for three-body reactions are cm 6 s −1 , T e is the electron temperature, in K, for reaction E5 and in eV for the remaining reactions, M represents any neutral species, and T g is the gas temperature in K.For other rate coefficients labelled by EEDF, the value is calculated from the EEDFs computed with the MCF.Note that these reactions are used in the fluid model.See text for a description of electron impact reactions included only in the MCF solver.(I25-28) are reported in table 5. Most of the rate coefficients for these reactions have been taken from the work by Beuthe and Chang [17] and by Viegas et al [23].

Associative ionization.
Associative ionization reaction rate coefficients are reported in table 6 (reactions A1-3).In particular, as described in section 4.2, associative ionization reactions are particularly important at high pressures (i.e.above 100 mbar), since they are the main sources of electron production.The expression for the rate coefficients has been taken from the work by Park et al [94] as where the values for the coefficients C and α and for the activation energy E a (in K) are given in [96], and T g [K] is the gas temperature in K.Note that these rate coefficients have positive pre-exponential powers (i.e.α 1).However, as stated by Park [96], it is unlikely that α would remain positive for T g above 6000 K. Therefore, in this work, α = 0 has been assumed.The resulting expression for k ai is then equivalent to an extrapolation of these reaction rate coefficients for gas temperatures above 6000 K, as recommended by Park [96].

Neutral species composition and gas temperature
In this subsection, calculated neutral species mole fractions and gas temperature along the tube radius, corresponding to the center of the discharge vessel in the axial direction, are shown.Results of the calculations have been compared with measurements obtained by means of Raman scattering [13,24].The input parameters for the model are given in table 1.
The measured and calculated gas temperatures for the 100 mbar and 150 mbar cases are shown in figures 4(a) and (b), respectively.These conditions have been chosen as examples of a diffuse and a contracted discharge.Specifically, the calculated T g is up to 4900 K and 7200 K for the 100 mbar and 150 mbar cases, respectively.As a result of the higher peak in the MW power density (P in (0)), the contracted case is characterized by (i) higher gas temperature in the center of the tube, and (ii) larger gradients in the radial profiles (figure 4(b)), than in the diffuse case.As shown in subsection 4.2, the nonuniformity of T g affects also the spatial distribution of charged species concentrations.Overall, a good agreement has been found between measurements and calculations of T g , within   2 + e 1.86 × 10 −11 exp(−80 600/T g ) [94] the first 10 mm.In this respect, note that the expression for turbulent viscosity that has been used in this work (equation ( 14)) is taken from the work by Wolf et al [30], in order to ensure that radial turbulent transport is representative of the experimental conditions.Measured and calculated mole fractions for the main neutral species (CO 2 , CO, O, O 2 , and C) for the 100 mbar and 150 mbar cases are shown in figures 5(a) and (b), respectively.The agreement is good.In particular, figure 5 shows that the transition from a diffuse (100 mbar) to a contracted (150 mbar) discharge is accompanied by a change in neutral species composition, from a partially dissociated to a fully dissociated core, near the center of the tube.This transition is a direct consequence of the increase in the peak T g showed in figure 4. The highest discrepancies between simulation results and experiments have been found for the O atom mole fractions at 150 mbar (figure 5(b)), for radial positions within 5 mm.These discrepancies can be due to (i) an underestimation of the O atom formation rates from thermal dissociation of CO and O 2 at high T g (reactions N8-11 from table 3), or (ii) the local conservation of 1:2 stoichiometric ratio of C and O atoms, that is assumed in experiments for the quantification of the O atom Raman cross section [24].Despite the simplifications in the fluid model for the chemistry and the gas flow, figures 4 and 5 show that calculations are overall in good agreement with Raman measurements.In particular, although the calculated peak T g reaches a maximum of about 7200 K at 150 mbar, C atom mole fractions are always negligible.This result is also in agreement with calculations by Sun et al [29], where C atom mole fractions do not exceed values of 0.02 in similar conditions.Nevertheless, as discussed in subsection 4.3, the presence of C and O atoms is fundamental for electron production through associative ionization, in the core of a CO 2 contracted plasma.Another assumption taken in the 1D model is that electronically excited states are negligible.Although these species can be important in pulsed MW conditions [97], there is no experimental evidence for a significant presence of those excited states in CW conditions that could potentially invalidate present measurements or calculations [24].
The main reaction rates for thermal dissociation of CO 2 (i.e.reverse of reactions N3-7 from table 3) calculated from the model are shown in figures 6(a) and (b) for the 100 mbar and 150 mbar cases, respectively.In both conditions, the CO 2 dissociation is mostly due to collisions with neutral species, through the following reactions (N3-5): where M = CO 2 , CO, and O 2 .The produced atomic oxygen atoms lead to further dissociation of CO 2 through the reaction (N6): These reactions are enhanced by the high gas temperature [10,29].Note that, even for the case at 100 mbar, the rate of electron-impact dissociation of CO 2 is at least two orders of magnitude lower than the ones for thermal dissociation.Indeed, electron-impact dissociation of CO 2 is found to be dominant only at lower pressures, around 60 mbar [10,23].For radial positions larger than 5 mm, towards the edge of the tube, the CO 2 content increases due to a drop in the gas temperature.Reaction rates for CO 2 dissociation are largely reduced at lower T g .
To summarize, for the high gas temperatures considered in this work, CO 2 dissociation is mainly driven by thermal reactions and not by electron impact reactions.Hence, in agreement with [10], experimental results of conversion and energy efficiency can be explained by considering only thermal composition and turbulent gas flow dynamics, without requiring additional vibrational non-equilibrium effects [30].

Charged species composition and electron temperature
In the previous subsection, it has been shown that CO 2 dissociation is mainly driven by collisions between neutral atoms and molecules.However, charged particle kinetics is crucial for understanding contraction phenomena in MW plasmas [33].An increase in T g and ionization degree in the plasma core in the contracted regime has been observed with respect to diffuse conditions [14].Since the input power is coupled to the electrons and transferred by collisions to the heavy species, plasma contraction is strongly related to a change in the spatial distribution of electron density.
Figure 7 shows the calculated electron and ions number densities for the diffuse (100 mbar) and contracted (250 mbar) cases.The input parameters for the model are given in table 2. At 100 mbar and 250 mbar and 1000 W, the calculated neutral species mole fractions are similar to the ones shown in figure 5 (100 mbar and 150 mbar with 700 W), with a transition from a partially dissociated to a fully dissociated plasma core.In figure 7, results of the calculations have been compared with measurements of electron number densities obtained with Thomson scattering.In figure 7(a), it is shown that O + 2 is the dominant ion.This result is in agreement with the one obtained by Wang et al [28] and by Viegas et al [23], where O + 2 has been found to be the dominant ion in the discharge.At  , respectively.Note that ionic conversion tendentiously leads to the formation of ions with lower ionization potential than the reactants [23].Excellent agreement can be observed between measured and calculated electron number densities for the 100 mbar case for the first 3 mm along the tube radius, within the experimental errors.The transition from diffuse to contracted regime is characterized by an increase in the peak electron number density, from around 3.5 × 10 18 m −3 at 100 mbar (figure 7   at 250 mbar point out that the plasma in contracted regime is not composed by a single ionic species, but multiple ions are formed, namely O + 2 , O + , and CO + , with comparable number densities.Specifically, two spatial regions along the tube radius can be identified in the contracted regime: (i) the inner region, within 1 mm from the center, in which O + and O + 2 ions are dominant and an important increase of C + and CO + number densities can be observed, compared to the 100 mbar case, and (ii) the outer region, from 1 mm up to the inner wall of the tube, where the dominant ion is again O + 2 and, as a consequence of quasi-neutrality, this ion sets the electron number density calculated from the model at the plasma edge.In particular, results of the model suggest that CO + is mainly produced by associative ionization (A2): In other words, collisions of O atoms with charged and neutral species determine the electron number density profile in the contracted regime.The model under use takes a profile of power density as input, with a Gaussian shape with FWHM defined from the measured atomic oxygen emission intensity (O(3p 5 P) → O(3s 5 S)).Indeed, it is assumed that atomic oxygen emission intensity, electron density and power density follow the same radial profile.As such, it is relevant to compare the measured and simulated radial profiles of electron density with the assumption based on atomic oxygen emission intensity.The normalized radial profiles of intensity measured from atomic oxygen emission and of calculated electron number density are shown in figures 8(a) and (b) from the 100 mbar and 250 mbar cases, respectively.
At 100 mbar, the normalized profile of electron number density is broader than the emission intensity profile, whereas, at 250 mbar, the normalized electron number density and emission intensity profile overlap.This can be explained by the fact that optical contraction occurs at 100 mbar, when the plasma is not fully contracted.Moreover, at 100 mbar, the radial gradient  of O-atoms mole fraction is comparable to the radial gradient of the electron number density profile (figures 5(a) and 7(a)), whereas, at 250 mbar, the O-atoms mole fraction is almost radially homogeneous in the region of high n e .This result is in agreement with the ones obtained by Viegas et al [38], where the optical contraction phenomenon has been predicted for gas pressures below 150 mbar and this phenomenon has been attributed to the inhomogeneous O-atoms mole fractions radial profile.
The change of ionic composition in the plasma core with pressure is shown in figure 9, where an overview of the peak electron number density measured and calculated at different pressures in the plasma core, from 100 mbar to 400 mbar, is presented.The peak number densities for the dominant ions calculated from the model have been plotted in the graph, as well.Specifically, the model reproduces the increase in electron number density in the plasma core with pressure, that is associated with contraction of CO 2 MW plasmas at 1000 W. Correspondingly, the measured ionization degree (n e /n) in the plasma core increases from 2.5 × 10 −5 at 100 mbar to 2 × 10 −4 at 400 mbar.Discrepancies exceeding experimental errors can be observed for the 400 mbar case, probably due to an underestimation of the absorbed MW power density used as input to the model.The sensitivity of the calculated n e with the input MWPD inp (0) is quantified in subsection 4.4.Note that, as pressure increases from 100 mbar to 400 mbar, a change in the ionic composition of the plasma core has been observed in the model, with significantly increasing number densities of O + and CO + .The number density of the negative ion O − is always below 10 17 m −3 , thus this ion is negligible compared to positive ions, mainly due to detachment processes (reactions I21-24 in table 5).Moreover, for pressures above 200 mbar, O + is the dominant ion in the plasma core.Thus, from this study, it appears that the contraction phenomenon in CO 2 MW plasmas is related to the formation of atomic ionic species in the plasma core.This result is in agreement with the study by Viegas et al [23], where the importance of ion composition for CO 2 MW plasma contraction has been studied with a 0D model.
Diffuse and contracted plasmas are also characterized by different electron temperatures.In particular, radial profiles of calculated and measured T e at 100 mbar and 250 mbar are shown in figures 10(a) and (b), respectively.Measurements are limited within a 3 mm range from the center of the tube, as the Thomson signal is reduced by a lower electron number density and by an increase in CO 2 number density at outer positions along the tube radius.Calculated electron temperatures (T e ) are derived from MCF as T e = 2 /3, where is the electron mean energy.
On the one hand, the diffuse case (figure 10(a)) shows that T e is constant within the first 3 mm along the tube radius.For this case, measurements and 1D model results are in good agreement, within the experimental errors.On the other hand, T e is found to increase outwards for the contracted case (figure 10(b)), from a value of 1 eV in the center of the tube to around 2 eV at 3 mm.This increase is not reproduced by the model, where a concave T e profile has been obtained.Note that a similar increase in T e near the wall has been measured also by Carbone et al [98] with Thomson scattering in contracted MW plasmas operated with argon at intermediate pressures, from 20 mbar to 88 mbar.This trend has been explained by Jimenez-Diaz et al [99] as a consequence of the non-homogeneity of the ionization frequency across the radius, using a 2D fluid model.However, note that the model by Jimenez-Diaz et al [99] considers Maxwell-Boltzmann EEDFs in the calculations.This assumption has been demonstrated to not be accurate for argon MW plasmas by Ridenti et al [34], where the increase in T e has been explained as a result of an implicit Maxwellian assumption for the EEDF in Thomson scattering.In [34], a Thomson temperature (T e,Th ) is defined as the temperature obtained by fitting the non-equilibrium EEDF with a Maxwell-Boltzmann distribution in the 0-2 eV energy range.As the low-energy part of the distribution (i.e. the head of the EEDF) gives most of the contribution to the measured Thomson signal, the T e,Th defined in this way is comparable with the electron temperature measured by Thomson scattering.Following the idea by Ridenti et al [34], the T e,Th values derived from the corresponding non-equilibrium EEDFs are shown in figure 10, for the cases at 100 mbar and 250 mbar.In particular, for the 100 mbar case (figure 10(a)), differences between calculated T e and T e,Th are within 15%.At 250 mbar 10(b)), the T e,Th is found to be in good agreement with experimental values, as the outward increase of measured T e has been reproduced.Therefore, the measurement of convex T at 250 mbar (figure 10(b)) are attributed to the Maxwellian assumption implicit in the Thomson scattering diagnostic and not to a higher mean electron energy.Note that results presented in this section have been obtained without the inclusion of electron-electron collisions.Indeed, with the addition of self-collisions of electrons in the MCF code, the EEDF is expected to approach the Maxwell-Boltzmann equilibrium distribution [100].The importance of this effect would be larger in the center of the tube than at the edge, since the ionization degree (n e /n) is highest at r = 0. Nevertheless, EEDF calculations using the two-term solver BOLSIG+ [101,102] confirmed that the addition of self-collisions of electrons do not have an appreciable effect on the previous results.In fact, rate coefficients for electron impact ionization of CO and O obtained including electron-electron collisions in BOLSIG+ are within 5% with respect to the ones obtained by the same solver by neglecting this contribution.Therefore, from the present considerations, it is expected that these collisions would not change the results in this section by more than a few percent, although this mechanism would contribute to increase ionization in the core, and thus would be an extra contribution to contraction.
Central E/n values calculated at different pressures are shown in figure 11, where E = E 0 / √ 2 is the mean square value of the MW field, with E 0 the maximum amplitude of the field in the course of a period oscillation [62].Figure 11 shows that E/n in the plasma core is decreasing with increasing pressure.The same result has been reported by Viegas et al [23], Pietanza et al [22], and Groen et al [103].In particular, Pietanza et al [22] obtained lower E/n values than the ones calculated in this work, in the range 48-50 Td and 25-42 Td for the 100 mbar and 250 mbar cases, respectively.Nevertheless, their calculated peak electron number densities are close to the ones reported in this work, within a few percent.The main differences between the present study and the one by Pietanza et al [22] are related to input parameters, different dimensionality of the models, gas temperature and mass transport that are considered in the models.In fact, in [22], a 0D model has been used in which the MW power density has been assumed to be uniform in the discharge, whereas a radial profile has been considered in this work (equation ( 26)).Moreover, the peak power densities values in [22] are lower than the ones used here.In Viegas et al [23], it has been found that E/n decreases from around 80 Td at 100 mbar to 67 Td at 300 mbar.However, in the same work, it has also been shown that this parameter is subjected to large uncertainties depending on mass transport and input gas temperature.Moreover, as described in subsection 3.5, different electron impact cross sections sets have been used by Pietanza et al [22] and by Viegas et al [23], with respect to this work.The central E/n values obtained in this work are in good agreement with estimations by Groen et al [103].However, note that different methods have been used for calculations of E/n.On the one hand, in [103], a 3D electromagnetic model has been used, where the field has been obtained by varying the MW power through impedance matching.On the other hand, in this work, the field has been considered as an internal parameter of the simulation  and has been obtained by an iterative procedure under the condition that the plasma power losses adjust to compensate for the total MW deposited power, as described in subsection 3.4.
Measured and calculated T e and T g at r = 0 are shown in figure 12, for different pressures from 100 mbar to 400 mbar.Opposite trends are found for T e and T g with increasing pressure.In particular, as a result of the decrease in the central E/n (figure 11) with pressure, T e decreases from 1.5 eV at 100 mbar to 0.6 eV at 400 mbar, whereas T g increases slightly from 0.5 eV to 0.6 eV.Note that this equilibration is not predicted by the 0D model results by Viegas et al [23] and Pietanza et al [22], where T e has been found to increase with increasing pressure.Sources for these discrepancies are discussed in subsection 4.4.
To summarize, the present comparison of measurements and calculations of discharge parameters in CO 2 MW discharges showed that (i) the diffuse regime is characterized by a single dominant ion species (O + 2 ) and T e > T g in the plasma core, (ii) the contracted regime is characterized by a multiionic composition, with higher values of O + ion number density than the contracted case and T e ≈ T g in the plasma core for p g > 200 mbar.

Insight into plasma contraction from model results
This section investigates the charged particle kinetics and highlights the main chemical reactions leading to the results in the previous section.In particular, in figure 9, it has been shown that the peak electron number density increases from 2.5 × 10 18 m −3 at 100 mbar to around 3.0 × 10 19 m −3 at 250 mbar.As described by Zhong et al [104,105], it is expected that higher n e would require higher T e values in the center of the tube, in order to sustain the discharge in the contracted regime.However, in figure 12, it is shown that the electron temperature at r = 0 decreases from 1.5 eV to 1.0 eV, as the pressure increases from 100 mbar to 250 mbar.In order to investigate the mechanisms leading to this decrease in T e , the main electron production rates are plotted in figure 13.At 100 mbar (figure 13(a)), the discharge is mostly sustained by electron detachment from CO and O − (reaction I21 in table 5) and by electron impact ionization of CO and O in the center of the tube (reactions E6 and E11 in table 4).At 250 mbar (figure 13(b)), the rates of associative ionization for electron production are several orders of magnitude larger than those  for electron-impact ionization and electron detachment.In particular, the dominant ionization reaction is the one producing CO + (reaction A2 of table 6): Moreover, note that models neglecting associative ionization reactions cannot predict accurate values of T e in both diffuse and contracted regimes.For example, the work by Pietanza et al [22] reports T e values in the range 0.86-1.20 eV at 100 mbar and 1.20-1.46eV at 250 mbar, that are not compatible with the present experimental measurements.Indeed, without associative ionization mechanisms, the only source for electron production included in that model is electron impact ionization, whose rate increases at higher T e values.Thus, it is possible to conclude that, as a consequence of associative ionization, higher electron production rates can be achieved in contracted conditions, without the need of higher T e values, with respect to the diffuse case.
The main loss rates for electrons at 100 mbar and 250 mbar are shown in figures 14(a) and (b), respectively.Figure 14 shows that diffuse and contracted regimes are characterized by a different electron kinetics.In particular, at 100 mbar (figure 14(a)), electron attachment to O 2 (E15) is the dominant electron loss mechanism in the center of the tube (within 4 mm radius), whereas electron attachment to CO 2 (E13) is the dominant electron loss mechanisms towards the edge.At 250 mbar (figure 14(b)), electron loss rates are dominated by dissociative recombination of electrons with CO + (E22) and O + 2 (E24) ions.Moreover, as pressure increases, attachment is suppressed.This result is in agreement with the work by Wolf et al [15] and by Fridman and Kennedy [106] stating that, in the contracted regime, CO 2 discharges are controlled by electron-ion recombination.Moreover, on the one hand, at 100 mbar, the radial electron diffusion rate is of the same order of magnitude as the electron-ion recombination rate.On the other hand, at 250 mbar, the radial electron diffusion rate in the center of the tube is a factor 5 lower than the rate of dissociative recombination of electrons with CO + .Hence, radial diffusion of electrons is suppressed at high pressures.
To summarize, the following considerations about contraction in CO 2 MW plasmas can be deduced.At higher pressures, the gas temperature increases in the plasma core.Hence, associative ionization of C and O atoms becomes the dominant mechanism for electron production.This mechanism leads to increasing peak n e , without the need of increasing T e at r = 0, hence generating a self-reinforcing cycle between heating and increasing n e in the plasma core.As a result of a higher T g at r = 0 with increasing pressure, a change of composition occurs in the plasma core.Hence, the plasma becomes fully dissociated into CO and O in the center and molecular ions dissociate into atomic ones.In particular, for pressures above 200 mbar at 1000 W, O + 2 is no longer the most abundant ion in the center, in favour of O + .The transition in dominant ion has also been noted by Martinez et al [33] and by Ridenti et al [34], for MW plasmas operated in argon at atmospheric pressure.Moreover, the importance of the change in neutral and charged species composition for discharge contraction has also been highlighted by Viegas et al [23].However, in [23], it is hypothesised that the change of electron impact ionization rate coefficients with composition is a key mechanism for the increase of n e in the plasma core.This hypothesis is excluded by the results of this work, as associative ionization rates are dominant at high pressure, with respect to electron impact ionization rates.From this study, it is possible to conclude that non-uniform gas heating along the tube radius is the responsible mechanism for discharge contraction in CO 2 .The same mechanism has been hypothesised by Viegas et al [23] for CO 2 MW plasmas.However, in [23], radial gradients in gas temperature and species mole fractions have not been taken into account, as the focus of that work was the simulations of the main processes in the plasma core using a 0D plasma chemistry model.Indeed, due to gradients in the radial profile of T g , molecular ions are efficiently dissociated in the center of CO 2 MW plasmas, therefore reducing the number of molecular ions for dissociative recombination [33].Additionally, as pressure increases, radial transport of electrons and attachment of electrons to O 2 molecules are suppressed.Therefore, in contracted regime, charged species are not easily transported towards the edge of the tube, as in the diffuse case, and the radial extension of the plasma is uniquely determined by a balance between electron-ion recombination and associative ionization.19)) has been estimated from experimental measurements of the 777 nm spectral line emission profile of atomic oxygen [14].The estimation of this parameter from the emission intensity profile is based on several assumptions [38].Specifically, (i) a Gaussian profile of emission intensity, electron number density and MW power density has been assumed, (ii) the gas composition and reduced electric field have been considered to be radially homogeneous, and (iii) proportionality between emission intensity and electron density has been assumed.Since these assumptions are not fully met for the conditions under investigation (see figure 8(a)), large uncertainties are associated to P in (r), up to 50% [23,38].In order to investigate the effect of these uncertainties on the model results, calculations have been performed increasing P in (0) by 50%, as well as decreasing this parameter by the same percentage.Since volume integrated MW power density is conserved, the axial and radial FWHM are multiplied by a factor 1.26 for a 50% decrease in P in (0) or by a factor 0.79 for a 50% increase in P in (0).Hence, the total (volume integrated) input MW power is unchanged.Results of the calculations of n e and T e at r = 0 as a function of pressure are shown in figures 15(a) and (b), respectively.
For a 50% increase in P in (0) with respect to the original value, n e increases by a factor between 3 and 10.In fact, higher values of P in (0) lead to higher peak T g through equation (19) and, in turn, to higher rates of associative ionization (i.e.reactions A1-3 in table 6), that are the main mechanisms for electron production.Conversely, a 50% decrease in P in (0) results in lower ionization rates from reactions A1-3, which in turn leads to lower electron number densities than the experimental ones.Nevertheless, the difference between n e obtained from the original model and the ones obtained by lowering P in (0) is not a factor 10 as before, but only a factor 3. Figure 15(b) shows that lower T e values are obtained by increasing P in (0) by 50%.Indeed, for higher P in (0) associative ionization becomes increasingly important.Hence, a lower T e value is sufficient to maintain the discharge and match P in (r) in the model.Therefore, it is possible to conclude that a change in the input peak value of the MW power density drastically affects the calculations of T g and, in turn, of n e .Moreover, simulation results point out that the original input P in (r) is the most appropriate for reproducing experimental values of n e , for all the conditions considered in this work.
Changing the value of P in (0) directly affects E/n calculations and the T e values that are computed from MCF.In order to avoid that uncertainties in P in (0) affect model results, a possible alternative solution is to take as input a radial n e profile, instead of P in (r).This procedure has been adopted by Viegas et al [38], where a 1D fluid model has been used to study optical contraction in CO 2 MW plasmas.Nevertheless, note that measurements of n e are limited to radial positions within 3 mm from the center of the tube and accurate radial profiles of electron density are not directly accessible from experiments.Therefore, in [38], a collisional-radiative model for O-atoms has been coupled to the fluid model, in order to calculate the electron properties of the discharge through an iterative procedure.Another possibility would be to include the solution of Maxwell's equations for the electromagnetic field and a momentum balance equation for the electrons.As such, the model would become fully self-consistent and the profile of the MW absorbed power density could be directly computed from the Ohmic heating equation [63].However, note that the forward-vortex reactor used at DIFFER breaks the axial symmetry due to the rectangular MW cavity.Hence, it would be better to use 2D or 3D self-consistent plasma models in this case, in order to better reproduce the experimental set up.

Effect of associative ionization.
As mentioned in the previous subsection, associative ionization from the interactions between C and O atoms is the dominant mechanism for electron production in the contracted CO 2 plasma.However, large uncertainties are associated with the rate coefficient for this process.The following three different expressions for associative ionization rate coefficients have been found in literature: • In this work, the rate coefficient for the reaction C + O → CO + + e has been assumed to be equal to where T g [K] is the gas temperature in kelvin.The pre-exponential factor and the activation energy in equation ( 39) above have been taken from the work by Park et al for the same reaction [94].However, as recommended in [96], the pre-exponential power has been assumed equal to zero (see subsection 3.5).Note that the expression of k C-O ai assumed by Park et al [94] is derived from the N + O associative ionization reaction that has similar activation energy and pre-exponential factor.In fact, associative ionization involves the interaction of excited C and O atoms [107] and this effective rate coefficient considers C and O atoms as a whole, with excited states populated according to a Boltzmann distribution.This formulation is consistent with the assumptions made in the model., where a Boltzmann population is assumed for the electronically excited states.The rate coefficient for this process has been assumed equal to Since the C + O associative ionization reaction and N + O reaction have comparable activation energies [94], it is common in chemistry models to assume that the rate coefficient for the C + O reaction is equal to the one for the N + O reaction.We should notice that, in the work by Popov [107], the rate constant for associative ionization of N + O involves exclusively electronically excited Natoms.This is why this rate constant is so high and has a very weak dependence on the gas temperature.Hence, it should be applied only in models solving the mass balance equation for separate electronically excited states, which is not the case in this work.
Figure 16 shows a plot of these rate coefficients as a function of the gas temperature, in the range 300-8000 K. k O-O ai is lower than k C-O ai by about two orders of magnitude at 7000 K, due to the higher activation energy associated to the O + O reaction than to the C + O one.Furthermore, k N-O ai is almost constant in the gas temperature range considered and it is more than two orders of magnitude higher than k C-O ai , at 7000 K. Since C + O associative ionization is the main mechanism for electron production, it is expected that the use of a very different rate coefficient for this process would lead to high deviations in n e and T e .Figures 17(a) and 17(b) show results of the model for n e and T e at r = 0, when different rate coefficients from equations ( 40) and ( 41) are used instead of equation (39), that is the one adopted in the original model.As expected, the use of k O-O ai instead of k C-O ai leads to lower n e values than the ones calculated from the original model, due to the fact that k O-O ai is much lower than k C-O ai .At the same time, higher T e values than the measured ones have been obtained, by about 1.0 eV.Similar values have been obtained by Viegas et al [23], where electron temperatures in the range 1.8-2.0eV, for pressures between 100 mbar and 300 mbar, have been obtained.Indeed, in [23], k O-O ai has been used also for the C + O associative ionization reaction and this assumption leads to large deviations from experimental measurements.Moreover, note that the 0D model results in [23] predict an increase in T e with pressure, that is in contrast with the Thomson scattering measurements reported in this work.These discrepancies are due to the fact that, in [23], both n e and T g from experimental measurements have been used as input for the model.In particular, the peak T g values used in the 0D model are lower than the ones calculated here, by about 800 K. Hence, lower associative ionization rates than the ones in this work have been obtained in [23].
Assuming that the associative ionization rate coefficient for the C + O reaction is the same as the one for the N + O reaction of Popov [107], higher peak n e values have been obtained (figure 17(a)), with respect to the ones calculated with the original model.In fact, the rate coefficient for associative ionization reported by Popov [107] (equation (41)) is exceedingly higher than the one used in this work (equation (39)), thus higher electron production rates have been obtained.Figure 17(b) shows also that T e is largely affected by using k N-O ai , instead of k C−−O ai .In particular, T e is almost constant with pressure when using the rate coefficient by Popov [107] in the model.This result is expected, since k N-O ai is almost constant in the range 300-8000 K (figure 16), hence this associative ionization mechanism would be the dominant electron production rate, even at lower pressures that are typical of diffuse plasmas.
To summarize, calculations of n e and T e in CO 2 MW discharges are largely affected by the choice of associative ionization rate coefficients.Indeed, by changing the rate coefficient for associative ionization from C and O atoms, discrepancies between calculations and measurements exceeding experimental errors have been found.Therefore, the present results point out that future plasma chemistry models for plasmas featuring gas temperatures of the order of thousands K should take into account accurate expressions for associative ionization rate constants.These results show that the expressions for rate constants of associative ionization reported in table 6 are the ones giving the best agreement between calculations and experiments.However, further validation of associative ionization rate coefficients against shock-tube measurements at high T g is needed.

Conclusions
A vortex-stabilized CO 2 MW plasma at a fixed MW input power of 700 W and 1000 W and different pressures, from 100 mbar to 400 mbar, has been studied through comparisons between simulations and experiments.For this purpose, a 1D radial fluid model that takes into account accurate calculations of chemical rates, thermodynamic, transport and radiation quantities has been developed.Moreover, the fluid model has been coupled to a MCF code for electron kinetics.Neutral species mole fractions, gas temperature, electron number density and electron temperature, calculated from the model, have been validated against spatially-resolved measurements obtained by means of Thomson and Raman laser scattering diagnostics.
Measurements and simulations show, with good agreement between them, that, as pressure increases, T g increases in the plasma core from 4900 K to 7200 K, hence the neutral species composition becomes increasingly more dissociated due to thermal chemistry.Both n e and n e /n increase in the plasma core with pressure, from 2 × 10 19 m −3 at 100 mbar to 6 × 10 19 m −3 at 400 mbar and from 2.5 × 10 −5 at 100 mbar to 2 × 10 −4 at 400 mbar, respectively.Moreover, T e decreases from 1.5 eV to 0.6 eV in the plasma core, with increasing pressure.
Results of the model show that, as pressure increases above 100 mbar, sharper radial gas temperature gradients are formed.Moreover, at higher pressures, radial diffusion of electrons and electron attachment to CO 2 and O 2 molecules are largely suppressed.Hence, for pressures above 100 mbar, charge losses are uniquely determined by electron recombination with CO + and O + 2 ions.As a consequence of the high T g , molecular ions tend to dissociate into atomic species in the center of the tube, thus the role of electron-ion recombination decreases in the core, while being maintained at the edges.The present results are in agreement with previous studies by Martinez et al [33], Ridenti et al [34], and Viegas et al [23] confirming that plasma contraction can be explained by non-uniform gas heating.In addition, this work highlights the role of associative ionization as dominant electron production mechanism in the plasma core, as pressure increases above 100 mbar.This mechanism is favoured by the high T g and by the consequent dissociation of CO into C and O atoms and it leads to an increase of core n e and a decrease of core T e , with increasing pressure.Thus, a self-reinforcing cycle between n e and heating in the core is established.Results of this work confirm the hypothesis by Wolf et al [15] of the importance of the coupling between thermal-ionization instability and plasma chemical kinetics in driving contraction in CO 2 MW plasmas.Nevertheless, as opposed to the conventional thermal-ionization instability based on electron impact ionization described by Zhong et al [104,105], here the instability is driven by associative ionization reactions.
This work has shown that the output of the fluid-MCF model is also largely sensitive to the choice of adjustable quantities, such as the input MW power density and associative ionization rate coefficients.In this respect, this work can be considered as a step forward with respect to 0D plasma chemistry models towards a systematic validation of multidimensional model results against experiments for MW discharges in CO 2 .

Figure 2 .
Figure 2. Visual appearance of CO 2 MW plasma detected using a Charged-Coupled Device (CCD) camera in diffuse (left) and contracted (right) regime, at 100 mbar and 250 mbar, respectively.

Figure 3 .
Figure 3. Workflow of the fluid-MCF model.See text for a detailed description of the iterative procedure.Note that MCF is run locally for each point along the radial coordinate.

Figure 6 .
Figure 6.Thermal dissociation rates for CO 2 as a function of the radial position, at (a) 100 mbar and (b) 150 mbar.

Figure 8 .
Figure 8. Measured intensity profile from atomic oxygen emission (O(3p 5 P) → O(3s 5 S)) (solid line) and calculated electron density profile (dashed line) at (a) 100 mbar and (b) 250 mbar.The profiles have been normalized to one in correspondence to the peak, at r = 0.

Figure 9 .
Figure9.Measured electron number density[25] (dots) and calculated electron and ion number densities (lines) in the center of the tube, for different pressures.

Figure 10 .
Figure 10.Radial profile of measured electron temperature from Thomson scattering [25] (dots), calculated electron temperature as T e = 2 /3 (solid line), and Thomson temperature T e,Th derived from fitting of the EEDFs with a Maxwell-Boltzmann distribution in the 0-2 eV energy interval (dashed line), for (a) 100 mbar and (b) the 250 mbar.

Figure 11 .
Figure 11.Calculated mean square values of the reduced electric field in the center of the tube, for different pressures.

Figure 12 .
Figure12.Measured[25] (dots) and calculated (lines) electron temperatures (black) and gas temperatures (red) in the center of the tube, as a function of pressure.

Figure 13 .
Figure 13.Reaction rates for the main electron production mechanisms calculated at (a) 100 mbar and (b) 250 mbar, as a function of the radial coordinate.

Figure 14 .
Figure 14.Main electron loss rates calculated at (a) 100 mbar and (b) 250 mbar, as a function of the radial coordinate.

Figure 15 .
Figure 15.Measured (dots) and calculated (lines) (a) number densities and (b) electron temperatures at r = 0, as a function of pressure.Effect of uncertainties in the peak absorbed MW power density (P in (0)) is investigated by comparing results of the original model (black line), with the ones obtained by increasing (red line) or decreasing (blue line) P in (0) by 50%.

4. 4 .
Possible sources of uncertainties 4.4.1.Effect of absorbed microwave power density.The MW power density profile (P in (r)) used as input to the model (equation (

Figure 16 .
Figure 16.Rate coefficient for associative ionization of C + O (k C-O ai ) used in this work (solid black line), for O + O (k O-O ai ) reported by Park et al [96] (dashed blue line), and for N + O (k N-O ai ) reported by Popov [107] (dashed-dotted red line).

Figure 17 .
Figure 17.Measured (dots) and calculated (lines) (a) electron number densities and (b) electron temperatures at r = 0 as a function of pressure.Effect of uncertainties in the associative ionization rate coefficient for C + O is investigated by comparing results of the original model (black line), with the ones obtained by using different expression for this rate coefficient, as given by Park [96] (blue line) and by Popov [107] (red line).

Table 5 .
List of charged particle reactions: ion transfers, detachment, ion-ion recombination.Units for rate coefficients for two-body reactions are cm 3 s −1 and for three-body reactions are cm 6 s −1 .M represents any neutral species and T g is the gas temperature in K.

Table 6 .
List of associative ionization reactions.T g is the gas temperature in K.

•
[23]6], Park report an expression for the rate coefficient for associative ionization between O atoms (i.e.O + O → Note that this rate coefficient has been used in the model by Viegas et al[23]for both O + O and C + O mechanisms.