Bridging the gap between global models and full fluid models: a fast 1D semi-analytical fluid model for electronegative plasmas

Analytical and numerical models allow investigation of complicated discharge phenomena and the interplay that makes plasmas such a complex environment. Global models are quick to implement and can have almost negligible computation cost, but provide only bulk or spatially averaged values. Full fluid models take longer to develop, and can take days to solve, but provide accurate spatio-temporal profiles of the whole plasma. The work presented here details a different type of model, analytically similar to fluid models, but computationally closer to a global model, and able to give spatially resolved solutions for the challenging environment of electronegative plasmas. Included are non-isothermal electrons, gas heating, and coupled neutral dynamics. Solutions are reached in seconds to minutes, and spatial profiles are given for densities, fluxes, and temperatures. This allows the semi-analytical model to fill the gap that exists between global and full fluid models, extending the tools available to researchers. The semi-analytical model can perform broad parameter sweeps that are not practical with more computationally expensive models, as well as exposing non-trivial trends that global models cannot capture. Examples are given for a low pressure oxygen CCP. Excellent agreement is shown with a full fluid model, and comparisons are drawn with the corresponding global model.


Introduction
Numerical and analytical models are invaluable tools to the plasma physics community. Detailed models are able to capture a large range of phenomena, and can provide information on quantities not accessible through experiment [1]. However, even though the capabilities of computer hardware and software have increased dramatically over the last few decades, creating and running a highly detailed model still poses a substantial challenge. Thus the community has developed models that use a variety of approximations, allowing one to reduce the computational complexity of modelling at the expense of some level of accuracy in the results achieved.
A common approach is to average the ensemble of particles over their thermal motion, resulting in fluid-like equations for the density, flux, temperature etc of each species, which can then be solved numerically. Detailed fluid models generally allow for variation in time as well as at least one spatial dimension [1][2][3][4]. Typically the equations are solved using a form of finite element analysis, which allows for simple boundary conditions and solutions to capture local effects. Results obtained from these models typically show good agreement with experiment for the limit of medium to high collisionality [1,4]. However, the computation time required to obtain these results can range from hours to days, depending on the system and the techniques employed.
In order to combat this long time to reach solutions, global models have been developed [5][6][7][8][9], which allow for the solution of bulk properties through a collection of approximations and empirical relations, including the neglection of all spatial derivatives. This allows either rapid convergence to an equilibrium [5] or a description of time evolution of bulk properties [6]. Despite the large number of assumptions made, these models can provide reasonable estimations of bulk values and system trends within a certain parameter space [10]. They are commonly used for systems containing complex chemistries, as their rapid solution time allows the inclusion of many different species and reaction pathways, which would be computationally infeasible with full fluid models [11,12].
Unfortunately, due to the large number of assumptions made, and despite their widespread use, global models struggle to provide good results for systems with high degrees of spatial non-uniformity, or atypical discharges where the important empirical relations break down. Their lack of spatial resolution means that they cannot accurately describe systems where a large fraction of the plasma bulk is not uniform, a common occurrence as spatial gradients often exist long before the development of a sheath. The empirical relations used to link bulk values with sheath edge properties fail to account for the different non-linear couplings that can occur for example between non-uniform densities and power deposition profiles.
The work presented here aims to improve the options available to researchers by bridging the gap that exists between global models and full fluid models. Computational effort is exchanged for analytical intricacy, and differential fluid equations are derived that can be solved in one spatial dimension through an initial value type scheme, thus avoiding the high computation times associated with finite element methods. Example results are given, and compared with a full fluid model [13]. Despite the stricter assumptions, and thus potential innaccuracies as well as a smaller region of applicability, the semi-analtical model is shown to agree well with time averaged results from the more detailed model. In order to quantify any improvement over a global model, one is also created, and the three models are compared together to show their agreements and limitations.

Analytical derivation
The model is based on an idealised radio frequency, capacitively coupled, infinitely planar discharge, as shown in figure 1. As such, spatial derivatives are expressed only perpendicular to the electrodes, and profiles are symmetric about the central plane. The plasma consists of four species: electrons (e), positive ions (i), negative ions (n), and neutral gas (g).

