Charged-particle transport models for global models

This work analyses the influence of different charged-particle transport models on the global modeling of low-temperature plasmas. The simulations use the LisbOn KInetics simulation tool, calculating the charged-particle loss frequency due to transport with various formulations, categorized into two large groups: ambipolar-based and h-factor transport models. The models are applied to the description of (i) a DC discharge in oxygen (as example of an electronegative multi-ion plasma), at low to intermediate pressures, adopting a validated kinetic scheme as reference model; (ii) a microwave discharge in helium (as example of an electropositive multi-ion plasma), from low to atmospheric pressures, proposing and validating an updated kinetic scheme that bridges two previous models, adjusted separately for low to intermediate pressures and for atmospheric pressure. The analysis shows that using different charged-particle transport models can result in uncertainties of 20%–60% and 8%–115% in the discharge characteristics of oxygen and helium, respectively, with larger dispersion at low pressure and low electron density. The spreading in the results is observed also in the densities of the main plasma species, corresponding to uncertainties up to 60% and within 50%–150% in oxygen and helium, respectively. Since transport accounts for more than 50% of total charge losses, this mechanism should always be part of the quantitative sensitivity analysis of a kinetic scheme, considering several models with different validity domains, according to the electro-positive/-negative single-/multi-ion plasma under study and the low/high pressure conditions considered.