Approximations
Further to the physical characteristics described above, the following approximations are made.
• The described plasma is quasineutral throughout.
• No external magnetic field is applied.
• The system is taken to be in equilibrium, so that time derivatives are removed from equations. • Viscosity effects are negligible.
• Reaction rates can be expressed as the product of a species density and a reaction rate coefficient. • Reaction rate coefficients can be formulated as functions of mean electron energy, taken from an EEDF as calculated by a two term Boltzmann solver [14]. • With the exception of electron reaction rate and electron transport coefficients, all fluid values can be sufficiently approximated by assuming particles have Maxwellian energy distributions. • Elastic collisions between charged species are considered negligible, and self elastic collisions have no effect on a one dimensional system. • The effect of inelastic collisions on directed velocity is small compared with that of elastic collisions. • Thermal gradients do not affect transport rates for species that are not electrons. • Neutral heat flux can be described by Fourier's Law of thermal conductivity. • Electron heat flux can be described by considering only collisions with the neutral gas. • All power deposition into the plasma is in the form of ohmic (collisional) heating of the electrons. • Ohmic power deposition can be calculated from the total current density and the plasma conductivity. • Positive ion velocity is monotonic throughout the discharge.
These characteristics limit the validity of the model to symmetric capacitively coupled systems, operating in a collisional pressure regime, but with a low enough pressure such that three body collisions are not important. It is also limited to systems where stochastic heating is negligible, as is secondary electron emission, so that the plasma is operating in 'α-mode'. Due to time averaging and the enforcement of quasineutrality, the model is precluded from resolving either temporal instabilities, or irregular behaviours of the electric potential, such as double layers. In terms of physical properties, this means that a modelled discharge should have a pressure-length product of roughly 0.1 Pa m-5 Pa m, depending on plasma density, and an ionisation fraction of less than 1%.

Boltzmann moments
The model is constructed from the zeroth to second moments of the Boltzmann equation, given in general form in equations (1)-(3) [15][16][17] with the listed approximations applied, representing conservation of particles, momentum, and energy respectively. In these and subsequent equations, α n , α u , α T , α m , and α Z refer to the density, flow velocity, temperature (in Kelvin), mass, and charge (in units of e), respectively, of species α being acted on by an electric field E. Each species also has a scalar pressure, heat flux, kinetic energy, and volumetric power deposition denoted by α p , α q , α K , and α S respectively. ∇ refers to a spatial gradient, and / δ δt indicates change of a quantity over time due to collisions.

Collision terms
The collision terms ( / δ δt) can be derived by considering the interaction of two particle species, and the differences between their distribution functions before and after the collision. By then taking velocity moments, in the same vein as for the Boltzmann equation, one arrives at a set of fluid collisional terms, given in equations (4)-(8) [15,17]. The change in density (4) is simply defined by the difference between particle gains and losses. For the first moment, velocity changes (5) are reformed as the difference between momentum ( ) α α u n and density changes [17]. The momentum changes (6) are given as the common Krook collision operator plus a consideration for the change in elastic collision rate over energy gradients, α g . For the second moment, three effects are accounted for when considering the collisional change in mean particle energy (7): kinetic energy changes due to elastic collisions, momentum changes, and changes brought about through the creation and loss of particles. For changes in kinetic energy due to elastic collisions, considerations are made for the change of both flow and random velocities of particles, leading to two terms in (8) describing both thermal and kinetic energy transfer between species. In (4)-(8) and following expressions, the summation over α R corresponds to the collection of all particle creation and destruction mechanisms for species α, where α G R is the number of species α created by reaction R, n R1 and n R2 are the reaction partners, and K R is the reaction rate coefficient. Summation over β refers to the sum over the given expression for all of the elastic collision partners β of species α, ν = αβ β αβ n K is the momentum transfer collision rate and /( ) κ = + αβ α β α β m m m m 2 2 is the kinetic energy transfer coefficient between species α and β.

Closure terms
As stated in section 2.1, the heat flux for the neutrals and electrons is estimated via Fourier's Law (9) and a weakly ionised approximation (10), respectively. In these equations, is the thermal conductivity of the neutral gas, where the assumption is made of a linear dependence on gas temperature, with the coefficients h a and h b . The fractional change in collision frequency over electron temperature gradients is represented by 2 e e e eg e e e B e e (10)

Normalisation scheme
The choice of normalisation scheme is highly influential on the functionality and capabilities of the resulting model. The scheme used, given below, builds on the work of Raimbault and Liard [18,19] and extends to an electronegative plasma with electron temperature gradients. In this scheme: n f is the gas fill density; 1 2 is the positive ion Bohm velocity; K 0 is a normalisation reaction rate coefficient, in this case = K K 0 i g , the momentum transfer collision rate coefficient between positive ions and neutrals; and T e0 is the central electron temperature.  A normalised electric field, ζ, falls naturally from this normalisation scheme, as given in (11). The normalisations chosen lead to a coupling of the pressure and system length, as is known to happen in physical systems, such that the pressure-length product of the model can be extracted independently of the pressure, as shown in (12), where L is the normalised system length. This value is useful for comparison with physical systems, but is not required for the operation of the model. One can also change between a power density in Wm −3 to the normalised value, Σ, as well as a heat flux q to the normalised value of Q . These derived normalisations are defined in (11)- (14).

General equations
By combining (1)-(10) and normalising as per the given scheme, one obtains a series of general expressions (15)- (18), which describe normalised density and flux gradients for all species, and the second differential of the temperature coefficient for neutrals and electrons. For more detail see appendix A.
Equation (15) comes directly from (1) and (4). The density gradient, (16), stems from (2) and the expansion of its col lision term using (4)- (6). The second differential of the temper ature coefficient, ″ γ , is a result of having to differentiate (9) and (10), before insertion into (3) along with the collection of collision terms (4)- (8). Equations (15) and (16) apply to all species, while in this system (17) applies only to the neutral species, and electron temperature gradients are described by (18). Note that the elastic collision energy dependency factors α g are assumed to be non-zero only for electrons.
As detailed in appendix A, there are a large number of cancellations, combinations, and simplifications that occur on the route from the Boltzmann moments to (15)- (18), and so precise physical interpretations for each term are not necessarily straightforward. However, it is still useful to understand the effects that are captured in each equation, and so a description of each term is given below.
In (15) the understanding is simple, as the gradient of the flux of each species is simply described by the net creation or loss of particles that occurs at a particular point.
The terms in (16) are a little more complex, as some have been combined. The first term on the right hand side is simple, and is the effect of the time averaged electric field on the species α if it carries a charge. The second term is used only for species with a temperature gradient, and gives the effect on species densities due to relative changes in temperature, including gradients created by a non-uniform elastic collision frequency, if applicable. Term three describes the changes in density due to the creation and loss of particles, equated to the flux gradient for brevity. The final term describes how elastic collisions affects the density of species α due to the resulting changes in momentum.
Both (17) and (18) detail the spatial change in the temperature gradient due to increases and decreases of the species energy density from various sources. The first term of (17) and the first two of (18) describe energy changes due to the gradient of the heat flux, calculated in (17) using Fourier's Law of thermal conductivity (9) and in (18) from the expression given by (10). Term two in (17) and term three in (18) give the changes in thermal energy density due to the changing particle density. The next term gives changes due to the flow of particles, including the effects of a non-uniform elastic col lision frequency in (18). Term 4 in (17) and term five in (18) describe how both the random and directed kinetic energies of particles of type α are changed by elastic collisions. The final term in (17) and the sixth in (18) comes from a combination of changes in energy density due to the creation and loss of particles, as well as additional energy changes due to the flow of particles. Term seven in (18) gives the changes in energy from external sources and sinks, including the ohmic power deposition and inelastic collisional losses (to be discussed).

Electric fields and power deposition
As the system is assumed to be quasineutral throughout, it is approximated that the use of Poisson's equation is not necessary. This is particularly advantageous as it is known to negatively affect the numerical complexity of fluid models by introducing a high level of stiffness into the equation set, which requires more complicated numerical solvers and smaller step sizes than a non-stiff equation set. However, because of this, a different term must be used for the time averaged electric field found in (16).
In this time averaged, symmetric, and fully quasineutral system, the electric field ζ can be obtained from a rearrangement of (16) and the derivative of the quasineutrality condition, ∑ = ′ α α α Z N 0. For brevity, terms in (16) that are independent of the electric field are collected into α C , so that one obtains From this it is straightforward to combine and rearrange these two expressions to solve for ζ. The resulting 'equilibrium' electric field term, given in (19), is similar in nature to an ambipolar field. It can be explained as the electric field necessary to counteract all of the forces acting on the charged species, in order to maintain quasineutrality.
In order to more closely represent a physical system, the model needs to take into account the non-uniform characteristics of the ohmic power deposition to the electrons, which is not possible in a global model. As stated in section 2.1, the assumption is made that the power can be calculated from the electron current and the plasma conductivity [16]. In order to find the time averaged power deposition, one must take the amplitude of the sinusoidally varying current density, given by (20). This is due to ohmic heating in an RF plasma being an intrinsically time dependent phenomenon, and so using only time averaged values in its calculation would result in an incorrect value. The time average power deposition is therefore given by (21) and can be normalised using (13) to obtain (22).

Computational considerations
By recasting each of the second order differential equations for temperature gradients as a pair of first order differential equations, one obtains a system of twelve strongly coupled, non-linear, ordinary, differential equations. This set of expressions fully describes the gradients of fluxes, densities, and temper atures, and is suitable for numerical integration as either a boundary value problem (BVP) or initial value problem (IVP). The system is assumed to be symmetrical, so that boundary conditions are to be given at the centre (X = 0) and sheath edge (X = L) of the plasma. The central conditions are largely controlled by the requirement for symmetry, and the edge conditions by contact with the sheath. These are summarised in table 1.
Due to the combination of boundary conditions, the only unknowns to be used as inputs to the model are the electron density,central electron temperature, and the normalised cur rent density. All other quantities are defined by one or more boundary conditions. The use of conservation laws, such as quasineutrality and current conservation, allows further reduction of the free parameters by, for example, introducing the electronegativity, α, as a control parameter for the ion densities. The model is converse to other models, part icularly global models, in that one specifies the plasma parameters such as electron density, and obtains system properties such as pres sure-length product. This does not preclude the model from describing the same systems as others, but simply requires a different method of thinking. Note that the inclusion of variable central values for neutral gas density and temperature, and their corresponding spatial evolution, allows for the inclusion of neutral gas depletion effects, known to be important particularly in high density discharges [19][20][21][22].