Introduction
Global (zero-dimensional) models are very popular in the modeling of low-temperature plasmas (LTPs), especially when the focus is on plasma chemistry. In these models, the dynamics of the plasma species s is described using an algebraic form for the spatial term of the rate-balance equations with n s the density of the species, S s its net creation rate due to the plasma chemistry, and ν transp s a spatial-averaged loss frequency due to transport to the walls. Often, the conservation equation (1) are reduced to a simple balance between gain/loss rates as defined by the chemical reactions, hence neglecting the transport term ν transp s n s [1], assuming uniform plasma. This procedure is because global models provide only a spatiallyaveraged description of the plasma, but this does not mean that the time-evolution of the species (average) densities is exclusively controlled by the chemistry and not influenced by transport effects. Unfortunately, and despite the existence of proposals in the literature to consider space-dependent phenomena in global models [2], the transport of species (in particular charged particles) is still often not considered and/or discussed in these models (namely in terms of validity limits), particularly when the study concerns a chemistry-oriented analysis.
The transport of charged particles is key to identify the operating point of a gas discharge, defined through a Schottkylike condition [1,3] that balances the net creation rate of electrons and ions due to kinetic mechanisms in the volume, with their loss rate to the walls. Satisfying this condition allows determining as eigenvalue the electric field that maintains the plasma, consistent with the required net creation-rate of charged particles, for given discharge parameters (e.g. gas density and plasma dimensions) that define the corresponding particle losses. The result can be represented graphically as discharge characteristics, by plotting the reduced electric field as a function of the reduced dimensions, both quantities normalized to the gas density. In the end, a direct effect of this eigenvalue on the plasma chemistry is expected, via the electron excitation and ionization frequencies influencing the kinetics. For this reason, the present work focuses on the study of charged-particle transport, very relevant to obtain accurate global model solutions.
The literature proposes several formulations to describe the transport of charged particles in the global modeling of LTPs. All formulations aim at obtaining a spatial-averaged loss-rate of charged species due to transport, which can be generally represented by the product of the species density by the corresponding loss frequency (see (1)). Furthermore, the different formulations can be classified into two broad categories, here termed as follows: (i) ambipolar-based models [4,5], where the loss frequency corresponds to the ratio between an ambipolar diffusion coefficient and a diffusion length. These quantities are estimated from the solution of the particle and momentum-transfer fluid equations, subject to suitable boundary conditions (usually zero-density at the walls), the latter equation corresponding to an ambipolar diffusion flux proportional to the density gradient; (ii) h-factor transport models [2,[6][7][8] that rely on scaling laws for the density gradient, via the ratio between the ion density at the sheath edge and at the plasma center (the so-called h-factor). This ratio is used to write the outward flux of positive ions at the sheath edge, as a function of the density at the center, the h-factor and the Bohm velocity at the sheath edge. The particle flux is finally used in the volume averaging of the ion-particle balance equation, to express the corresponding loss frequency in terms of the hfactor.
Both ambipolar-based and h-factor models can be easily implemented in global model calculations, for example using some of the most popular numerical codes available within the LTP community: ZDPlasKin, a freeware code developed by Pancheshnyi et al [9]; GlobalKin, developed by Kushner et al [10,11] and available upon request; Quantemol-P, developed with a graphical user interface [12], corresponding to a commercial application that extends GlobalKin to the automatic generation of a plasma chemistry from a set of userdefined species; and the global model within the PLASIMO commercial software, developed by van Dijk et al [13]. More details on global models and the codes used for solving them can be found in [1,14].
In this work, we analyze the influence of several chargedparticle transport models on the global modeling of LTPs. The study considers various formulations of ambipolar-based and h-factor transport models, which are applied to the modeling of (i) a DC discharge in oxygen, as example of an electronegative (EN) gas, at low to intermediate pressures, adopting a validated kinetic scheme as reference model; (ii) a microwave discharge in helium, as example of an electropositive (EP) gas, spanning low to atmospheric pressures, proposing and validating an updated kinetic scheme that bridges two previous models, adjusted separately for low to intermediate pressures and for atmospheric pressure. The global modeling of these discharges uses the LisbOn KInetics (LoKI) simulation tool [15][16][17], running for the different kinetics schemes and charged-particle transport models considered. Modeling results are compared with experimental data for the discharge characteristics and the densities of plasma species, quantifying the deviations between calculations and measurements and between simulation results obtained with different transport models.
The organization of this paper is as follows. Section 2 presents the general plasma configuration adopted and introduces the fundamental quantities used in the formulation of the transport models. Sections 3 and 4 review the formulations of the ambipolar-based and h-factor transport models, respectively, obtaining the expressions of the loss frequency due to transport to the walls ν transp , as proposed in the framework of these models for various electronegativities, pressures and single/multi-ion scenarios. Here, an effort is made to organize the different theories, paying particular attention to the physical assumptions and validity limits associated with the different formulations. Section 5 briefly introduces the LoKI simulation tool, summarizes the charged-particle transport models implemented therein, and presents and discusses the results obtained in oxygen and helium. Section 6 presents final remarks and concludes.

General configuration and notation
The following sections will briefly review several descriptions of charged-particle transport in LTPs, which are among the most popular choices in kinetic simulations that adopt global models. The transport descriptions span different types of plasmas and pressure regimes, in several cases bridging various limiting scenarios, from electro-positive (EP) to electro-negative (EN) plasmas, and from low-pressure (LP) to high-pressure (HP) regimes.
The usual configuration corresponds to a cylindrical plasma with radius R and length L, generally composed by electrons e and several positive ions i and negative ions i n (eventually clustered into single positive and negative species + and −, respectively), with densities n k , masses M k , temperatures T k , mobilities µ k , and free-diffusion coefficients D k (k = e, i, i n , +, −) generally defined as with k B the Boltzmann constant, and ν k = e/(M k µ k ) the ionneutral collision frequency (note that k B T k /e = D k /µ k ). The plasma is further characterized by the following quantities: representing the plasma electronegativity, which is defined at the center of the discharge as α 0 ≡ n −0 /n e0 ; the temperature ratios and the positive ion mean-free-path (k = i, +) with v th k the ion thermal velocity, defined within each transport model. Note that the electron temperature is defined here as T e ≡ 2/3ε (in eV), where ε is the electron mean-energy obtained by integration of the kinetic energy over the non-equilibrium electron energy distribution function (EEDF).

Ambipolar-based transport models
In dense plasmas, characterized by λ D ≪ Λ (with λ D the Debye length and Λ the characteristic length for species transport), the transport of charged particles takes place under ambipolar conditions, and can be studied using a generalized diffusion description. In this framework, the flux of each positive ion i writes where D amb i represents an ambipolar diffusion coefficient for ion-species i, satisfying the condition with D ambe an ambipolar diffusion coefficient for electrons.
The ion flux (3a) can be used in the corresponding continuity equation, to introduce an effective frequency for the losses of positive ions due to the transport to the walls where the characteristic length can be estimated assuming a separation of variables in the axial and radial directions, and further considering zero-density boundary conditions at the discharge walls (identified with the sheath edges), yielding Equations (4a) and (4b) correspond to an algebraic form of the transport loss-term in the ion fluid-equations, very convenient for global models. This general formulation can be realized in various models, specifying different ambipolar diffusion coefficients in order to allow a transition between various pressure and electronegativity regimes.

Effective ambipolar diffusion coefficient
The ambipolar diffusion coefficient D amb i can be defined for HP/LP regimes, by expressing the results of the unified theory proposed by Self and Ewald [18] in the convenient form of an abacus for an effective (ambipolar) diffusion coefficient D eff i as a function of pressure, as proposed by Ferreira and Ricard [19].
The theory derived in [18] provides a smooth transition between the classical ambipolar diffusion description [20] (valid at HP), and the free-fall description [21] (valid at LP), accounting for the ion inertia under the presence of the ambipolar electric field. The formulation considers the continuity and the momentum-transfer equations for a single type of positive ions, the momentum-transfer equation for electrons (neglecting its inertia term), and the quasi-neutrality condition, and it was applied to an EP plasma composed of electrons and a single type of positive ions.
The original concept of an effective diffusion coefficient was later extended by Coche et al [5] to plasmas with several positive ions and a single negative ion with low density (α ≲ 0.1), assuming that: (i) all positive ions have similar temperatures, drift velocities, ionization frequencies, and effective momentum-transfer collision frequencies; and (ii) all positive ions have approximately the same density profile, satisfying ⃗ ∇n i /n i ≃ ⃗ ∇n j /n j . The final abacus corresponds to the functional relation, for electrons, D effe /D ae vs. Λ/λ + , namely [5] assuming that the ratio D eff i /D a i is the same for all positive ions. Here, the clustering of positive ions introduces densityweighted averages for any quantity X, defined according to [the exception being the ion mass, for which one introduces the subscript 'Self' tags the results obtained in [18], and equation (5c) corresponds to the original abacus proposed in [19] and updated in [5]; D ae and D a i are the classical ambipolar diffusion coefficients for electrons and ions in a plasma, under the quasi-neutrality condition n + /n e = 1 + α, respectively [4] λ + is defined according to (2e) with v th+ = (3k B T + /M + ) 1/2 [22] corresponding to the root-mean-square velocity. Note that, when applied to a low-EN plasma (α ≲ 1), equations (5a)-(5c) are valid only at LP (Λ/λ + < 1), to ensure vanishingly small collision frequencies [5]. However, for EP plasmas, the equations are valid across the entire pressure range; indeed, in the HP limit of Λ/λ + ≫ 1, equations (5a)-(5c) yield D effe /D ae ≃ D eff i /D a i ∼ 1, corresponding to the classical ambipolar diffusion model.
Note finally that, similarly to equation (3b), the effective coefficients also satisfy the ambipolar condition

Ambipolar diffusion coefficient for EN plasmas
The theory of ambipolar diffusion can be generalized to account for the presence of negative ions, by locally balancing the flux of positive ions to the fluxes of electrons and negative ions, under the assumption of quasi-neutrality. The general formulation was developed by Rogoff [4], yielding the following ambipolar diffusion coefficients in a three component plasma which are used to describe the charged-particle transport adopting different approximations.
Equations (9a) and (9b) can be used in plasmas with multiple ions, if their fluxes are clustered into two equations for single positive + and negative − components. The formulation requires adopting the assumption of proportionality for each type of ion, ⃗ ∇n i /n i ≃ ⃗ ∇n j /n j and ⃗ ∇n in /n in ≃ ⃗ ∇n jn /n jn , and the density-weighted averages (6) to define D + and µ + as a function of the individual coefficients D k and µ k [4,5].
The theory described here was adopted in the QGM [25]. In QGM, the positive ions loss-frequency ν QGM transp i is written according to the formulation of GlobalKin [26], which considers not only the effects of the ambipolar transport (similarly to (4a)) but also the thermal losses at the wall, described by characteristic frequencies estimated from the volume average of the corresponding particle balance equation (cf section 4.1) where A = 2π RL + 2π R 2 is the total area of ion collection and V = π R 2 L is the plasma volume. Therefore, in this formulation the ion losses at the walls are due to ambipolar transport and thermal losses, with frequencies where γ ri is the recombination probability of the i positive ion. Finally, the total loss frequency is obtained from the total characteristic time for both mechanisms [25] ν QGM where in general we take γ r i = 1. The ambipolar diffusion coefficient in (11d) is given by (10a), with D + calculated from the density-weighted average (6). Therefore, the implementation of the model requires information about the free-diffusion coefficients D i for individual positive ions, which in this case are calculated from Here, the mean-free-path is given by with σ ij the momentum-transfer cross section for the collision of the ion i with species j (neutral, positive ion, or negative ion) and δ ij the Kronecker delta that filters out self-collisions; v th i corresponds to the ion mean velocity and the ion temperature is written from guidelines given in [8], In equation (12b), the momentum-transfer cross section is approximated with the results for the hard-sphere model in the case of ion-neutral collisions, and for Rutherford scattering [20] in the case of ion-ion collisions, respectively where σ LJ i is the Lennard-Jones potential for the i-species, and b 0 ij is the classical distance of closest approach given by with m R ij = m i m j /(m i + m j ) the reduced mass, and v R ij ≃ v th i the relative speed.

Eigenvalue solution for the diffusion length.
The ambipolar coefficients (9a)-(9c) can be simplified without any assumption regarding the charged-particle density profiles, contrary to what is proposed in section 3.2.1.
For moderate-to-small electronegativity α(µ +,− /µ e ) ≪ 1 (as termed in [4]), and by further considering µ +,− /µ e ≪ 1 together with the cold-ion approximation γ, γ − ≪ 1, equations (9a)-(9c) yield Equations (14a)-(14c) can replace D amb i in (3a), to write the fluxes of positive ions, negative ions and electrons, respectively, and the resulting expressions can then be used in the corresponding continuity equations, to be solved for adequate boundary conditions. Note that the negative-ion ambipolar coefficient is negative, thus yielding a negative flux that reveals, in this case, a dominant drift under the influence of the space-charge electric field.
This formulation was initially proposed by Ferreira et al [27], to describe the radial transport in a DC axiallyhomogeneous positive column with three types of charged particles (+, − and e, as in (14a)-(14c)), considering electronimpact ionization, attachment and detachment as creation/destruction kinetic mechanisms, and was later extended by Guerra and Loureiro [28] for application to an oxygen plasma with several positive ions i, where electron-ion recombination and charge-transfer reactions were further considered. Note that the original expressions (14a)-(14c) in [27,28] replace the electron temperature in γ and γ − by the corresponding characteristic energy u ke /e ≡ D e /µ e (in eV) [16], accounting for non-equilibrium situations, in which case γD + = e u ke µ + and The final set of equations corresponds to an overconstrained two-point boundary-value problem involving two eigenvalues, α 0 and λ eigen , the latter defined as [28] Here, ν I i represents an effective ionization frequency for the creation of positive ion i by direct and stepwise electron impact, as well as collisions between heavy species, such that the effective average rate is written as ν I i n e ; C r i is the electron-ion recombination rate coefficient; and Z ij is the charge-transfer average rate for collisions of ion i with a neutral j-species. In [28], the single positive-ion formulation of [27] is extended to a multi-ion scenario by multiplying both sides of (15a) by n e = ∑ i n i /(1 + α) λ eigen and by assuming the equality of all individual components in the summations, such that This somewhat arbitrary procedure [28] allows writing (15c) as diffusion equations similar to (4a) for the individual positive ions i, where we have taken D a i ≃ γ i D i (valid at low electronegativity, see (14a)) and introduced the effective diffusion length with λ eigen given by an expression fitting the results of figure  6 of [27] λ eigen = 2.
In (17a), P and Q are the reduced attachment and detachment frequencies, respectively, defined as with ν a and ν d the total frequencies for electron attachment and detachment, respectively; and the coefficients a 1 − a 3 and t 1 − t 3 are given in table 1. Equation (16b) reveals a dual behavior for Λ eigen as a function of the plasma electronegativity. Indeed, on the one hand, an increase in α yields a direct increase in Λ eigen , as evidenced by (16b), in order to reduce the losses of positive ions as to compensate the increase in the number of negative ions; but, on the other hand, λ eigen increases for larger attachment frequency P and smaller detachment frequency Q (see (17a)), to account for deviations from Bessel-shaped profiles for the chargedparticle densities, which yields a decrease in Λ eigen according to (16b). For smaller electronegativity, equations (16b) and (17a) give Λ eigen ≃ R/2.405, thus retrieving the classical ambipolar diffusion length for EP plasmas.

The h-factors transport models
The transport of charged-particles in LP discharges is often described in global models using scaling laws for the density gradient, that rely on the flux of positive ions at the sheath edge, expressed for a single-ion in an EP plasma as with n +s the positive-ion density at the sheath edge, and u B the Bohm speed. This flux can be written as a function of the edge-to-center density ratio (designated h-factor), and can then be used in the ion-particle balance equation to express the corresponding gain/loss frequency in terms of the h-factor. The expression of this loss frequency, initially obtained for LP/EP discharges, has been generalized to different pressure and electronegativity regimes, by extending/updating the original formulation. In relation to pressure, these extensions are often made by simply patching together the limiting results for LP/HP, which can introduce inconsistencies if these results are expressed as a function of different (though related) physical quantities.
In EP discharges, the positive-ion densities are assumed to have an uniform profile throughout the discharge except near the wall, where the density is assumed to drop sharply to a sheath-edge density. In EN discharges, the electron density is assumed to be uniform throughout the discharge except near the sheath edge; the negative-ion density is assumed to be parabolic, dropping to zero at the sheath edge; and the positive-ion density is n + = n e + n − , with n +s ≃ n es at the sheath edge.

h-factor formulation for the loss frequency due to transport
The effective frequency for the losses of positive ions due to the transport to the walls can be estimated from the volume average of the corresponding particle balance equation where ν I represents an effective ionization frequency; Γ s,axial = n +s,axial u B and Γ s,radial = n +s,radial u B are the positiveion fluxes at the edges of the sheaths, along the axial and the radial directions, respectively, with the Bohm speed given by and A L = 2π RL and A R = 2π R 2 are the ion-collection areas along the axial and the radial directions, respectively (note that A L + A R = A). Strictly, the electron density in the denominator of (19) should correspond to an average over the plasma volume, but in general its value at the discharge center is taken instead.
The previous relations can be used in (19) to yield where is the effective area for ion loss, featuring the edge-tocenter positive-ion density ratios along the axial and the radial directions, h L+ ≡ n +s,axial /n e0 and h R+ ≡ n +s,radial /n e0 , respectively.
The h-factors in (21) can be estimated by equating the flux of positive ions at the sheath edge (18) to its expression obtained from different formulations. The next sections briefly review some of the most popular models.

Drift-dominated models
Drift-dominated models assume that the ion drift-velocity due to the electric field dominates over the velocity due to the pressure gradient (which corresponds to assume the cold-ion approximation, setting T + = 0). The models can be organized according to their validity limits in terms of the LP/HP and EP/EN nature of the discharges.

LP solutions for EP discharges.
At LP, the ion driftvelocity ⃗ u + = µ + ⃗ E (with ⃗ E the space-charge electric field) can be written using a drift-dependent mobility [20] which assumes that ion-neutral encounters are dominated by charge-transfer collisions. In this case, by further assuming that the electrons are in Boltzmann equilibrium with the spacecharge potential, one obtains a non-linear form for the steadystate ion continuity equation, where the ionization frequency ν I is assumed to be homogeneous. The problem has been solved by Godyak [6] and by Lichtenberg et al [7] for the density profile in EP discharges, who have also extended heuristically the solution obtained to the collisionless limit obtained by Tonks and Langmuir [21], when the ion drift-velocity is governed by the energy conservation equation.
The results of this study, valid in the low to intermediate pressure (LP) regime where L/(2λ + ) ⩽ γ and R/λ + ⩽ γ, allowed obtaining the h-factors h Godyak for a plane-parallel discharge and an infinitely-long cylindrical discharge, respectively.

Analytical solution for plane-parallel EP discharges.
The particular case of a plane-parallel EP discharge can be solved analytically for a single positive-ion, if the momentumbalance equation describing the drift-velocity is written neglecting the contribution of charge-transfer collisions and by assuming a constant ion-neutral collision frequency for momentum-transfer ν + = u B /λ + . Note that this expression is equivalent to (2e), if the thermal velocity v th+ is identified with the Bohm speed. As in 4.2.1, the electrons are assumed in Boltzmann equilibrium with the space-charge potential, but here the ion momentum-balance equation further includes the corresponding inertia term. Czarnetzki and Alves [29] solved the resulting continuity and momentum-balance equations in Cartesian coordinates, obtaining the following implicit expression for the ion lossfrequency ν h−fact transp + , due to transport across any x, y, z direction of a cubic discharge with edge length L where β = 1, 1/2, 1/3 for transport along one, two and three (symmetric) Cartesian dimensions, respectively. Explicit inverted solutions can only be obtained in the LP (ν ≫ 1) and the HP (ν ≪ 1) limiting cases, which motivates an heuristic formula spanning the whole pressure range. For the one-dimensional situation the bridging formula is [29] yielding the following h-factors in the LP/HP limits (see (21)) for L ≫ λ + (HP), corresponding to a factor π/2 ∼ 1.57 difference with respect to the classical HP situation (see (28b) below).
Note that the extension of this formulation to cylindrical coordinates and/or a multi-ion scenario prevents the analytical resolution sketched above [30].

Drift-diffusion models for EP/EN discharges
The drift-dominated models in section 4.2.1 are limited to LP scenarios, with negligible pressure gradients. These gradients should be additionally considered when extending the models to HP, where diffusion phenomena become dominant.

4.3.1.
Transition from EP to moderate EN discharges. Lee and Lieberman [8] generalized the theory of section 4.2.1 to include transitions from EP to EN regions, also for a HP regime. The generalization was based in the one-dimensional analytical oxygen model developed by Lichtenberg et al [7], for HP discharges with plane-parallel geometry, where three different plasma regions were identified: • A EN region, characterized by a negative-ion density n − ≫ n e , confined to the center of the discharge and exhibiting a parabolic profile; • Two EP regions, with n − ≪ n e , which develop between the EN region and the two sheath edges; • Two thin sheath regions, where n + ≫ n e and n − ∼ 0.
In this case, the relation between the width of the EN region, 2l, and the chamber length L = 2l p , is found to be with γD + /(u B l p ) ≪ 1 at HP, which gives l ≃ l p for α 0 ≫ 1; thus, at high electronegativity, the EN region fills the entire volume, where the electron density is approximately constant. Moreover h HP,parabola which is similar to the HP diffusion solution for a planeparallel EP discharge (thus for α 0 = 0, with a cosine density distribution) h HP,EP,cosine Heuristically matching (24a), (28a), and (28b), one obtains a general h L+ factor that can be used for transitions from LP to HP, and from EP to EN regimes which indeed allows retrieving (24a) for LP/EP discharges (at α 0 = 0 and γD + /(u B L) ≫ 1), and a mixed form of (28a) and (28b) for HP/EN discharges (with α 0 ̸ = 0 and γD + /(u B L) ≪ 1). A similar treatment can be performed in the radial direction, to obtain a general h R+ factor with J 1 the first-order Bessel function, and where one retrieves (24b) for LP/EP discharges (at α 0 = 0 and We underline that the expressions (29) and (30) are valid if there is only one type of positive ion present. For multiple ions, (they are) valid if we make the assumption that the ratio of the sheath edge density to the bulk density is independent of the type of ion [8].

4.3.2.
Transition from EP to strong EN discharges. Usually, the model in section 4.3.1 is applied to the one-dimensional description of plasmas with moderate electronegativity. Kim et al [31] have attempted to generalize these results to a planeparallel plasma at a broader range of electronegativities and pressures, combining three EN discharge regimes (a) A LP regime, with a parabolic-profile EN core and two EP edges; (b) A LP regime at high α 0 , with a one-region parabolic EN plasma; (c) A HP regime at high α 0 , with a one-region flat-topped EN plasma.
The single description that covers the full parameter range (a)-(c) is obtained by stitching together [31] the edge-tocenter density ratios, using a linear ansatz that sums the h L results for the different regimes yielding: Here, where v th+ corresponds to the ion mean velocity (12c), and k rec denotes the (effective) positive-negative ion recombination coefficient calculated as with k recr the rate coefficient of recombination reaction r, and n (r) k the densities of the corresponding reactants. In principle, equations (31a)-(32c) are for a single type of positive ions and a single type of negative ions. However, Kim et al [31] generalize the formulation to multiple positive or negative ion species, introducing density-weighted averages for the ion temperature, the ion mean-free-path, the ion mass and the recombination coefficient, according to the rules described by (6).
Note that , yielding the results of Godyak [6] for a LP and EP plasma in a onedimensional slab geometry; In this case, n +0 = (1 + α 0 )n e0 increases in the EN core, but the loss flux in the EP edge remains nearly constant, accounting for the decrease of h Kim L+ with (1 + α 0 ) −1 (particularly of its h a component). As α 0 increases, the central EN core fills most of the plasma, such that the loss flux is governed primarily by a diffusion solution in an EN plasma with a parabolic profile; In this case, the central EN core also begins to flatten because the positive-negative ion recombination loss competes with the positive ion diffusion, to make the particle balance more local, leading to a steepening of the edge gradient. At HP, these effects are observed at larger EN, to ensure • the linear ansatz (31a) was verified by Monahan and Turner [32], using particle-in-cell simulations.

Revised formulation at LP and extension to cylindrical
geometry. The effects of electronegativity in chlorine discharges (diluted with oxygen) were investigated by Thorsteinsson and Gudmundsson [33,34] adopting a modified form of the edge-to-center positive-ion density-ratio (31a) proposed by Kim et al [31]. They have join the solutions (31b) and (31d), for the LP parabolic-profile EN core with EP edges (h a ), and for the HP at high α 0 one-region flat-topped EN profile (h c ), proposing the general ansatz h 2 = h 2 a + h 2 c [35] and generalizing the formulation to the axial and radial scaling factors for each Note that the original expressions (31b) and (31d), for h a and h c respectively, were modified: (i) to account for diffusion at HP, by incorporating the diffusion terms of (29) and (30) [8]; (ii) by replacing γD + with the more general expression of the ambipolar diffusion coefficient [34] (see (10a)) with α s the electronegativity at the sheath edge; and (iii) by replacing in (31d) the densities n +0 and n −0 at the center of the discharge, with the corresponding average quantities n + and n − , respectively. Note further that, according to [34], equations (33a)-(33d) and (34) are written for each positive-ion i, using individual ambipolar diffusion coefficients D a i (instead of considering the clustered description of the original theory, see section 3.2), electron-to-ion temperature ratios γ i ≡ T e /T i , mean-free-paths λ i (defined similarly to (12b) as λ −1 i ≃ ∑ j Nσ ij , adopting constant momentum-transfer cross sections σ ij ∼ 10 −19 − 10 −17 m 2 [33]), and ion sound-speeds u B i = √ k B T e /M i . A more general expression for the ion soundspeed was adopted in the previous paper [33], considering the effect of the negative ions in the structure of the plasma-sheath: . Instead, here the scaling factors (33a)-(33e) capture the modification to the Bohm speed due to the presence of negative ions near the sheath edge [31,34].
Note finally that equations (33a)-(33d) and (34) correspond to a typical patch of the limiting results for LP/HP, where the LP result is expressed using the ion mean-free-path λ i , whereas the HP result uses the ion diffusion coefficient D i , which in general is directly obtained from the literature rather than being estimated from the elementary quantity λ i (or σ ij ).
Equations (33a)-(34) require information on the plasma electronegativity at the discharge center α 0 , and at the sheath edge α s . The (total) relative density of negative ions at the discharge center can be calculated as a function of the average relative density α assuming a parabolic profile α( In weakly EN discharges α s ∼ 0, but in strongly EN discharges the sheath electronegativity can be comparable to the electronegativity at the plasma-presheath edge α b , according to [31] The numerical solution of equation (36) can be fit to a function of α s and γ − (for γ − > 10) as [34] with α b ≃ α and the coefficients a 1 − a 6 given in table 2. Note that equations (37a) and (37b) may yield negative results for α s /α b , which should be set to zero.
The results of Kim et al [31] were also revisited by Chabert [36] to propose a corrected h L -formula for LP and EN discharges, using as starting point equation (24a) for LP and EP discharges. The formulation assumes cold positive ions, Boltzmann electrons and Boltzmann negative ions to model a steady-state, planar (one-dimensional) EN discharge, by solving the positive ion particle and momentum equations, together with Poisson's equation.
The result yields where the second term is similar to that in (29), and where the new factor (1 + α 0 ) 1/2 in the denominator of (38) ensures a quicker transition to the constant collision frequency regime as α 0 increases. As intended, the corrected equation (38) allows retrieving (24a) at low EN and for LP conditions, i.e. for α 0 ≪ 1 and L/λ + ≪ 1. Note that (38) is not equal to the positive ion edge-to-center density ratio, since the ion velocity at the sheath edge is not equal to u B (equation (20)) when α 0 ̸ = 0. Note further that the fluid model in [36] is formulated for a single positive ion and a single negative ion. Using (38) for multi-ion scenarios is either nonapplicable, or subjected to the assumption that the ions have similar edge-to-center density ratios. There is general good agreement between (38) and (31a), at various pressures and electronegativities [36], but the following should be noted: • Equation (38) is more accurate at very LPs (L/λ + ≪ 1), since it was deduced to bridge the EP and the EN regimes, following (24a) and (29). Indeed, in the LP limit the results shown in [36] are notably different, even if the solution of (31a) does not diverge at low L/λ + , contrary to what is claimed in [36]; • The fluid model used in [36] assumes Boltzmann negative ions, which is valid only at LP, when the negative ions are not rapidly destroyed (either by recombination with positive ions or by detachment on neutral excited states), at the location where they have been created (by attachment). Therefore, the model is appropriated for LP discharges under the various plasma-density regimes of interest; for high EN and HP discharges, it is recommended using equation (31a).
In the latter regimes, the dynamics of the negative ions is governed by internal processes, and the h L factor becomes a function of the negative-ion density [36].

Results and discussion
The models presented in the previous sections were implemented in the simulation tool LoKI [15] to describe the chargedparticle transport and losses at the walls, in EP/EN plasmas for various LP/HP regimes. LoKI couples two main calculation blocks: a Boltzmann solver (LoKI-B) for the electron Boltzmann equation [16,17], released as open-source code under the GNU general public license [37]; and a Chemical solver (LoKI-C) for the global kinetic models of pure gases or gaseous mixtures [38].
The coupled solution of LoKI-B and LoKI-C is as follows. For a given set of discharge parameters (gas pressure p and temperature T g ; discharge dimensions R and L; reduced electric field E/N and reduced excitation frequency ω/N, where N = p/(k B T g ) is the gas density; electron density n e ; and gas mixture composition), LoKI-B calculates the non-equilibrium EEDF and the corresponding electron rate coefficients and transport parameters, using complete sets of validated electron-impact cross sections, obtained from the IST-Lisbon database with the open-access website LXCat [39]. The electron parameters are calculated by appropriate integrations of the electron cross sections over the EEDF, and are used in LoKI-C to obtain the self-consistent solution of the plasma chemistry. LoKI-C solves the zerodimensional (volume average) rate balance equations for the most relevant charged and neutral species in the plasma, and for the gas temperature (if a gas-heating thermal model is activated), knowing the kinetic scheme for the gas/plasma system. The initial gas pressure is then adjusted, considering the effects of dissociation/recombination mechanisms, in order to obtain the prescribed working condition p at steady state. A test upon charged-particle neutrality calculates a new E/N value, which is introduced in LoKI-B along with the new concentrations of the different heavy-species, in an iterative procedure that continues until convergence. Typical criteria impose relative errors of ∼10 −5 − 10 −2 for the different convergence cycles. As output, LoKI-C provides the densities of species, the gas temperature, the creation/destruction rates for the different reactions, and the self-consistent value of the reduced maintenance electric field. The coupled solution of LoKI-B and LoKI-C allows taking into account the impact of the plasma chemistry and the maintenance electric-field on the electron energy distribution.
Each transport model was used (when applicable) in the global modeling of helium plasmas (as example of an EP system) and oxygen plasmas (as example of an EN system).
For each system, the simulations started from a validated kinetic scheme, taken as reference, and were further pursued by checking the effects in the results of adopting different transport models.
The results of the global model ensure the balance between the net creation of charged particles due to kinetic mechanisms in the volume, and their destruction due to transport and losses at the walls, for example as formulated by equations (4a) or (19). These equations correspond to generalized forms of the Schottky condition, and thus can be formally expressed as functional relationships of the reduced electric field with the discharge parameters [40] or graphically as representations of E/N vs. NΛ, termed discharge characteristics, parameterized in the remainder working conditions. As mentioned, the reduced electric-field is self-consistently calculated by the global model to satisfy charged-particle neutrality. For given electron density (or discharge current), the variation of E/N influences the creation-destruction rates of the charged particles (by changing the values of the electron parameters via the non-equilibrium EEDF), allowing the equality in the corresponding densities Therefore, it is expected that the self-consistent E/N values (i) are sensitive to variations in the ion losses due to transport, which can be captured directly by changes in the discharge characteristics calculated for various transport models; (ii) are responsible for changes in the densities of species, which are related indirectly with modifications in the loss rates due to transport.

Summary of the transport models implemented in LokI
The transport models presented in the previous sections are implemented in LoKI as described below. For ambipolar-based transport models, the loss frequency ν amb transp i is calculated from the general expressions (4a) and (4b), with where D eff i is obtained from (5a) to (5c). This model is valid for multiple positive ions and a single negative ion with low density. Note that for EP/HP plasmas, with D eff i ≃ D a i and Λ given by (4b), one retrieves the classical ambipolar diffusion model (cf (7a) and (7b)); where Λ eigen is obtained from (16b) considering (17a)-(17c). This model can be applied to describe the radial transport of charged particles in moderate EN plasmas with multiple positive ions and a single negative ion; or is identified with ν QGM transp i given by (11d), for the QGM (see section 3.2.1), with D a+ obtained from (10a), considering (6) and (12a)-(13b). This model is valid for multiple positive and negative ions.
For h-factors transport models, the loss frequency ν h−fact transp i is generally given by (21), when explicit expressions for h L and h R are available. In this case, we adopt (see section 4.3.3) • The LP h-factor model, summarized in Annušová et al [41], which corresponds to the following loss frequency for each positive ion calculated adopting: equation (33b) for the edge-to-center density-ratio along the radial direction, taking (33d), (33e), (34), (35) and (37a)-(37b); and equation (38) for the edgeto-center density-ratio along the axial direction, more accurate at LP for the various EP/EN plasma-density regimes; • The alternative HP/EN h-factor model, preferred at HP and for high EN plasmas, according to the recommendations of [36], Equations (41a) and (41b) were used in plasmas with multiple positive and negative ions. Note that, in the HP and low EN limit, both frequencies (41a) and (41b) yield , (42) differing from the result of the classical ambipolar diffusion for EP plasmas, in a scenario of multiple positive ions, by factors of √ π/2 ∼ 1.25 and √ 2.405/[2 J 1 (2.405)] ∼ 1.52, for the axial and the radial diffusion lengths Λ L and Λ R , respectively.
Finally, we consider also the results of the analytical EP h-factor model of section 4.2.2, using the one-dimensional bridging formula (26) (corresponding to a situation where R → ∞), to obtain ν h−fact transp + for an EP plasma with a single positive-ion. Here, and similarly to the other h-factors models, the individual mean-free-paths λ i are defined similarly to (12b), adopting constant momentum-transfer cross sections σ ij ∼ 10 −18 m 2 .
In the next sections, we will apply these models to oxygen and helium, assuming the same temperature for all heavy species (i.e. T + ≃ T − ≃ T g ), corresponding to γ ≃ γ − . Moreover, we will adopt the same type of elementary data (mean-freepaths, collision cross sections and/or collision frequencies) mentioned in the references where the models were originally introduced. Simulation results are expected to depend on these data, and in some cases we will demonstrate and quantify this dependence.

Results for oxygen
The simulations in pure oxygen plasmas consider a cylindrical DC glow discharge (ω = 0) with radius R = 1 cm and length L = 52.5 cm, produced at low-to-intermediate pressures p ∼ 25-1300 Pa, for a low gas flow of 2-10 sccm and discharge currents of 10, 20, 30, 40 mA, corresponding to electron densities and gas temperatures in the ranges n e ∼ 10 15 − 8 × 10 15 m −3 and T g ∼ 320-580 K, respectively.
The model adopts the reaction mechanism presented and discussed in [42], for the ground-states O 2 (X,v = 0-41), O( 3 P) and O 3 (X); the electronic excited states O 2 (a 1 ∆), ; the positive ions O + 2 and O + ; and the negative ion O − . In [42], the simulations are compared with the measurements of Booth et al [43][44][45] for the different working conditions considered, and the reader should refer to this paper for details about the model validation (unfortunately, experimental uncertainties are not disclosed/quantified in these papers). We highlight the relevance in the simulation results of the O( 3 P) wall recombination probability γ[O( 3 P)], which depends on the wall temperature, the gas pressure and the discharge current. Following [42], we have considered here the values of γ[O( 3 P)] obtained from the measurements of Booth et al [44], for the corresponding working conditions, which are found in the ranges 3.0 × 10 −3 − 7.5 × 10 −4 (for p = 30-500 Pa) and 8.7 × 10 −3 − 1.7 × 10 −3 (for p = 30-1300 Pa), for discharge currents of 10 mA and 40 mA, respectively. For the purpose of the present work, the reaction mechanism of [42], which considers the eigenvalue diffusion length model to describe the charged-particle transport, is taken as reference. Figures 1(a) and (b) shows the discharge characteristics of E/N vs. NΛ, calculated in oxygen for the working conditions indicated above, at discharge currents I dc = 10, 40 mA, respectively, adopting the following charged-particle transport models: effective diffusion with D eff i = D a i , corresponding to classical ambipolar diffusion; QGM; and LP h-factor, taking the same momentum-transfer cross section σ ij = 10 −18 m 2 for all i − j collisions; in addition to the eigenvalue diffusion length reference model. For comparison, this figure also shows the experimental values reported in [42,43] (the working pressure was converted into the gas density N, using the ideal gas law with the measured values of T g ). The simulations involve the self-consistent calculation of (i) the electron density, as to obtain the input value of the discharge current by using I dc ≃ en e v d π R 2 , with v d = µ e E the electron drift-velocity also calculated by the model; (ii) the average gas temperature, from the solution to a gas-heating thermal model [42] considering a wall temperature of 320 K according to measurements.
Note that some of the transport models presented in sections 3 and 4 could not be adopted here, namely the analytical EP h-factor which is limited to EP plasmas (with a single ion and slab geometry), and the effective diffusion applicable only to plasmas with low EN. Indeed, for the working conditions adopted here, the plasma electronegativity is found always above 0.1 (reaching more than 0.8 at LP and small discharge currents), as can be confirmed from figures 2(a) and (b).  [42,43], and calculated (curves) for a DC discharge with R = 1 cm and L = 52.5 cm at I dc = 10 mA (a) and 40 mA (b), using the following charged-particle transport models: classical ambipolar diffusion (solid), eigenvalue diffusion length (dashed), LP h-factor (dotted), and QGM (dash-dotted).
An observation of figure 1 shows a dispersion in the values of E/N between 20% and 60% (being larger at lower discharge currents), further revealing that the LP predictions of the LP h-factor model are closer to the measurements than those of the eigenvalue diffusion length reference model. Moreover, the predictions using ambipolar-based transport models are systematically overestimated with respect to the measurements, whereas the opposite is observed for the h-factors transport model. The exception is the ambipolar-based QGM, which also underestimates the measurements. Recall that QGM follows a 'mixed' formulation that estimates characteristics frequencies considering ambipolar transport and thermal losses at the wall, which are combined by taking the volume average of the particle balance equation, as in the h-factors models (see section 3.2.1). Note finally that even for the maximum pressure considered of~1300 Pa the ambipolar-based models and the h-factors model yield different E/N results; this was predicted in section 5.1 (for the case of EP discharges), but the differences are likely to be attenuate at HPs (see section 5.3). One might question if the dispersion found in figure 1 for the E/N values, when using different charged-particle transport models, translates into relevant changes in the EEDF, calculated by solving the electron Boltzmann equation, and the densities of species, calculated from the chemistry model. Figure 3 shows the EEDFs calculated at I dc = 10 mA and p = 133 Pa (i.e. NΛ ∼ 10 20 m −2 ), for the E/N values obtained using the different charged-particle transport models considered here. As expected, the non-equilibrium EEDFs are directly affected by the reduced electric fields, whose increase leads to an enhancement in the tail of the distributions. The influence of transport in the kinetics of oxygen ions is confirmed in figures 4(a) and (b) that shows, as a function of NΛ, the relative total ion losses due to the transport of O + and O + 2 , calculated for the same working conditions as before, using the different charged-particle transport models considered here. As observed, the destruction of positive ions is clearly controlled by transport, which accounts for more than 70%-90% of the total losses, for high and low discharges currents, respectively, whatever the working pressure. The large influence of charged-particle transport on the kinetics of oxygen plasmas extends also to neutral species. Simulation tests reveal that the use of different charged-particle transport models may lead to uncertainties of 60% in the densities of O 2 (a 1 ∆), O 2 (b 1 Σ + g ) and O( 3 P). Note that a more detailed analysis of the ion kinetic-paths reveals that the transport accounts essentially for the losses of O + 2 , not O + . This is because the atomic ion is rapidly converted into the molecular ion, mainly via charge-transfer with the oxygen ground-state, according to O + + O 2 (X,v = 0) → O + 2 + O( 3 P) that becomes the dominant loss mechanism for O + . Indeed, transport is responsible for >90% of O + 2 losses, while charge-transfer accounts for 40%-95% of O + destruction (between LP and HP whatever the value of I dc ), transport being the second most important loss mechanism of the atomic ion. For this reason, the population of the atomic ion remains 2-4 orders of magnitude below that of the molecular ion, making O + 2 responsible for over 95% of the total positive charge in oxygen. This can be observed in figures 5(a) and (b) that plots n O + 2 /n + × 100 (with n + ≡ n O + 2 + n O + ), as a function of NΛ, calculated for the same working conditions as before, using the different charged-particle transport models considered here. The results in this figure do not mean that one can ignore O + in the oxygen kinetics, they simply confirm that the molecular ion accounts for most of the positive charge in oxygen.
In figure 1, the results of the LP h-factor model were obtained for σ ij = 10 −18 m 2 , since the original [33,34], that regard the modeling of chlorine-oxygen mixtures, provide general estimations for these cross sections, which are reasonably found in the range 10 −19 − 10 −17 m 2 . As an example for I dc = 10 mA, figure 6 illustrates the influence of this parameter in the results of the LP h-factor model for the discharge characteristics, and it further compares the predictions of the models LP h-factor and HP/EN h-factor with varying pressure.
The results in figure 6 show that larger collision cross sections are responsible for smaller reduced maintenance fields (due to smaller transport losses), yielding a dispersion in the E/N values between 2% and 20%, which narrows down to 2%-12% for larger discharge currents (not shown here). For the σ ij values considered here, the calculated discharge characteristics remain below the measured ones, but this result is highly dependent of the collisional data. Moreover, the changes in the collisional cross section imply a modification of the ion mean-free-path λ i only in the LP terms of equations (33c) and (33d), where the corresponding HP terms remain unchanged. As mentioned in section 4.3.3, this is because the value of the ion diffusion coefficient D i , used in equation (34), is obtained from the literature, instead of being coherently estimated from λ i , which explains the convergence of the h-factors model results at HP. Figure 6 also shows that equations (41a) for LP h-factor and (41b) for HP/EN h-factor yield coincident results (the curves for the two models are indistinguishable), which is not surprising since L ≫ R for the current working conditions.
In this case, and contrary to oxygen, there is no single reference paper that we can easily adopt as starting point for the simulations in this work, especially when aiming the modeling of both LP/HP regimes. Indeed, although our previous models [46,47] are very similar, they present relevant differences in the set of species considered (as [47] further includes excimer species), in the electron-impact cross sections adopted (which are extensively updated in [47]), and in the rate coefficients of the associative ionization (AI) mechanisms, which are expected to depend on the gas temperature (although [46,47] present no variation law in T g for these rate coefficients, which are given as constant values instead).
Therefore, in this work we have started by reviewing the models of [46,47], to propose (and validate) a unified kinetic description for helium plasmas that bridges the LP (lower T g ) and HP (higher T g ) limits. The resulting reaction mechanism is essentially that described in [47], for both electron and heavy-species collisions. In particular, we have used the electron-scattering data proposed in [47], already available in the IST-Lisbon database on LXCat [39] for direct collisions under the datagroup He (the data for collisions with excited states is available as supplementary material in [47], and will also be published in LXCat). However, the present work introduces the following updates with respect to [46,47]: • A typo was found in table 2 of [47], in the rate coefficient of reaction (12b) for the dissociation of He 2 + , where the exponent of Tg(K) should be +0.67 instead of −0.6. Moreover, the following typos were found in the supplementary material with [47]: (i) the collision strength Ω(x) for dipole-forbidden transitions (equation (2b)) should be further multiplied by x 2 ; (ii) in equation (5a) for G(x), the square parenthesis [.] should not include the expression ln(x + 0.2). This error, which is already present in [48], can significantly change the shape of the cross section; (iii) in table 2, the electron excitation cross section from He(2 1 S) to He(3 3 S) is identified as σ 2 3 S 2 1 S instead of σ 3 3 S 2 1 S . • The set of heavy-species includes all the He(n 2s+1 l) states with principal quantum number n ⩽ 7 (l and s are the orbital and the spin quantum numbers, respectively), the atomic and molecular ions He + and He + 2 , and the excimer state He ⋆ 2 . Here, and contrary to what is done in [46,47] where the states He(7 2s+1 l) are considered as loss reservoir for the remaining states, we have defined an ensemble of kinetic mechanisms to describe their creation/destruction, similar to that of states He(6 2s+1 l); • Compared to the LP scheme of [46], several mechanisms were extended/revised: (i) Penning ionization reactions are considered also for the He(2 1 P) state; (ii) the dissociative recombination rate coefficients are calculated from the corresponding electron-impact cross sections, and this mechanism is extended to the production of He(2 1 S), He(2 3 P) and He(2 1 P), in addition to He(2 3 S); (iii) the threebody recombination of atomic and molecular ions is also considered; • The rate coefficients for AI, He(1 1 S) + H(n 2s+1 l) → e + He + 2 , have been extensively revised to account for their variation with gas temperature and to avoid tuning with the exclusive aim of improving the agreement between simulations and measurements [46,47]. According to [49,50], it is expected that these rate coefficients vary significantly with T g , which can explain (at least partially) the differences in the values adopted in [46] at T g = 1000 K, and [47] at T g = 1800 K. However, because [47] adopts the AI rate coefficients proposed in [50] for T g = 2450 K, and because [46,47,50] have used these coefficients as modeling parameters, it is difficult to take further conclusions especially while attempting to unify the LP/HP models. In this work, we used the original AI cross sections calculated by Cohen [49] for the triplet states He(3 3 S,P,D) and He(4 3 S,P,D) to evaluate the corresponding rate coefficients by adequate integration over a Maxwellian distribution function for the neutral atoms at T g . These results are then used to obtain the coefficients for the reactions involving other states. A presents details about the calculation procedure, including the cross sections dat in [49]. Table 3 presents the results obtained at T g = 1000, 1800 K, which were adopted in this work without any adjustment; this table presents also the relative deviations with respect to the values proposed for the AI rate coefficients in [46,47]. • The charged-particle transport is described using the effective diffusion model, which yields the classical ambipolar solution in the HP limit (see section 3.1).
The predictions of the revised model, with the updates described above, are presented in figures 7 and 8, which show the discharge characteristics and the LP/HP excitation spectra, respectively. In this case, the discharge characteristics correspond to curves of E rms /N vs. NΛ, where E rms ≡ E/ √ 2 is the root-mean-square of the HF electric field E cos(ωt) and Λ ≃ Λ R (see (4b)). Note that, in [46], the discharge characteristics are plotted introducing an effective electric-field, defined as E e ≡ ν ce / √ ω 2 + ν 2 ce , with ν ce /N = 6.8 × 10 −14 m 3 s −1 a representative value of the electron-neutral momentum-transfer rate coefficient. Note that the reduced effective field E e /N represents the relevant quantity to describe the power absorption by plasma electrons, and in fact values of E rms /N = 1 − 10 4 Td translate into E e /N ∼ 1 − 200 Td, for the pressure range considered here. The excitation spectra are presented as Boltzmann plots of ln[n j g 0 /(Ng j )] ≃ −u j /(k B T excj ) vs. u j , where the u j , T excj and g j represent, respectively, the energies, the excitation temperatures and the statistical weights of the different excited states j. For comparison, these figures show also the experimental values reported in [46,47], as well as the calculation results obtained with the AI rate coefficients proposed in these references.
We conclude that the revised model (i) yields results similar to those obtained with the original AI rate coefficients; (ii) reproduces the experimental trend for the E rms /N vs. NΛ curves, underestimating the measured E rms /N values bỹ 30%-45% (see figure 7); (iii) predicts excited states densities in agreement with the measurements (see figures 8(a) and (b)), within average errors <4% at LP and <2% at HP, calculated for the excitation spectra as ∥T excj − T exp excj ∥/T exp excj × 100 (see figures 9(a) and (b)). These observations contribute to the validation of the revised model, which will be taken as reference to analyze the variations of the results with modifications in the description of charged-particle transport.
As in the previous section, we have calculated the discharge characteristics in helium adopting different charged-particle transport models: effective diffusion with D eff i = D a i , corresponding to classical ambipolar diffusion; eigenvalue diffusion length; QGM, adopting for He + 2 the same Lennard-Jones parameters as for He + ; LP h-factor with σ ij = 2.6 × 10 −19 m 2 for all i − j collisions, according to Phelps [51] for He + -He collisions at k B T g /e ∼ 0.1 eV (for T g = 1000 K); and the effective diffusion reference model. The calculations adopted the working conditions indicated above with the AI rate coefficients proposed in this work, and the results are shown in figure 10 together with the experimental values reported in [46].
An observation of figure 10 shows a dispersion in the values of E rms /N between 8% (at HP) and 115% (at LP). In general, the lowest results are given by QGM and the highest ones by the pure ambipolar models, similarly to what was already observed in oxygen (see figure 1). As expected, the classical ambipolar diffusion and the eigenvalue diffusion length models yield coincident results, since (16a) and (16b) for EP plasmas (with α = 0 and P = Q = 0) correspond to (4a) and (4b) at D amb i = D a i . Because the classical ambipolar diffusion coefficient D a i is larger than the effective diffusion coefficient D eff i , the classical ambipolar diffusion model enhances the transport losses, yielding E rms /N values larger than the ones predicted by the effective diffusion model, as seen in figure 10. Indeed, the calculations with classical ambipolar diffusion (or eigenvalue diffusion length) could not go below NΛ ∼ 2 × 10 19 m −2 , because of the extremely high (non-physical) E rms /N values obtained at LP with these models. At HP, the Rel. deviation compared to [46] (%) Rate coefficient at Tg = 1800 K (m 3 s −1 ) Rel. deviation compared to [47] (%) He ( ion-transport losses are strongly reduced (see also figures 13 and 14). In this case, all transport theories converge to classical ambipolar diffusion losses, described by the same diffusion coefficients D a i for an EP plasma, as can be confirmed from the ambipolar formulations presented in section 3, and despite the small factor differences highlighted in section 5.1 between the ambipolar-based models and the h-factors model. As in oxygen, the differences in the E/N values obtained with different charged-particle transport models (see figure 10) yield changes also in the EEDFs, as confirmed in figure 11 calculated at p = 158 Pa (i.e. NΛ ∼ 10 19 m 2 ). As observed, the EEDFs deviate from equilibrium, exhibiting a depletion that starts around the first excitation potential of helium (∼20 eV).
The results in figure 10 depend not only on the transport models adopted, but also on the data used in these models. Figure 12 illustrates the influence in the results of the molecular-ion reduced mobility µ He + 2 N, used in the effective diffusion model. In [47], the calculations at atmospheric pressure adopt µ He + 2 N = 2.6 × 10 21 V −1 m −1 s −1 , following the suggestion of Belmonte et al [50] that considers the fielddependent mobility at E e /N ∼ 0.4 Td, according to [52]. In figure 12, the effective diffusion calculations used either the µ He + 2 N value of [47] or a value~10 times higher, obtained from a temperature-dependent function fitting the results of [53] which yields µ He + 2 N = 5.1 × 10 22 V −1 m −1 s −1 at T g = 1000 K. The higher molecular-ion mobilities are confirmed by various authors [54][55][56][57], whose reported values reveal a dispersion in this parameter of at least 25%, partly due to the difficulty in identifying the specific He + 2 state detected in the measurements. As observed in figure 12, the calculations at higher mobility yield larger reduced maintenance fields due to the enhanced transport losses. In this case, the calculated E rms /Ns are 40% larger for intermediate values of NΛ, which allows better agreement with the measurements. Note that the results of the effective diffusion model become independent of the ion mobility both in the LP limit, where transport is the most relevant loss mechanism, and in the HP limit, where transport losses are strongly reduced.
As in oxygen, simulations results depend also on the ioncollision cross sections adopted in the LP h-factor model. Our numerical tests confirm that larger σ ij yield smaller reduced maintenance fields at LP and similar results at HP, as the changes in the collisional cross section affect only the LP terms of equations (33c) and (33d).
The relevance of transport losses in the kinetics of helium ions is confirmed in figures 13(a) and (b) which represents, as a function of NΛ, the relative losses of He + and He + 2 , respectively, due to their most relevant destruction channels, as predicted by the effective diffusion model. As observed, the destruction of positive ions is controlled by transport for NΛ ≲ 5 × 10 20 m −2 (corresponding to p ≲ 5000 Pa, where transport accounts for more than 50% of the total losses of ions), matching the LP region where the different transport theories are responsible for a dispersion in the discharge characteristics (see figure 10). In the HP region (for NΛ ≳ 10 21 m −2 ) the atomic ion is converted into the molecular ion by a charge-transfer reaction, and the molecular ion is destroyed by electron dissociative recombination reactions yielding either He(2 3 P) or He(2 1 S). The results presented in figures 13(a) and (b), obtained with the effective diffusion model, remain valid also for other transport models, namely when identifying the main mechanisms responsible for charge losses. In particular, figure 14 compares the predictions of the different transport models considered here, concerning the relative total ion losses due to the transport of He + and He + 2 . As shown, charge transport accounts for more than 50% of the total losses for NΛ ≲ 5 × 10 20 − Figure 10. Discharge characteristics of Erms/N vs. NΛ in helium, measured (points) [46], and calculated (curves) for the same conditions of figure 7 with the AI rate coefficients proposed in this work, using the following charged-particle transport models: effective diffusion (solid), classical ambipolar diffusion (shortdashed), eigenvalue diffusion length (dotted), LP h-factor for σ ij = 2.6 × 10 −19 m 2 (long-dashed), and QGM (dash-dotted). Short-dashed and dotted curves are coincident. 5 × 10 21 m −2 (corresponding to p ≲ 5 × 10 3 − 5 × 10 4 Pa), depending of the transport model adopted.
Note that the loss of charges at HP is exclusively due to the molecular ion, which then becomes the dominant ion, as already suggested by the results in figure 13(a). This is because the charge-transfer mechanism He + + 2He → He + 2 + He is responsible for the very effective conversion of the atomic ion into the molecular ion. As result, He + 2 accounts for at least 80% of the total positive charge when NΛ ≳ 5 × 10 20 − 10 21 m −2 , as can be observed in figure 15 that plots n He + 2 /n + × 100 (with n + ≡ n He + 2 + n He + ), as a function of NΛ for the different charged-particle transport models considered here. Again, the dispersion in the results of figure 15 (revealing uncertainties of 0.5%-150%, depending on the pressure) demonstrates the influence of the chargedparticle transport models on the plasma kinetic simulations, which are not exclusive of ion densities. In fact, numerical tests yield relative variations of the neutral species densities at LP between 60% and 500% (the largest spreading in the results being observed for the metastables He(2 1 S) and He(2 3 S) and for the He(6 2s+1 l) states), and at HP in the range 1%-50% (the largest spreading in the results being observed for the metastable states).
In helium, we also performed simulations adopting the analytical EP h-factor model described in section 4.2.2. In this case, we changed the working conditions adopted because the model is valid for Cartesian coordinates in EP plasmas with only one type of positive ion. In particular, we have considered a radio-frequency (ω/(2π) = 13.56 MHz) slab-plasma with L = 0.3 cm, at low electron-density n e = 5 × 10 16 m −3 , intermediate-to-HPs p ≳ 100 Pa, and the same gas temperature T g = 1000 K as in previous simulations. The choice of a radiofrequency plasma with lower density is related to the capacitive excitation. The choice of the pressure region is twofold: (i) because the charged-particle model is valid for a single ion, Figure 13. Relative ion losses for He + (a) and He + 2 (b), as a function of NΛ, calculated for the same working conditions of figure 7 using the effective diffusion model. The curves correspond to the following loss channels: in (a), transport (solid), 3-body charge-transfer He + + 2He → He + 2 + He (dashed); in (b), transport (solid), electron dissociative recombination e + He + 2 → He ⋆ + He (dashed).
we have limited the simulations to the HP region where He + 2 is the dominant ion (see figure 15); (ii) for the values p, T g considered, we have ω/N ≲ 10 −14 m 3 s −1 ≪ 10 −13 m 3 s −1 ∼ ν ce /N, which allows calculating the EEDF using the stationary form of the electron Boltzmann equation [17]. Figure 16(a) shows simulation results, as a function of NΛ, for the relative population of the molecular ion in helium, calculated using the analytical EP h-factor model neglecting the presence of the atomic ion, which in practice was done by removing its losses due to transport and by disregarding its influence in the transport of He + 2 . For comparison purposes, the figure shows also simulations obtained with the h-factor model, considering or disregarding He + in charged-particle transport mechanisms. An observation of Figure 14. Relative total ion losses due to transport in helium, as a function of NΛ, calculated for the same working conditions and charged-particle transport models of figure 10. The short-dotted curve corresponds to the effective diffusion model with the temperature-dependent He + 2 mobility of [53], as in figure 12. Solid, short-dashed and dotted curves are coincident. figure 16(a) reveals that the three model versions yield similar results for NΛ ≳ 5 × 10 20 m −2 (corresponding to p ≳ 6500 Pa), which confirms that the influence of the atomic ion in charge transport can indeed be neglected in this HP region. This conclusion is confirmed in figure 16(b), representing the Simulation results in helium, as a function of NΛ, for (a) the relative population of the molecular ion, and (b) the discharge characteristic, calculated for a slab-plasma at ω/(2π) = 13.56 MHz, L = 0.3 cm, Tg = 1000 K and ne = 5 × 10 16 m −3 , using the following charged-particle transport models: h-factor (solid), h-factor disregarding He + in charged-particle transport mechanisms (dotted), analytical EP h-factor disregarding He + in charged-particle transport mechanisms (dashed). discharge characteristics for the same working conditions and charged-particle models as in figure 16(a), where the simulations were limited to NΛ > 10 20 m −2 . As before: (i) the analytical EP h-factor and the h-factor models, both disregarding He + in charged-particle transport mechanisms, yield similar results; (ii) the analytical EP h-factor model, disregarding He + in charged-particle transport mechanisms, and the h-factor model, considering the presence of the atomic ion, yield similar results for NΛ > 5 × 10 20 m −2 . These observations contribute to validate the analytical formulation.

Final remarks
In this work, we have analyzed the influence of different charged-particle transport models on the global modeling of LTPs. The simulations used the LoKI simulation tool adopting various formulations to calculate the charged-particle loss frequency due to transport ν transp , namely • Ambipolar-based transport models (see section 3): classical ambipolar diffusion, for EP/HP plasmas; effective diffusion, for multiple positive ions and a single negative ion with low density; eigenvalue diffusion length, describing the radial transport in moderate EN plasmas with multiple positive ions and a single negative ion; QGM, for multiple positive and negative ions; • h-factor transport models (see section 4): LP h-factor, in radial direction joining the LP parabolic-profile for a EN core with EP edges, to the HP flat-topped EN profile, and in axial direction considering the corrected h L -formula proposed in [36] at LP for the various EP/EN plasma-density regimes; HP/EN h-factor, preferred at HP and for high EN plasmas, as suggested in [33,34]; analytical EP h-factor describing the one-dimensional axial transport in an EP plasma with a single positive-ion.
These models were applied to the description of (i) a DC discharge in oxygen, at low to intermediate pressures (p ∼ 25-1300 Pa), adopting a validated kinetic scheme as reference model [42] (see section 5.2); (ii) a microwave discharge in helium, from low to atmospheric pressures (p ∼ 5 − 10 5 Pa), proposing and validating an updated kinetic scheme that bridges two previous models [46,47], adjusted separately for low to intermediate pressures and for atmospheric pressure (see section 5.3). The updated scheme considers the full kinetic description of states He(7 2s+1 l); extends/revises several mechanisms: Penning ionization, dissociative recombination, three-body recombination, and especially AI for which coherent temperature-dependent rate coefficients are proposed; and considers the effective diffusion model to describe the transport of charged particles. Figures 17(a), (b) and 18(a), (b) summarize the simulation results obtained in oxygen at I dc = 10, 40 mA, respectively, for the discharge characteristics of E/N vs. NΛ and the densities of the main plasma species O 2 (a 1 ∆), O 2 (b 3 Σ + g ) and O( 3 P) as a function of NΛ. For helium, similar representations are presented in figures 19(a) and (b) showing, respectively, E rms /N vs. NΛ and the densities of metastables He(2 3 S) and He(2 1 S) as a function of NΛ. In these figures, the results obtained with different charged-particle transport models are depicted by shaded bands, which reveal the impact of these models in the simulations. In oxygen, we observe uncertainties of 20%-60% in the discharge characteristics, which can  [42,43]; bands: regions of values calculated adopting the various charged-particle transport models considered in this work. also reach 60% in the species densities, with larger dispersion at LP and low discharge current. In helium, the uncertainties are in the ranges 8%-115% and 50%-150% for the discharge characteristics and the densities of metastable species, respectively (again, larger dispersion is observed at LP). Note that the modeling results may scatter further, by changing the values of the elementary data (collision cross sections, mean-free-paths, ion mobilities) used to calculate the transport loss-frequencies.
In any case, the simulation results reproduce fairly well the available measurements of the discharge characteristics and the densities of species across the large pressure domain considered, which confirms the predictive power of global models.
We have also shown that, in oxygen and helium plasmas, transport accounts for more than 50% of total charge losses (usually remaining above 80% in oxygen). This result holds for other plasmas, where transport often corresponds to the , O 2 (a 1 ∆) (medium-grey) and O 2 (b 3 Σ + g ) (dark-grey), calculated adopting the various charged-particle transport models considered in this work. dominant loss-mechanism for charged particles, especially at low to intermediate pressures. Naturally, the specific influence of charged-particle transport on the description of a plasma depends on the working conditions (e.g. pressure and dimensions) and on the kinetic features of each plasma species, but in general transport phenomena should be considered and assessed in plasma chemistry studies like any other kinetic mechanism. For example, charged-particle transport should be included as part of the quantitative sensitivity analysis of a kinetic scheme, possibly adopting several models with different validity domains, according to the EP/EN single-/multi-ion plasma under study and the LP/HP conditions considered. Following these practices is essential if the study spans different pressures, as illustrated here when analyzing the analytical EP h-factor (single-ion) model in helium (see figure 16).
To this date, there is no general formulation describing the transport of charged particles in plasmas that can be used in  [46]; bands: regions of values for Erms/N (in (a)), He(2 3 S) (dark-grey in (b)) and He(2 1 S) (light-grey in (b)), calculated adopting the various charged-particle transport models considered in this work.
global models for all types of plasmas and working conditions. Most likely, such a general formulation is out of reach, considering the variety of specific features found in LTPs. Instead, efforts should focus on translating accurate fluid/PIC-MCC simulation results into fitting/abacus formula, as illustrated in the examples of formulations presented in sections 3 and 4, and on revising some of the formulations to express the main kinetic coefficients (mobilities, diffusion coefficients, loss frequencies) as a function of measured/calculated elementary data (ion-scattering cross sections).

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: LXCat (http://www. lxcat.net/).

Appendix. AI rate coefficients in helium
Temperature-dependent rate coefficients for AI processes, He(1 1 S) + He(n 2s+1 l) → e + He + 2 , can be evaluated from kinetic information by proper integration of the corresponding cross section, σ n 2s+1 l (u). However, due to the lack of cross section data, this is only done for triplet states He(3 3 S,P,D) and He(4 3 S,P,D) with the data obtained from [49]. For all other states, the rate coefficients are estimated from those of He(3 3 S,P,D) and He(4 3 S,P,D) as follows: (i) For the singlet states He(3 1 S,P,D) and He(4 1 S,P,D), by scaling the cross sections of the corresponding triplet states considering the ratios of the singlet-to-triplet rate coefficients estimated in [50] from measurements reported in the literature at T g = 300 K (or 2450 K); (ii) For the states He(4 1,3 F), by considering the cross sections of He(4 3 D) scaled with the ratios of the rate coefficients for the states He(4 1,3 F) to He(4 3 D) proposed in [50]; (iii) For the higher levels He(n 2s+1 l) with n = 5, 6, by scaling the cross sections for the corresponding He(4 2s+1 l) states, considering the ratios of the rate coefficients for the levels He(n 2s+1 l) to He(4 2s+1 l) proposed in [47]; (iv) And for the states He(7 2s+1 l), by taking the same cross sections as for the corresponding states He(6 2s+1 l).
In general, for an AI process involving the state He(n 2s+1 l), the rate coefficient is calculated assuming the temperature dependence of another state He(n ′2s ′ +1 l ′ ) and scaling the absolute value using a ratio of rate coefficients at a certain temperature T g 0 : with f Mxw corresponding to a Maxwellian distribution function at temperature T g and Θ n 2s+1 l n ′2s ′ +1 l ′ representing the scaling factor given by C n 2s+1 l σ n ′2s ′ +1 l ′ Θ n 2s+1 l n ′2s ′ +1 l ′ C 3 3 S σ 3 3 S 1.00 × 10 0 C 3 1 S σ 3 3 S 1.00 × 10 1 C 3 3 P σ 3 3 P 1.00 × 10 0 C 3 3 D σ 3 3 D 1.00 × 10 0 C 3 1 D σ 3 3 D 6.25 × 10 0 C 3 1 P σ 3 3 P 2.51 × 10 −1 C 4 3 S σ 4 3 S 1.00 × 10 0 C 4 1 S σ 4 3 S 9.00 × 10 −1 C 4 3 P σ 4 3 P 1.00 × 10 0 C 4 3  Note that the scaling of cross sections considering the ratios of the rate coefficients estimated in [46] or in [50] yields similar results. Table A1 summarizes the scale factors adopted, and table A2 presents the cross section data obtained graphically from [49]. All rate coefficients for AI were obtained by numerical integration of (A.1) with the data in tables A1 and A2, considering an appropriated energy grid and linear extrapolation of the data in table A2. Cross sections for AI with He triplet states He (3 3 S,P,D) and He(4 3 S,P,D). Data obtained from [49]