Numerical considerations
As stated above, the system of equations is suitable for integration as either a BVP or IVP. Formulation as a BVP would allow a reasonably simple solution through discretisation, but this faces the same long execution time as a full fluid model. Solving the system as an IVP gives the potential for a greatly reduced integration time, at the expense of a more complicated numerical scheme. As one of the main motivations for this model is fast computation, the ability to solve as an IVP To be calculated based on edge constraints. is valuable. This section details the numerical considerations and algorithms required to solve the system.

Numerical integration
The boundary conditions specified in table 1 indicate that the integration as an IVP is best performed spatially from the centre of the discharge to the edge. This is performed with the ode113 numerical integration routine in MATLAB 7.14 [23], which uses a predictor-corrector, linear, variable order, multistep solver (Adams-Bashforth-Moulton method [24]). An unfortunate side effect of the normalisation scheme is that, due to the decoupling of the physical discharge parameters, the spatial extent of the plasma is not known until the edge boundary conditions are met. The normalised positive ion velocity, / = Γ V N i i i , is monotonic over the discharge, and has defined central and edge values for all discharge parameters, so can be used as an integration coordinate. As shown in table 1, the sheath boundary is determined by the point at which the positive ions reach the Bohm velocity. As velocities are normalised to this value, the integration bounds are determined by = V 0 i at the centre, and = V 1 i at the edge. This change is effected by simply dividing the calculated gradients by ′ V i and including an extra variable to track the true spatial coordinate. It is worth noting that despite there existing a modified Bohm criterion for electronegative plasmas [25], it is derived by assuming that both negative ions and electrons are in Boltzmann equilibrium with the plasma potential, and that negative ions are present at the sheath edge. In the semianalytical model, it is assumed that there is no negative ion production in the sheath. Thus the flux, and so too the density, of negative ions at the sheath boundary must be zero. In such a situation the modified Bohm criterion reverts to that of the electropositive case, so this value is used to estimate the position of the sheath edge.
A further change to the integration scheme is made by transforming the equations to describe the natural logarithm of density values. This is done to improve the numerical stability of the integration due to the occasionally large differences between density values, and also prevents overshoot to negative density if the automatic step size is too large. To make this change, (16) is divided by the species density to provide the logarithmic derivative, and the other equations are updated to accept ( ) α N ln as arguments.

Electronegativity minimisation
As shown in table 1, the negative ion flux ( ) Γ n is required to be zero at the system edge, that is Γ = 0 n L , , where L refers to edge values. As symmetry dictates also that Γ = 0 n,0 , it is the species densities that have control over whether or not this edge boundary condition is met. Specifically, Γ n L , is controlled indirectly by the central electronegativity, α 0 , which determines the central ion densities for a given central electron density. It is therefore necessary to repeat the integration with different values of α 0 and minimise the value of Γ n L , iteratively. Unfortunately, due to the highly non-linear characteristics of the equation system, the parameter space of α 0 and Γ n L , is not trivial, and contains steep gradients and discontinuities. Thus standard minimisation routines often fail to converge, or are impractically slow. In order to have an automatic solution, a custom minimisation algorithm has been created, using minimisation by bisection, with a variety of integration outputs being used to indicate in which direction the minimisation should progress. This proves to be a robust method for a large range of input conditions.

Perturbations
Despite the minimisation routine being able to give the required α 0 to the limit of double precision (approximately 15 significant figures), there are still cases where the highly non-linear nature of the equations prevents the condition Γ = 0 n L , from being met to an acceptable level, through an inability to specify a precise enough α 0 . In order to access the trajectories that meet the boundary conditions, numerical perturbations are applied partway through the integration. An example of this process is shown in figure 2. Figure 2(a) shows the original integration output with α 0 specified to the limit of double precision. A perturbation is applied to the positive and negative ion densities of order / ∆ ≈ − N N 10 7 partway along the trajectory and the integration continued from this point, which then meets the boundary conditions, as shown in figure 2(b). The position and magnitude of the perturbation are found through bisection. The magnitude of the perturbations is significantly smaller than the relative derivatives of the ion densities. Values between the beginning of the integration and the point of perturbation are therefore still accurate to the numerical precision of the system.

Neutral properties minimisation
Similarly to the negative ion flux, the neutral density and temperature are also specified at the edge of the integration, but unlike Γ n , they are free parameters at the centre. As Γ n L , is largely independent of both N g, 0 and T g, 0 , with the exception of the weak effect of elastic collisions, the neutral parameters can be solved separately to α 0 . However as the density and temperature are closely coupled, it is logical to minimise for them simultaneously. Once again, however, the shared para meter space of N g, L and T g, L is non-trivial for the inputs N g, 0 and T g, 0 , and contains discontinuities and inaccessible regions. A custom minimisation routine has been created, using a combination of simple independent linear extrapolation, 2D co-dependent linear regression, and Monte Carlo techniques, depending on the current knowledge of the parameter space. A record is kept of each point in parameter space that is tested, so that the region bounding the correct solution can be found. Extrapolation is used until enough points have been tested to provide a suitable data set for the regression. If the regression fails, then a Monte Carlo technique is used to find more suitable points in the parameter space.

Results
In order to present results of the semi-analytical model, oxygen has been used as an example of an electronegative gas, although application to other gases is possible. For different gases the equation set will be similar, but rate coefficients and therefore results will differ, as well as there being the possibility of different numerical behaviours. In order to keep computation time as short as possible, a reduced reaction set for oxygen was chosen, that also enables comparison with a time and space resolved fluid model of oxygen [13]. This set, detailed in table 2, consists of eight particle creation/ destruction mechanisms, plus two that are considered only for electron energy loss, as well as elastic (momentum transfer) collisions between charged particles and neutrals.
For this reaction set, all reaction rate coefficients that do not have an explicit form can be well approximated, across electron energies of 0-50 eV, to within about 1% by the form of (23), where the coefficients a R, n are found through linear regression in logarithmic space. This form also conveniently provides expressions for g e and / × ∂ ∂ T g T e e e through differentiation of the expression for K eg .  × − 2 10 15 [32] a Recombination to excited state and subsequent de-excitation is considered but not explicitly included. b Density of ( ) ∆ O 2 1 g estimated using empirical fit to T e0 from [33] c Used only in calculating electron energy loss. Note: T g in Kelvin. ( ) f T e indicates RRC estimated from tabulated energy versus cross section data, via a two term Boltzmann approximation [14].
Also required are the thermal conductivity parameters for the neutral energy balance. The values for these are obtained using data provided by the National Institute of Standards and Technology, accessible in an online database [34]. From a linear fit to this data, the coefficients are =

Example outputs
The definition of the reaction scheme completes the picture needed to create a model of an oxygen discharge. Following implementation of all inputs, equations, and numerics in Matlab, one can generate results using just the three inputs previously mentioned: normalised electron density, electron temperature, and the normalised current density. An example set of results is given in figure 3. Shown are species densities, fluxes, and temperatures from the centre to the edge of the discharge, plotted on normalised axes, with the exception of temperatures, which have been denormalised for ease of understanding. Also given is the profile of the reaction rate for ionisation (reaction I in table 1). Figure 3(a) shows the charged species densities, where particularly apparent is the transition between the bulk region and the presheath. This can also be seen in the fluxes, given in figure 3(b), which additionally give a clear visual indication of the current and particle flux conservations that result from the equation set. The steady increase of T e through the discharge is shown in figure 3(c), as is the spatially resolved ionisation rate. These three plots together show how the system behaviour is affected by the sharp dependence of, in particular, K I on T e in this energy range. As T e increases, reaction rate coefficients alter rapidly due to their individual nonlinear trends, and species behaviour can change quickly in space. This is the cause of the sudden transition into a presheath. Close to the edge, the dropping electron density causes a turnover in the ionisation rate, despite the continued increase in T e . In a physical system, T e would reach a peak just inside the sheath, then tend toward zero as one approaches the wall. However the lack of sheath in this model means that if one were to continue the integration beyond the Bohm criterion, T e would keep increasing, as the relationship of deposited power with − N e 1 means both become singular as → N 0 e . This cessation of the numerical integration is also the reason for the non-zero ion and electron density at the integration edge. Figure 3(d) shows how the neutral properties have roughly parabolic profiles of both density and temperature, as their spatial gradients are slow to change. In this example, over 99% of the central density depletion is due to the increase in temperature. In this model, the rest is due to frictional forces from the ions and electrons. The wide range of conditions possible in the semi-analytical model is demonstrated in figure 4, which shows electronegativity results for a parameter sweep across central electron temperatures of 3 eV-4.8 eV (linearly) and relative electron densities of 10 −9 to 10 −4 (logarithmically). To ease comprehension, the parameters of pressure-length scale and current density have been used as axes instead (assuming a 4 cm plasma length for the denormalisation). In order to generate such a broad parameter sweep, a relationship between electron density and power deposition was found such that Σ = N e e . This means that, from (13), on average each electron gains one unit of energy ( ) k T B e0 per time unit ( / ) n K 1 f 0 , and so the normalised power per electron is kept constant. From this, the normalised current density factor in parentheses in (22) is estimated. . Shown are (a) charged particle densities, (b) particle fluxes, (c) electron temperature and ionisation rate profiles, and (d) neutral properties. In (a) and (b): --denotes positive ion properties; ---negative ions; · · · · · · electrons; and -· -neutrals. Transition from bulk to presheath is seen in (a) and accompanied by a peak in ionisation rate in (c). Current and mass conservation is clear in (b). Neutral gas depletion is evidenced in (d).
It is clear that α 0 is largely independent of plasma density at low discharge current densities, but at higher current densities a dependence is seen. This can be explained through the connection of plasma current density with the ionisation fraction; at the lowest current densities in figure 4, the ionisation fraction is of the order of 10 −8 , but is closer to 10 −4 at the highest current densities. At low ionisation fractions, the negative ion destruction pathways involving neutral species dominate many times over those between charged particles only. However, as one increases the current density, and so too the ionisation fraction, reactions between charged species start to become important. This then elicits a mode transition such that destruction of negative ions through collisions with positive ions and electrons starts to noticeably alter the dominant reactions of the system, at an ionisation fraction of around 10 −6 , causing a reduction in electronegativity with increasing current density. Figure 5 gives examples of two systems with different input conditions, showing qualitative differences that can occur between points in parameter space. Figures 5(a)-(c) show a system with a relatively low electron temperature (i.e. higher system pressure). In figure 5(a), ion densities actually increase slightly from the centre until T e builds to a point where reaction rate coefficient ratios cause a shift in density gradients, and a collapse through the presheath towards the sheath. Similar to the example given in figure 3, the ionisation rate, shown in figure 5(c), increases away from the centre, but the starting value is much lower. This is linked to the more uniform behaviour in the bulk when compared with figure 3. At the point at which the plasma transitions into the presheath, the gradient in both T e and ionisation rate is much steeper than that of figure 3, so that the change in T e across the system is much greater and the edge value is comparable with that in figure 3(c). This in contrast to the example given in figures 5(d)-(f), where T e is already high at the centre (thus representing a lower pressure discharge), and ion densities start decreasing immediately from the centre. The ionisation rate is noticeably smoother in figure 5(f) than in figure 5(c), and although there is still a marked transition into the presheath, it is much gentler than in the high pressure example. The shallower gradients mean that T e actually reaches a lower value at the system edge than in figure 5(c).

Comparison
In order to test the model, comparisons were drawn between the one presented here and a time and space resolved fluid model [13], in order to see how the neglection of time dependencies, sheath effects, and wall properties affect the model. It is known that wall interactions play a significant role in oxygen plasmas [35], particularly for the dynamics of ( ) ∆ O 2 1 g , which takes part in a number of the dominant reactions [36]. It is an important reaction partner for the − O negative ion, and an additional ionisation pathway for the creation of + O 2 . Therefore of particular interest is how the empirical inclusion of ( ) ∆ O 2 1 g in the semi-analytic model compares to the self consistent inclusion of ( ) ∆ O 2 1 g in the full fluid model. The central electron temperature, relative electron density, and normalised current density was calculated for time averaged results from the full fluid model and input to the semi-analytical model. The results were then denormalised, and plotted together in figure 6.
As demonstrated, there is excellent quantitative and qualitative agreement between the bulk density profiles of the two models, however the difference in behaviour due to the lack of sheath is clear. The full fluid model can support a deviation from quasineutrality through the solution of Poisson's equation, whereas in the semi-analytical model, the charged particle densities collapse rather than create a net space charge. The full fluid model transitions smoothly into a sheath, and indeed it is difficult to determine where the bulk ends and sheath begins without a rigorous definition [37], whereas the transition between bulk and sheath is significantly sharper in the semi-analytic model.
Further differences can be seen in the position of the transition, as well as in the exact value of the electronegativity. These can both be explained by differences in not only the spatial profiles but also the time dependencies of power deposition and electron temperatures. Non-ohmic heating mechanisms have been shown to effect a modulation of T e in time [38][39][40], which is not captured in the semi-analytical model. This temporal behaviour of T e affects the reaction pathways through the non-linear dependence of reaction rate coefficients. Neither the temporal modulation or non-ohmic heating are included in the semi-analytical model. Additionally, as a side effect of the enforcement of quasineutrality, the integration terminates with a smaller spatial extent than the full fluid model. This is fully expected as no sheath can exist in a quasineutral system. One can estimate the true pressure-length product if one has knowledge of the approximate fractional sheath size of the physical system being modelled. The comparison between the semi-analytical and full fluid model demonstrates that there is a difference in behaviour between the two, as is to be expected, however it is important to know if the considerable increase in computational efficiency is worth the discrepancy in outputs. For this purpose a global model has been created, similar to that developed by Monahan and Turner [41], but adjusted for a radio frequency plasma, and using the same reaction set as given in table 2.
See appendix B for more details on the creation of this model. Results for all three models were calculated for a range of system pressures, to see how the two models developed here changed in comparison to the full fluid model. These are provided in figure 7, which gives detail of the outputs of each model for each set of inputs derived from the full fluid model.
As an indication of the underlying reaction set, figure 7(a) shows the very similar central electronegativity results from  d) and (e), --denotes positive ion properties; ---negative ions; · · · · · · electrons; and -· -neutrals. Differences in behaviour due to T e0 are clear. Lower central T e0 ((a)-(c)) causes low ionisation and more uniform bulk, with sudden turnover into presheath. Higher T e0 ((d)-(f)) results in more bulk ionisation, and a smoother transition toward the sheath.  In (a), semi-analytical model positive and negative ion densities are shown by --and -· -respectively, whereas ---and · · · · · · are the time averaged positive and negative ion densities, respectively, from the full fluid model. In (b) --is the semi-analytical model electron density, and the time averaged value from the full model is represented by ---. Density profiles show excellent qualitative and quantitative agreement between the models, though the transition to sheath is notably sharper in the semi-analytical model. each of the three models. Clear are slight differences in the α 0 of the global and semi-analytical models when compared to the full fluid model, particularly at the lower pressures. The existence of such a discrepancy is not unexpected, and the worse behaviour at lower pressures can be attributed to the increasing sheath fraction seen in the full fluid model with decreasing pressure-length product. The encroachment of the sheath onto the bulk will reduce the effectiveness of the empirical relations of the global model, and increase the impact of effects in the sheath that the semi-analytical model cannot capture.
For the semi-analytical and global models, the pressurelength product reported is that of the bulk (and presheath) only. One also can calculate a time averaged bulk width, and thus a bulk pressure-length product, for the full fluid model [37]. The ratio of this value to the full system pressure-length product, as specified by the full fluid model, is given in figure 7(b). This data shows the expected trend of decreasing sheath width with increasing pressure for the fluid model. For the semi-analytical model, the quantitative agreement is good with the full fluid model, but there is a slight underestimation of the spatial extent of the system as one increases the pressure. This is not the case for the global model, which has a distinct inability to accurately specify the pressure-length product when programmed to provide a specific absolute electron density and electron temperature. This is a result of the changes described in appendix B that are required to give the same electron properties as the full fluid model, particularly the impact of these parameters on the electron power balance, which has the greatest control of the relative electron density. The power balance is strongly affected by approximations such as the edge to centre density ratio as well as the sheath voltage. Figures 7(c) and (d) show the ionisation fraction for each of the three models, defined as the ratio of positive ion to neutral particle densities, as a function of pressure-length product (c) and current density (d). These plots capture the reduced plasma density in the global model when compared to the other two, and demonstrate the close matching of the semianalytical and full fluid models. Without the spatially resolved reaction rates, such as those in figures 3(c), 5(c) and (f), the global model struggles to account for the correct amount of species creation and destruction, resulting in an underestimation of the plasma density.

Discussion
Looking first at figure 6, it is clear that although the qualitative and quantitative agreement between the semi-analytical and full fluid models is very good, there are still notable differences. While the spatially dependent power deposition profile allows, for example, the capture of non-monotonic density profiles, the precise spatial positioning of features differs between the two models. This is due to the neglection of space charges in the semi-analytical model, and thus performance decreases as one progresses toward the sheath, due to the breakdown of quasineutrality. As convenient as it is to imagine a distinct boundary between the sheath and the plasma bulk, this is not the case in a physical system, and so a gradual degradation in performance is to be expected as the space charge density builds. The neglection of the sheath affects the system not only through a change in edge characteristics, but also by preventing the flow of information from any wall interactions back into the bulk, making it impossible to self consistently portray the behaviour of species dominated by wall effects, such as ( ) ∆ O 2 1 g [35]. Despite this, figure 6 shows that the bulk is largely unaffected, particularly the value of the electronegativity which is mostly dependent on the reaction set used by the model. The assumption of quasineutrality also forbids the creation of stratified pre-sheath structures (double layers), known to appear under certain conditions [5,41,42]. Despite these issues, the assumption of quasineutrality, and thus the removal of Poisson's equation from the system, prevents the creation of a stiff set of differential equations, and thus improves both the model simplicity and computation time dramatically. The improvement in computational performance achieved by discarding Poisson's equation is so great that it is deemed a necessary sacrifice to improve the usefulness of the model.
A further known loss of information comes from the time averaging of equations. It is known that electronegative plasmas exhibit temporal instabilities under certain conditions [6,[43][44][45]. The combined loss of stratified pre-sheaths and temporal instabilities may explain some of the region where no solutions are possible, seen in figure 4. As shown in [6,43] the appearance of instabilities occurs as one increases the plasma power toward the transition to γ-mode, which itself cannot be captured due to the lack of wall interactions. [45] also reports that instabilities are more frequent at higher pressures, possibly explaining the shape of the solution boundary shown in figure 4.
The self consistent inclusion of neutral species has a clear effect on discharges with high ionisation fractions ( − 10 6 ), as seen in figure 3, and in figure 4 where some of the dependence of α 0 on current density is due to neutral gas depletion. This shows that in highly ionised plasmas the effect of neutrals cannot be neglected, as there is a feedback from the neutral species onto the ions and electrons, affecting plasma properties.
The comparisons between the three types of model is particularly valuable in evaluating the relative accuracy and usefulness of the different models, although the reversal of inputs and outputs when compared to more conventional global models, as discussed at the end of section 2, does require a different viewpoint. As seen in figure 7, while the identical reaction set and the specification of T e0 and N e dictates that the electronegativity is comparable between the three models, the lack of spatial information in the global model causes large errors in the pressure-length product and the ionisation fraction. The discrepancies between the semi-analytical and full fluid models are within what would be expected due to the lack of consideration of time dependencies and sheaths, and are significantly smaller than those between the global and full fluid models. The inclusion of spatial resolution of ionisation profiles in the semi-analytical model allows for much better estimation of plasma densities than the approximations of a zero dimensional model.

Conclusions
The results presented show that the semi-analytical model detailed here provides data of far greater utility than a typical global model, while maintaining the seconds to minutes execution time, as well as the simple user interface once the model is constructed. The drastically increased underlying complexity comes at a cost of flexibility, with a reduction in input parameter range when compared to a global model. However the remaining parameter range and normalisation scheme mean that the semi-analytical model is still applicable to a large range of systems. The applicability of the semianalytical fluid model is limited mainly by the specification of symmetry and a planar geometry; constraints that are not even considered in the global model. In comparison to the full fluid model, the applicability of the semi-analytical model is reduced due to the lack of temporal information. However within the available range of parameter space, the results for the time averaged plasma bulk are quantitatively very similar, but the time cost to the user is hundreds of times smaller for the semi-analytical model. The potential applications for this type of semi-analytical model are broad, and would be particularly useful for fast characterisation of systems with a high degree of non-uniformity.
Potential future improvements to this model include the introduction of a method for treatment of wall-dominated species that does not rely on empirical relations, and the removal of numerical singularities that may allow solutions to extend to the point at which charged species densities become zero.
The transformation of the general Boltzmann moments into normalised equations suitable for numerical integration is not always trivial, particularly for the higher moments. This section details the steps required between (1)- (3) and (15)- (18).
The transformation of the flux gradient into the normalised form is the least complicated, and simply involves the combination and normalisation of the Boltzmann moment and collision term. From (1) and (4): The importance of the normalisation scheme is particularly clear in (16), where the form of u B is particularly useful. From (2) and (4) The creation of the equations for the temperature coefficient are the most complex in the model, and are best tackled in discrete sections of Boltzmann moment, collision terms, and the two heat flux gradients.
From (3): From (4)- (8): From the neutral heat flux given by (9): From the electron heat flux given by (10): With the assumptions that   For the neutral species, assembling (A.9), (A.14), and (A.17) leads, after a small amount of manipulation, to (17). A slightly more in depth series of rearrangements, collections, and cancellations leads from (A.9), (A.14) and (A.20) to (18). The chosen normalisation scheme means that there are no residual constants in the final equations, and one is left only with the relationships between species.

Appendix B. Global model
A global model of an oxygen plasma consisting of one positive ion species, one negative ion species, and electrons, has been developed based on the work of Monahan and Turner [41], which in turn builds on the work of Kim et al [5,10] and Monahan [46]. Changes have been made to allow for a different reaction set and reactor parameters.
Global model derivation begins with equations for the particle (B.1) and energy (B.2) balances, which are derived from the zeroth and second moments of the Boltzmann equation, respectively, through the removal of spatial gradients and terms involving fluid flow. In these expressions ( ) α α u n L is the flux of species α at the wall, A is the reactor wall surface area, and V is the reactor vessel volume. T e K represents electron temperature measured in degrees Kelvin, P abs is the absolute power deposited into the plasma, P colls is the total power loss due to collisions, and E ei is the energy lost per electron-ion pair crossing the sheath. Dividing by reactor volume allows the particle balance to be expressed as: In (B.5) the first term on the right hand side represents elastic collisional energy losses, the second inelastic collisional energy losses, and the third temperature changes due to the creation and destruction of electrons.
The ion flux to the wall is estimated as the Bohm velocity multiplied by the ion sheath edge density = n h n i s l i , , 0 , where h l is the centre to edge density ratio, a recurring parameter found across the literature. For an electronegative plasma [41] uses an ansatz comprised of three empirical relations for different pressure regimes, as proposed by [5]. These empirical relations have been created to fit a limited data set from more detailed models. For the assumption of equal positive and negative ion temperatures ( ) ≈ = + − T T T i , the three components are as given below: In (B.6)-(B.9), l p is the discharge length, λ i is the positive ion mean free path, and K MN is the positive ion-negative ion recombination rate coefficient. These three components are added in quadrature to give an estimate of h l such that:  The energy lost per ion-electron pair requires an estimation of the sheath voltage, V s . For an RF CCP, this differs from that used by [41], and expressions are taken from Lieberman and Lichtenberg [16] and rearranged to find V s as a function of input power and frequency as given below in (B.12), where ω is the driving frequency in Rad s −1 and ν = n K m e, g eg is the electron elastic collision rate. For the limit of an infinite planar discharge, the area to volume ratio of the reactor simplifies to 2/l p . This with (B.5)-(B.12), leave the input power density, the plasma length, and system pressure (via neutral particle density) as input parameters to the model. The three density equations from (B.3) are combined with quasineutrality to leave one density equation and one equation describing the time evolution of electronegativity.
The system of three differential equations (electron density, electronegativity, and electron temperature) are evolved numerically until an approximate steady state is reached, though this is frequently subject to small scale oscillations. For the operating regime of interest, it is found that the equilibrium point reached is independent of the initial conditions used. Nevertheless, the starting conditions of the global model are chosen to be the same as the system of interest, for ease of comparison. As this model uses different inputs and generates different outputs to the semi-analytical model (electron temperature and density and outputs as opposed to inputs), it was run inside of a minimisation routine that allowed one to specify the desired electron properties and find the plasma length and system pressure that provide these. Otherwise a direct comparison between the global and semi-analytic models would be difficult.